Main Content

Design Optimization with Uncertain Variables (Code)

Open Script

This example shows how to optimize a design when there are uncertain variables. You optimize the dimensions of a Continuously Stirred Tank Reactor (CSTR) to minimize product concentration variation and production cost in case of varying, or uncertain, feed stock.

Continuously Stirred Tank Reactor (CSTR) Model

Continuously Stirred Tank Reactors (CSTRs) are common in the process industry. The Simulink® model, sdoCSTR, models a jacketed diabatic (i.e., non-adiabatic) tank reactor described in [1]. The CSTR is assumed to be perfectly mixed, with a single first-order exothermic and irreversible reaction, $A\rightarrow B$. $A$, the reactant, is converted to $B$, the product.

In this example, you use the following two-state CSTR model, which uses basic accounting and energy conservation principles:

$$\frac{d C_A}{dt} = \frac{F}{A * h}(C_{feed} - C_A)-r*C_A$$

$$\frac{d T}{dt} = \frac{F}{A * h}(T_{feed} - T) - \frac{H}{c_p \rho} r -\frac{U}{c_p*\rho*h}(T - T_{cool})$$

$$r = k_0*e^{\frac{-E}{R*T}}$$

  • $C_A$, and $C_{feed}$ - Concentrations of A in the CSTR and in the feed [kgmol/m^3]

  • $T$, $T_{feed}$, and $T_{cool}$ - CSTR, feed, and coolant temperatures [K]

  • $F$ and $\rho$ - Volumetric flow rate [m^3/h] and the density of the material in the CSTR [1/m^3]

  • $h$ and $A$ - Height [m] and heated cross-sectional area [m^2] of the CSTR.

  • $k_0$ - Pre-exponential non-thermal factor for reaction $A\rightarrow B$ [1/h]

  • $E$ and $H$ - Activation energy and heat of reaction for $A\rightarrow B$ [kcal/kgmol]

  • $R$ - Boltzmann's gas constant [kcal/(kgmol * K)]

  • $c_p$ and $U$ - Heat capacity [kcal/K] and heat transfer coefficients [kcal/(m^2 * K * h)]

Open the Simulink model.

open_system('sdoCSTR');

The model includes a cascaded PID controller in the Controller subsystem. The controller regulates the reactor temperature, $T$, and reactor residual concentration, $C_A$.

CSTR Design Problem

Assume that the CSTR is cylindrical, with the coolant applied to the base of the cylinder. Tune the CSTR cross-sectional area, $A$, and CSTR height, $h$, to meet the following design goals:

  • Minimize the variation in residual concentration, $C_A$. Variations in the residual concentration negatively affect the quality of the CSTR product. Minimizing the variations also improves CSTR profit.

  • Minimize the mean coolant temperature $T_{cool}$. Heating or cooling the jacket coolant temperature is expensive. Minimizing the mean coolant temperature improves CSTR profit.

The design must allow for variations in the quality of supply feed concentration, $C_{feed}$, and feed temperature, $T_{feed}$. The CSTR is fed with feed from different suppliers. The quality of the feed differs from supplier to supplier and also varies within each supply batch.

Specify Design Variables

Select the following model parameters as design variables for optimization:

  • Cylinder cross-sectional area $A$

  • Cylinder height $h$

p = sdo.getParameterFromModel('sdoCSTR',{'A','h'});

Limit the cross-sectional area to a range of [1 2] m^2.

p(1).Minimum = 1;
p(1).Maximum = 2;

Limit the height to a range of [1 3] m.

p(2).Minimum = 1;
p(2).Maximum = 3;

Specify Uncertain Variables

Select the feed concentration and feed temperature as uncertain variables. You evaluate the design using different values of feed temperature and concentration.

pUnc = sdo.getParameterFromModel('sdoCSTR',{'FeedCon0','FeedTemp0'});

Create a parameter space for the uncertain variables. Use normal distributions for both variables. Specify the mean as the current parameter value. Specify a variance of 5% of the mean for the feed concentration and 1% of the mean for the temperature.

uSpace = sdo.ParameterSpace(pUnc);
uSpace = setDistribution(uSpace,'FeedCon0',makedist('normal',pUnc(1).Value,0.05*pUnc(1).Value));
uSpace = setDistribution(uSpace,'FeedTemp0',makedist('normal',pUnc(2).Value,0.01*pUnc(2).Value));

The feed concentration is inversely correlated with the feed temperature. Add this information to the parameter space.

%uSpace.RankCorrelation = [1 -0.6; -0.6 1];

The rank correlation matrix has a row and column for each parameter with the (i,j) entry specifying the correlation between the i and j parameters.

Sample the parameter space. The scatter plot shows the correlation between concentration and temperature.

rng('default'); %For reproducibility
uSmpl = sdo.sample(uSpace,60);
sdo.scatterPlot(uSmpl)

