# Continuously Stirred Tank Reactor (TACO)

Jump to: navigation, search

This page contains a model of the Continuously Stirred Tank Reactor problem in AMPL format, making use of the TACO toolkit for AMPL control optimization extensions. The original model can e.g. be found in [Diehl2001]Author: M. Diehl
School: Universit\"at Heidelberg
Title: Real-Time Optimization for Large Scale Nonlinear Processes
Url: http://www.ub.uni-heidelberg.de/archiv/1659/
Year: 2001 . Note that you will need to include a generic AMPL/TACO support file, OptimalControl.mod. To solve this model, you require an optimal control or NLP code that uses the TACO toolkit to support the AMPL optimal control extensions.

### AMPL

This is the source file cstr_taco.mod

```# ----------------------------------------------------------------
# Continuously stirred tank reactor using AMPL and TACO
# (c) Christian Kirches, Sven Leyffer
#
# Source: M. Diehl, 2001
# ----------------------------------------------------------------
include OptimalControl.mod;

# time and fixed end time [in s]
var t;
var tf := 1500.0, >= 1500, <= 1500;

# molar concentration of species A [mol/l]
var cA := 2.14, >= -0.02, <= 6.0;

# molar concentration of species B [mol/l]
var cB := 1.09, >= -0.02, <= 4.0;

# reactor temperature [degrees Celsius]
var Tr := 114.2, >= 50.0, <= 160.0;

# jacket temperature [degrees Celsius]
var Tj := 112.9, >= 50.0, <= 160.0;

# feed flow control [1/h]
var ff := 14.19, >= 3.0, <= 35.0;
let ff.type := "u0";

# cooling rate control [kJ/h]
var cr := -1113.5, >= -9000.0, <= 0.0;
let cr.type := "u0";

# pre-set deviation of initial state from steady state
# set this to something different from 1.0
param alpha := 0.8;

# various chemistry parameters
param k10 :=      1.287E+12;
param k20 :=      1.287E+12;
param k30 :=      9.043E+09;
param cA0 :=      5.1;
param E1  :=  -9758.3;
param E2  :=  -9758.3;
param E3  :=  -8560.0;
param theta0 := 104.9;
param rho :=      0.9342;
param Cp  :=      3.01;
param H1  :=      4.2;
param H2  :=    -11.0;
param H3  :=    -41.85;
param kw  :=   4032.0;
param AR  :=      0.215;
param VR  :=     10.0;
param mK  :=      5.0;
param CPK :=      2.0;

# steady state values
param cAs      :=     2.14021053017;
param cBs      :=     1.09030436131;
param theta_s  :=     1.14191084421E+02;
param theta_Ks :=     1.12906592910E+02;
param FFs      :=    14.19;
param Qdot_Ks  := -1113.5;

# objective weights
param Q11 := sqrt (0.2);
param Q22 := sqrt (1.0);
param Q33 := sqrt (0.5);
param Q44 := sqrt (0.2);
param R11 := sqrt (0.5);
param R22 := sqrt (5.0E-7);

# least-squares deviation from steady-state
minimize Deviation:
integral (   ((cA - cAs)*Q11)^2
+ ((cB - cBs)*Q22)^2
+ ((Tr - theta_s)*Q33)^2
+ ((Tj - theta_Ks)*Q44)^2
+ ((ff - FFs)*R11)^2
+ ((cr - Qdot_Ks)*R22)^2, tf );
let Deviation.scale := 100.0;

var k1 = (k10*exp(E1/(273.15+Tr)));
var k2 = (k20*exp(E2/(273.15+Tr)));
var k3 = (k30*exp(E3/(273.15+Tr)));

# dynamics
subject to

ODE_cA: diff(cA,t) = (1.0/3600.0) * (ff*(cA0 - cA) - k1*cA - k3*cA^2);

ODE_cB: diff(cB,t) = (1.0/3600.0) * (- ff*cB + k1*cA - k2*cB);

ODE_Tr: diff(Tr,t) = (1.0/3600.0) * ( ff * (theta0 - Tr)
- 1.0/(rho*Cp) * (k1*cA*H1 + k2*cB*H2 + k3*cA^2*H3)
+ kw*AR / (rho*Cp*VR) * (Tj - Tr) );
ODE_Tj: diff(Tj,t) = (1.0/3600.0) * ( 1.0/(mK*CPK) * (cr + kw*AR * (Tr - Tj)) );

# initial value constraint
IVC_cA: eval(cA,0) =              (1-alpha)*cAs;
let IVC_cA.type := "dpc";

IVC_cB: eval(cB,0) =              (1-alpha)*cBs;
let IVC_cB.type := "dpc";

IVC_Tr: eval(Tr,0) = alpha*85.0 + (1-alpha)*theta_s;
let IVC_Tr.type := "dpc";

IVC_Tj: eval(Tj,0) = alpha*85.0 + (1-alpha)*theta_Ks;
let IVC_Tj.type := "dpc";

option solver ...;

solve;```