Optimal Control for a Continuously Stirred Tank Reactor

Introduction

Find an optimal control signal for a continuously stirred tank reactor (CSTR). The temperature in the tank influences a reaction rate, while the temperature in the cooling jacket acts as a control signal. The goal is to go from one concentration to another as quickly as possible, while obeying safety limitations in temperature.

Explore the available scenarios by evaluating each cell in order. Evaluation is done by selecting a cell and pressing Shift+Enter.

In order to get the full experience of this example, you need the following:

Wolfram SystemModeler; a free trial download is available at www.wolfram.com/system-modeler

The model that came in the same bundle as this notebook

For more information about Wolfram SystemModeler, visit www.wolfram.com/system-modeler.

Initialize

Define the model name:

mmodel="EducationExamples.ChemicalEngineering.ChemicalTankReactor";

Optimal Control Design

Discretize the equations from the model, using forward difference for differentiation. Optimize the cooling jacket temperature control signal to achieve the target concentration as fast as possible while obeying 
limitations in tank temperature.

Set parameters used in the solver process:

steps =150; (* The number of steps *)
stepsize=150/steps;(* The step length *)
Tc0A=250; (*Cooling temperature at point A*)
TrefA=250.052; (*Temperature at point A, start temperature*)
crefA=956.209; (*Concentration at point A, start concentration*)
TcrefA=Tc0A;
Tc0B=280; (*Cooling temperature at point B*)
TrefB=280.099; (*Reference temperature at point B*)
crefB=338.691; (*Reference concentration at point B*)
TcrefB=Tc0B;

Retrieve system equations for the model:

eqns=SystemModel[mmodel,"SystemEquations",t];

Replace the parameters with their numerical values:

neqns=(eqns/. SystemModel[mmodel,"ParameterValues"])/.QuantityVariable[var:_,_]:>Symbol[var]

Replace the derivatives with forward differences:

{ceqn,teqn}=
neqns/.{c'[t]->(c[t+1]-c[t])/stepsize, T'[t]->(T[t+1]-T[t])/stepsize}

Create constraints by discretizing into steps and applying And for each step:

cTconstraints = And[And@@Table[ceqn,{t,2,steps-1}],And@@Table[teqn,{t,1,steps-1}]];

Define the cost function to use:

costFn=Sum[(crefB-c[i])^2+(TrefB-T[i])^2+(TcrefB-Tc[i])^2,{i,steps}];

Create constraints for the temperatures, requiring them to stay inside given safety limits:

TLim=And@@Table[T[i]<=350,{i,2,steps}];
TcLim=And@@Table[230≤Tc[i]<=370,{i,1,steps}];

Define and apply initial conditions:

c2 =Solve[ceqn,c[1+t]]/.{c[t]->crefA , T[t]->TrefA}  /.t->1

initialConditions={T[1]->TrefA,c[1]->crefA,c2}//Flatten

costFn=costFn/.initialConditions;

constr= Join[cTconstraints,TLim,TcLim]/.initialConditions;

Find the optimal control signal with FindMinimum. The call to FindMinimum takes about a minute, depending on your hardware.

allVars= Join[Table[{Tc[i],250},{i,1,steps}], Table[{c[i],956},{i,3,steps}], Table[{T[i],253},{i,2,steps}]];

First[AbsoluteTiming[forwardm=FindMinimum[{costFn,constr },allVars][[2]]]]

Collect the results in three signal variables:

forwardfc=Table[c[i],{i,3,steps}]/.forwardm;

forwardfT=Table[T[i],{i,2,steps}]/.forwardm;

forwardfTc=Table[Tc[i],{i,1,steps}]/.forwardm;

Show the results in three plots:

forwardplota=Show[ListLinePlot[forwardfc,PlotLabel->"Concentration",PlotStyle->Thick],Plot[{crefB,crefA},{t,0,steps},PlotStyle->{ColorData[86][3],ColorData[86][2]}]]

forwardplotb=Show[ListLinePlot[forwardfT,PlotLabel->"Temperature, T",PlotStyle->Thick],Plot[{TrefB,TrefA},{t,0,steps},PlotStyle->{ColorData[86][3],ColorData[86][2]}],PlotRange->{240,All},AxesOrigin->{0,240}]

forwardplotc=Show[ListLinePlot[forwardfTc,PlotLabel->"Cooling Temperature, Tc",PlotStyle->Thick],Plot[{TcrefB,TcrefA},{t,0,steps},PlotStyle->{ColorData[86][3],ColorData[86][2]}],PlotRange->{230,All},AxesOrigin->{0,240}]

Interactive Control Design

Create an interactive interface to design the cooling temperature control signal by dragging points on a plot:

optimalSolutionPlot=ListLinePlot[forwardfc,PlotStyle->Directive[Opacity[0.5],Thick,Dashed]];
simulateTank[]:=
Module[{},pts=Sort[pts];
pts=pts+Prepend[Differences[pts[[All,1]]]/. {0.-> 0.1,_?(#>0&)-> 0},0];
intfn=Interpolation[pts,InterpolationOrder->1];
sol=NDSolve[{neqns/.Tc->intfn,T[1]==TrefA,c[1]==crefA},{c,T},{t,1,150}];
maxT=Max[Table[First[T[t]/.sol],{t,1,150}]];];
ctrlPoints=8;
pts=Table[{i,forwardfTc[[i]]},{i,Round[Subdivide[1,150,ctrlPoints]]}];
simulateTank[];

Dynamic[
Grid[{
{
(* The control plot *)
LocatorPane[
Dynamic[pts,
{
(pts=ReplacePart[#,{{1,1}->1,{-1,1}-> 150}])&,
(
simulateTank[]
)&
}],
ListLinePlot[intfn,PlotRange->{{1,150},{230,370}},ImageSize->Medium]]},
{
(* The plot with the tank temperature *)
Plot[Evaluate[{T[t]/.sol,350}],{t,1,150},PlotRange->{230,365},AxesOrigin->{0,240},PlotLabel->"Temperature",ImageSize->Medium,PlotLegends->Placed[{"Temperature Trajectory","Max Temperature"},Below]],
(* The resulting concentration plot *)
Show[
Plot[
Evaluate[{c[t],crefA,crefB}/.sol],{t,1,150},
PlotRange->{230,1000},AxesOrigin->{0,240},PlotLabel->"Concentration",ImageSize->Medium,
PlotLegends-> Placed[LineLegend[{
ColorData[86][1],Directive[ColorData[86][1],Dashed], ColorData[86][2],ColorData[86][3]},{"Concentration Trajectory","Optimal Concentration Trajectory","Initial Concentraion","Reference Concentration"},LegendLayout-> "Column"],Below]],
optimalSolutionPlot]
}
}]
]

Terms and Conditions of Use

This domain example is an informational resource made freely available by Wolfram Research.

Use of This Example

• You may not use this example for any purpose that is unlawful or dangerous.
• You assume total responsibility and risk for your use of this example.
• You may not present this example with any alteration to its trade dress, trademarks, or design.

License

All content in this bundle is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License. By accessing the content or using it in any way, you accept and agree to be bound by the terms of this license. If you do not agree to these Terms of Use, you may not use this content. Wolfram Research reserves the right to change, modify, add to, or remove portions of these Terms of Use at any time without notice. Please refer back to www.wolfram.com for the latest Terms of Use.

A summary of the licensing terms can be found at:
http://creativecommons.org/licenses/by-nc-sa/3.0

The full legal code can be found at: 
http://creativecommons.org/licenses/by-nc-sa/3.0/legalcode