Ideally you want to evaluate the design for every combination of points in the design and uncertain spaces, which implies 30*60 = 1800 simulations. Each simulation takes around 0.5 sec. You can use parallel computing to speed up the evaluation. For this example you instead only use the samples that have maximum & minimum concentration and temperature values, reducing the evaluation time to around 1 min.

[~,iminC] = min(uSmpl.FeedCon0);
[~,imaxC] = max(uSmpl.FeedCon0);
[~,iminT] = min(uSmpl.FeedTemp0);
[~,imaxT] = max(uSmpl.FeedTemp0);
uSmpl = uSmpl(unique([iminC,imaxC,iminT,imaxT]) ,:);

Specify Design Requirements

The design requirements require logging model signals. During optimization, the model is simulated using the current value of the design variables. Logged signals are used to evaluate the design requirements.

Log the following signals:

  • CSTR concentration, available at the second output port of the sdoCSTR/CSTR block

Conc = Simulink.SimulationData.SignalLoggingInfo;
Conc.BlockPath               = 'sdoCSTR/CSTR';
Conc.OutputPortIndex         = 2;
Conc.LoggingInfo.NameMode    = 1;
Conc.LoggingInfo.LoggingName = 'Concentration';

  • Coolant temperature, available at the first output of the sdoCSTR/Controller block

Coolant = Simulink.SimulationData.SignalLoggingInfo;
Coolant.BlockPath               = 'sdoCSTR/Controller';
Coolant.OutputPortIndex         = 1;
Coolant.LoggingInfo.NameMode    = 1;
Coolant.LoggingInfo.LoggingName = 'Coolant';

Create and configure a simulation test object to log the required signals.

simulator = sdo.SimulationTest('sdoCSTR');
simulator.LoggingInfo.Signals = [Conc,Coolant];

Create Objective/Constraint Function

Create a function to evaluate the CSTR design. This function is called at each optimization iteration.

Use an anonymous function with one argument that calls the sdoCSTR_design function.

evalDesign = @(p) sdoCSTR_design(p,simulator,pUnc,uSmpl);

The evalDesign function:

  • Has one input argument that specifies the CSTR dimensions

  • Returns the optimization objective value

The sdoCSTR_design function uses a for loop that iterates through the sample values specified for the feed concentration. Within the loop, the function:

  • Simulates the model using the current iterate, feed concentration, and feed temperature values

  • Calculates the residual concentration variation and coolant temperature costs

To view the objective function, type edit sdoCSTR_design.

Evaluate Initial Design

Call the evalDesign function with the initial CSTR dimensions.

dInit = evalDesign(p)

dInit = 

  struct with fields:

              F: 11.3356
       costConc: 6.4390
    costCoolant: 4.8965

Plot the model response for the initial design. Simulate the model using the sample feed concentration values. The plot shows the variation in the residual concentration and coolant temperature.

sdoCSTR_plotModelResponse(p,simulator,pUnc,uSmpl);

The sdoCSTR_plotModelResponse function plots the model response. To view this function, type edit sdoCSTR_plotModelResponse.

Optimize Design

Pass the objective function and initial CSTR dimensions to sdo.optimize.

pOpt = sdo.optimize(evalDesign,p)

 Optimization started 2023-Aug-19, 15:59:35

                               max                     First-order 
 Iter F-count        f(x)   constraint    Step-size    optimality
    0      4      5.17535            0
    1      8      3.81522            0         2.01         7.83
    2     12      2.63963            0         0.57         3.03
    3     16      2.53426            0        0.159        0.292
    4     20      2.50995            0        0.017        0.432
    5     24      2.50266            0        0.125         1.44
    6     28      2.50097            0       0.0966        0.234
    7     32      2.47473            0       0.0135       0.0719
    8     37      2.47473            0     0.000955       0.0694
Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.
 
pOpt(1,1) =
 
       Name: 'A'
      Value: 2
    Minimum: 1
    Maximum: 2
       Free: 1
      Scale: 0.5000
       Info: [1x1 struct]

 
pOpt(2,1) =
 
       Name: 'h'
      Value: 2.2850
    Minimum: 1
    Maximum: 3
       Free: 1
      Scale: 2
       Info: [1x1 struct]

 
2x1 param.Continuous

Evaluate Optimized Design

Call the evalDesign function with the optimized CSTR dimensions.

dFinal = evalDesign(pOpt)

dFinal = 

  struct with fields:

              F: 2.4747
       costConc: 1.4505
    costCoolant: 1.0242

Plot the model response for the optimized design. Simulate the model using the sample feed concentration values. The optimized design reduces the residual concentration variation and average coolant temperature for different feed stocks.

sdoCSTR_plotModelResponse(pOpt,simulator,pUnc,uSmpl);

Related Examples

To learn how to use sensitivity analysis to explore the CSTR design space and select an initial design for optimization, see Design Exploration Using Parameter Sampling (Code).

References

[1] Bequette, B.W. Process Dynamics: Modeling, Analysis and Simulation. 1st ed. Upper Saddle River, NJ: Prentice Hall, 1998.

% Close the model
bdclose('sdoCSTR')

See Also

| |

Related Topics