Semibatch esterification problem (TACO)

This page contains a model of the Semibatch esterification problem in AMPL format, making use of the TACO toolkit for AMPL control optimization extensions. The original model is due to [Kuehl2005]Author: P. K\"uhl; A. Milewska; M. Diehl; E. Molga; H.G. Bock
Booktitle: Proc. Int. Workshop on Assessment and Future Directions of NMPC
Pages: 467--474
Title: NMPC for runaway-safe fed-batch reactors
Year: 2005
and [Milewska2006]Author: A. Milewska
School: Warsaw University of Technology
Title: Modelling of batch and semibatch chemical reactors -- safety aspects
Year: 2006
. 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 semibatch_taco.mod

# ----------------------------------------------------------------
# Semibatch esterification problem using AMPL and TACO
# (c) Christian Kirches, Sven Leyffer
# ----------------------------------------------------------------
include OptimalControl.mod;
 
var t;
 
# propionic anhydride (species A, sd0) [mol]
var nA := 0.0, >= 0.0, <= 10.0;
 
# 2-butanol (species B, sd1) [mol]
var nB := 6.893955747436589, >= 0.0, <= 10.0 suffix scale 10.0;
 
# propionic acid (species C, sd2) [mol]
var nC := 0.0, >= 0.0, <= 10.0;
 
# reactor temperature change (sd3) [K]
var Tr := 0.0, >= -23.0, <= 30.0 suffix scale 10.0;
 
# accumulated dosing of propionic anhydride [mol]
var nK1 := 0.0, >= 0.0, <= 7.0 suffix scale 5.0;
 
# theoretical adiabatic temperature rise (DAE state) [degrees Celsius]
var xA := 293.15, >= 270.15, <= 363.15 suffix type "dae" suffix scale 300.0;
 
# dosing feed [kg/s]
var dF := 0.0002472, >= 0.0, <= 0.0005 suffix type "u0" suffix scale 1.0e-2;
 
# dosing temperature [K]
var Tdos := 298.15, >= 298.15, <= 298.15 suffix type "u0" suffix scale 100.0;
 
# ambient temperature [K]
var Tamb := 298.15, >= 298.15, <= 298.15 suffix type "u0" suffix scale 100.0;
 
# jacket temperature change [K]
var dTj := 0.0, >= 0.0, <= 0.0 suffix type "u0";
 
param TJ    := 313.15;
param Tr_sp := 313.15;
param MA    :=   0.13014;
param MB    :=   0.07412;
param MC    :=   0.07408;
param MD    :=   0.13011;
param U1    := 195; 
param U2    := 155;
param V1u   :=   0.8;
param V2u   :=   1.6;
param nS0   :=   5.01/98.08;
param N1    :=   1.58;
param N2    :=   0.76;
param EA    :=   9.14E+4;
param A1    :=   9.35E+10;
param A2    :=   9.22E+10;
param A3    :=   9.78E+10;
param R     :=   8.314;
param V1    :=   1;
param V2    :=   2;
param CpI1  := 117.3;
param CpI2  := 198.4;
param Amin  :=   0.0113;
param Vmin  :=   0.124;
param d     :=   0.155;
param alpha :=   0.1;
param HA    := 59458;
 
var Tr_   = Tr + Tr_sp;
var dens0 = 1000*0.13014/(0.46199*0.2166^((1-Tr_/630.0)^0.2857));
var dens1 = 1000*0.07412/(0.3199 *0.2088^((1-Tr_/536.01)^0.2857));
var dens2 = 1000*0.07408/(0.2848 *0.1993^((1-Tr_/611.0)^0.2857));
var dens3 = 1000*0.13011/(0.5078 *0.2265^((1-Tr_/578.01)^0.2857));
var V     = (nA*MA/dens0 + nB*MB/dens1 + nC*MC/dens2 + nC*MD/dens3) * 1000;
var Ua    = U1 + (U2-U1)/(V2u-V1u)*(V-V1u);
var xA_   = nA/(nA+nB+2*nC+nS0);
var r     = ( A1*exp(-EA/(R*Tr_)) + A2*exp(-EA/(R*Tr_)) 
          * (nC/V)^N1 + A3*exp(-EA/(R*Tr_)) * (nS0/V)^N2 ) * (nA/V) * (nB/V);
var Cpi   = CpI1 + (CpI2 - CpI1)/(V2-V1) * (V-V1);
var Area  = (Amin + 4*(V-Vmin)/(1000*d));
var cp0   = 683.7305 - 3.7009 * Tr_ + 9.6E-3 * Tr_^2 - 7.458E-6 * Tr_^3;
var cp1   = 242.5725 - 2.3128 * Tr_ + 1.07E-2 * Tr_^2 - 1.159E-5 * Tr_^3;
var cp2   =  94.8369 - 0.2432 * Tr_ + 2.2E-3 * Tr_^2 - 2.51E-6 * Tr_^3;
var cp3   = 206.162  + 0.1994689 * Tr_ - 3.55512E-4 * Tr_^2 - 1.03625E-6 * Tr_^3;
var cpm   = nA*cp0 + nB*cp1 + nC*cp2 + nC*cp3;
var qdil  = -(-2232.74201/0.13963 * exp(-xA_/0.13963)) * (dF/MA-r*V);
 
minimize Tracking: 
	integral (nB^2,3600) suffix scale 1.0E+3;
 
subject to 
 
ODE_nA:  diff(nA,t) = dF/MA - r*V;
ODE_nB:  diff(nB,t) = -r*V;
ODE_nC:	 diff(nC,t) = r*V;
ODE_Tr:  diff(Tr,t) = ( r*HA*V - qdil - Ua*Area*(Tr_-(TJ+dTj)) 
                        - alpha*(Tr_-Tamb) - dF/MA*cp0*(Tr_-Tdos) ) / (Cpi+cpm);
ODE_nK1: diff(nK1,t) = dF/MA;
 
DAE_xA: 0 = xA - (Tr_ + nA * 1000 * HA/(cpm * V * 900));
 
IC_nA:  eval(nA,0) = 0;
IC_nB:  eval(nB,0) = 6.893955747436589;
IC_nC:  eval(nC,0) = 0;
IC_Tr:  eval(Tr,0) = 0;
IC_nK1: eval(nK1,0) = 0;
 
option solver ...;
 
solve;

Other Descriptions

Other descriptions of this problem are available in

• Mathematical notation at Semibatch esterification problem

References

[Kuehl2005]	P. K\"uhl; A. Milewska; M. Diehl; E. Molga; H.G. Bock (2005): NMPC for runaway-safe fed-batch reactors. %publisher%, Proc. Int. Workshop on Assessment and Future Directions of NMPC
[Milewska2006]	A. Milewska (2006): Modelling of batch and semibatch chemical reactors -- safety aspects. Warsaw University of Technology