This example shows how to use multiple MPC controllers to control a nonlinear continuous stirred tank reactor (CSTR) as it transitions from low conversion rate to high conversion rate.

Multiple MPC Controllers are designed at different operating conditions and then implemented with the Multiple MPC Controller block in Simulink. At run time, a scheduling signal is used to switch controller from one to another.

A Continuously Stirred Tank Reactor (CSTR) is a common chemical system in the process industry. A schematic of the CSTR system is:

This is a jacketed non-adiabatic tank reactor described extensively in Seborg's book, "Process Dynamics and Control", published by Wiley, 2004. The vessel is assumed to be perfectly mixed, and a single first-order exothermic and irreversible reaction, A --> B, takes place. The inlet stream of reagent A is fed to the tank at a constant volumetric rate. The product stream exits continuously at the same volumetric rate and liquid density is constant. Thus the volume of reacting liquid is constant.

The inputs of the CSTR model are:

and the outputs (y(t)), which are also the states of the model (x(t)), are:

The control objective is to maintain the concentration of reagent A, at its desired setpoint, which changes over time when reactor transitions from low conversion rate to high conversion rate. The coolant temperature is the manipulated variable used by the MPC controller to track the reference. The inlet feed stream concentration and temperature are assumed to be constant. The Simulink model `mpc_cstr_plant`

implements the nonlinear CSTR plant.

It is well known that the CSTR dynamics are strongly nonlinear with respect to reactor temperature variations and can be open-loop unstable during the transition from one operating condition to another. A single MPC controller designed at a particular operating condition cannot give satisfactory control performance over a wide operating range.

To control the nonlinear CSTR plant with linear MPC control technique, you have a few options:

If a linear plant model cannot be obtained at run time, first you need to obtain several linear plant models offline at different operating conditions that cover the typical operating range. Next you can choose one of the two approaches to implement MPC control strategy:

(1) Design several MPC controllers offline, one for each plant model. At run time, use Multiple MPC Controller block that switches MPC controllers from one to another based on a desired scheduling strategy, as discussed in this example. Use this approach when the plant models have different orders or time delays.

(2) Design one MPC controller offline at a nominal operating point. At run time, use Adaptive MPC Controller block (updating predictive model at each control interval) together with Linear Parameter Varying (LPV) System block (supplying linear plant model with a scheduling strategy). See Adaptive MPC Control of Nonlinear Chemical Reactor Using Linear Parameter Varying System for more details. Use this approach when all the plant models have the same order and time delay.

If a linear plant model can be obtained at run time, you should use Adaptive MPC Controller block to achieve nonlinear control. There are two typical ways to obtain a linear plant model online:

(1) Use successive linearization. See Adaptive MPC Control of Nonlinear Chemical Reactor Using Successive Linearization for more details. Use this approach when a nonlinear plant model is available and can be linearized at run time.

(2) Use online estimation to identify a linear model when loop is closed. See Adaptive MPC Control of Nonlinear Chemical Reactor Using Online Model Estimation for more details. Use this approach when linear plant model cannot be obtained from either an LPV system or successive linearization.

To run this example, Simulink® and Simulink Control Design® are required.

if ~mpcchecktoolboxinstalled('simulink') disp('Simulink(R) is required to run this example.') return end if ~mpcchecktoolboxinstalled('slcontrol') disp('Simulink Control Design(R) is required to run this example.') return end

First, a linear plant model is obtained at the initial operating condition, CAi is 10 kgmol/m^3, Ti and Tc are 298.15 K. Functions from Simulink Control Design such as `operspec`

, `findop`

, and `linearize`

, are used to generate the linear state-space system from the Simulink model.

Create operating point specification.

```
plant_mdl = 'mpc_cstr_plant';
op = operspec(plant_mdl);
```

Feed concentration is known at the initial condition.

op.Inputs(1).u = 10; op.Inputs(1).Known = true;

Feed temperature is known at the initial condition.

op.Inputs(2).u = 298.15; op.Inputs(2).Known = true;

Coolant temperature is known at the initial condition.

op.Inputs(3).u = 298.15; op.Inputs(3).Known = true;

Compute initial condition.

```
[op_point,op_report] = findop(plant_mdl,op);
% Obtain nominal values of x, y and u.
x0 = [op_report.States(1).x; op_report.States(2).x];
y0 = [op_report.Outputs(1).y; op_report.Outputs(2).y];
u0 = [op_report.Inputs(1).u; op_report.Inputs(2).u; op_report.Inputs(3).u];
```

Operating point search report: --------------------------------- Operating point search report for the Model mpc_cstr_plant. (Time-Varying Components Evaluated at time t=0) Operating point specifications were successfully met. States: ---------- (1.) mpc_cstr_plant/CSTR/Integrator x: 311 dx: 8.12e-11 (0) (2.) mpc_cstr_plant/CSTR/Integrator1 x: 8.57 dx: -6.87e-12 (0) Inputs: ---------- (1.) mpc_cstr_plant/CAi u: 10 (2.) mpc_cstr_plant/Ti u: 298 (3.) mpc_cstr_plant/Tc u: 298 Outputs: ---------- (1.) mpc_cstr_plant/T y: 311 [-Inf Inf] (2.) mpc_cstr_plant/CA y: 8.57 [-Inf Inf]

Obtain linear model at the initial condition.

plant = linearize(plant_mdl,op_point);

Verify that the linear model is open-loop stable at this condition.

eig(plant)

ans = -0.5223 -0.8952

You design an MPC at the initial operating condition.

Ts = 0.5;

Specify signal types used in MPC. Assume both reactor temperature and concentration are measurable.

plant.InputGroup.UnmeasuredDisturbances = [1 2]; plant.InputGroup.ManipulatedVariables = 3; plant.OutputGroup.Measured = [1 2]; plant.InputName = {'CAi','Ti','Tc'}; plant.OutputName = {'T','CA'};

Create MPC controller with default prediction and control horizons

mpcobj = mpc(plant,Ts);

-->The "PredictionHorizon" property of "mpc" object is empty. Trying PredictionHorizon = 10. -->The "ControlHorizon" property of the "mpc" object is empty. Assuming 2. -->The "Weights.ManipulatedVariables" property of "mpc" object is empty. Assuming default 0.00000. -->The "Weights.ManipulatedVariablesRate" property of "mpc" object is empty. Assuming default 0.10000. -->The "Weights.OutputVariables" property of "mpc" object is empty. Assuming default 1.00000. for output(s) y1 and zero weight for output(s) y2

Set nominal values in the controller. Note that nominal values for unmeasured disturbance must be zero.

mpcobj.Model.Nominal = struct('X',x0,'U',[0;0;u0(3)],'Y',y0,'DX',[0 0]);

Set scale factors because plant input and output signals have different orders of magnitude.

Uscale = [10;30;50]; Yscale = [50;10]; mpcobj.DV(1).ScaleFactor = Uscale(1); mpcobj.DV(2).ScaleFactor = Uscale(2); mpcobj.MV.ScaleFactor = Uscale(3); mpcobj.OV(1).ScaleFactor = Yscale(1); mpcobj.OV(2).ScaleFactor = Yscale(2);

The goal will be to track a specified transition in the reactor concentration. The reactor temperature will be measured and used in state estimation but the controller will not attempt to regulate it directly. It will vary as needed to regulate the concentration. Thus, set its MPC weight to zero.

mpcobj.Weights.OV = [0 1];

Plant inputs 1 and 2 are unmeasured disturbances. By default, the controller assumes integrated white noise with unit magnitude at these inputs when configuring the state estimator. Try increasing the state estimator signal-to-noise by a factor of 10 to improve disturbance rejection performance.

```
Dist = ss(getindist(mpcobj));
Dist.B = eye(2)*10;
setindist(mpcobj,'model',Dist);
```

-->Converting model to discrete time. -->The "Model.Disturbance" property of "mpc" object is empty: Assuming unmeasured input disturbance #1 is integrated white noise. Assuming unmeasured input disturbance #2 is integrated white noise. Assuming no disturbance added to measured output channel #2. Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

All other MPC parameters are at their default values.

"mpc_cstr_single" contains a Simulink® model with CSTR and MPC Controller blocks in a feedback configuration.

```
mpc_mdl = 'mpc_cstr_single';
open_system(mpc_mdl)
```

Note that the MPC Controller block is configured to look ahead at (preview) setpoint changes in the future; that is, anticipating the setpoint transition. This generally improves setpoint tracking.

Define a constant setpoint for the output.

CSTR_Setpoints.time = [0; 60]; CSTR_Setpoints.signals.values = [y0 y0]';

Test the response to a 5% increase in feed concentration.

set_param([mpc_mdl '/Feed Concentration'],'Value','10.5');

Set plot scales and simulate the response.

open_system([mpc_mdl '/Measurements']) open_system([mpc_mdl '/Coolant Temperature']) set_param([mpc_mdl '/Measurements'],'Ymin','305~8','Ymax','320~9') set_param([mpc_mdl '/Coolant Temperature'],'Ymin','295','Ymax','305') sim(mpc_mdl,10)

-->Converting model to discrete time. Assuming no disturbance added to measured output channel #2. Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

The closed-loop response is satisfactory.

First, define the desired setpoint transition. After a 10-minute warm-up period, ramp the concentration setpoint downward at a rate of 0.25 per minute until it reaches 2.0 kmol/m^3.

CSTR_Setpoints.time = [0 10 11:39]'; CSTR_Setpoints.signals.values = [y0(1)*ones(31,1),[y0(2);y0(2);(y0(2):-0.25:2)';2;2]];

Remove the 5% increase in feed concentration used previously.

set_param([mpc_mdl '/Feed Concentration'],'Value','10')

Set plot scales and simulate the response.

set_param([mpc_mdl '/Measurements'],'Ymin','300~0','Ymax','400~10') set_param([mpc_mdl '/Coolant Temperature'],'Ymin','240','Ymax','360')

Simulate model.

sim(mpc_mdl,60)

The closed-loop response is unacceptable. Performance along the full transition can be improved if other MPC controllers are designed at different operating conditions along the transition path. In the next two section, two additional MPC controllers are design at intermediate and final transition stages respectively.

Create operating point specification.

op = operspec(plant_mdl);

Feed concentration is known.

op.Inputs(1).u = 10; op.Inputs(1).Known = true;

Feed temperature is known.

op.Inputs(2).u = 298.15; op.Inputs(2).Known = true;

Reactor concentration is known.

op.Outputs(2).y = 5.5; op.Outputs(2).Known = true;

Find steady state operating condition.

```
[op_point,op_report] = findop(plant_mdl,op);
% Obtain nominal values of x, y and u.
x0 = [op_report.States(1).x; op_report.States(2).x];
y0 = [op_report.Outputs(1).y; op_report.Outputs(2).y];
u0 = [op_report.Inputs(1).u; op_report.Inputs(2).u; op_report.Inputs(3).u];
```

Operating point search report: --------------------------------- Operating point search report for the Model mpc_cstr_plant. (Time-Varying Components Evaluated at time t=0) Operating point specifications were successfully met. States: ---------- (1.) mpc_cstr_plant/CSTR/Integrator x: 339 dx: 3.42e-08 (0) (2.) mpc_cstr_plant/CSTR/Integrator1 x: 5.5 dx: -2.87e-09 (0) Inputs: ---------- (1.) mpc_cstr_plant/CAi u: 10 (2.) mpc_cstr_plant/Ti u: 298 (3.) mpc_cstr_plant/Tc u: 298 [-Inf Inf] Outputs: ---------- (1.) mpc_cstr_plant/T y: 339 [-Inf Inf] (2.) mpc_cstr_plant/CA y: 5.5 (5.5)

Obtain linear model at the initial condition.

plant_intermediate = linearize(plant_mdl,op_point);

Verify that the linear model is open-loop unstable at this condition.

eig(plant_intermediate)

ans = 0.4941 -0.8357

Specify signal types used in MPC. Assume both reactor temperature and concentration are measurable.

plant_intermediate.InputGroup.UnmeasuredDisturbances = [1 2]; plant_intermediate.InputGroup.ManipulatedVariables = 3; plant_intermediate.OutputGroup.Measured = [1 2]; plant_intermediate.InputName = {'CAi','Ti','Tc'}; plant_intermediate.OutputName = {'T','CA'};

Create MPC controller with default prediction and control horizons.

mpcobj_intermediate = mpc(plant_intermediate,Ts);

-->The "PredictionHorizon" property of "mpc" object is empty. Trying PredictionHorizon = 10. -->The "ControlHorizon" property of the "mpc" object is empty. Assuming 2. -->The "Weights.ManipulatedVariables" property of "mpc" object is empty. Assuming default 0.00000. -->The "Weights.ManipulatedVariablesRate" property of "mpc" object is empty. Assuming default 0.10000. -->The "Weights.OutputVariables" property of "mpc" object is empty. Assuming default 1.00000. for output(s) y1 and zero weight for output(s) y2

Set nominal values, scale factors and weights in the controller.

mpcobj_intermediate.Model.Nominal = struct('X',x0,'U',[0;0;u0(3)],'Y',y0,'DX',[0 0]); Uscale = [10;30;50]; Yscale = [50;10]; mpcobj_intermediate.DV(1).ScaleFactor = Uscale(1); mpcobj_intermediate.DV(2).ScaleFactor = Uscale(2); mpcobj_intermediate.MV.ScaleFactor = Uscale(3); mpcobj_intermediate.OV(1).ScaleFactor = Yscale(1); mpcobj_intermediate.OV(2).ScaleFactor = Yscale(2); mpcobj_intermediate.Weights.OV = [0 1]; Dist = ss(getindist(mpcobj_intermediate)); Dist.B = eye(2)*10; setindist(mpcobj_intermediate,'model',Dist);

-->Converting model to discrete time. -->The "Model.Disturbance" property of "mpc" object is empty: Assuming unmeasured input disturbance #1 is integrated white noise. Assuming unmeasured input disturbance #2 is integrated white noise. Assuming no disturbance added to measured output channel #2. Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Create operating point specification.

op = operspec(plant_mdl);

Feed concentration is known.

op.Inputs(1).u = 10; op.Inputs(1).Known = true;

Feed temperature is known.

op.Inputs(2).u = 298.15; op.Inputs(2).Known = true;

Reactor concentration is known.

op.Outputs(2).y = 2; op.Outputs(2).Known = true;

Find steady-state operating condition.

```
[op_point,op_report] = findop(plant_mdl,op);
% Obtain nominal values of x, y and u.
x0 = [op_report.States(1).x; op_report.States(2).x];
y0 = [op_report.Outputs(1).y; op_report.Outputs(2).y];
u0 = [op_report.Inputs(1).u; op_report.Inputs(2).u; op_report.Inputs(3).u];
```

Operating point search report: --------------------------------- Operating point search report for the Model mpc_cstr_plant. (Time-Varying Components Evaluated at time t=0) Operating point specifications were successfully met. States: ---------- (1.) mpc_cstr_plant/CSTR/Integrator x: 373 dx: 5.57e-11 (0) (2.) mpc_cstr_plant/CSTR/Integrator1 x: 2 dx: -4.6e-12 (0) Inputs: ---------- (1.) mpc_cstr_plant/CAi u: 10 (2.) mpc_cstr_plant/Ti u: 298 (3.) mpc_cstr_plant/Tc u: 305 [-Inf Inf] Outputs: ---------- (1.) mpc_cstr_plant/T y: 373 [-Inf Inf] (2.) mpc_cstr_plant/CA y: 2 (2)

Obtain linear model at the initial condition.

plant_final = linearize(plant_mdl,op_point);

Verify that the linear model is again open-loop stable at this condition.

eig(plant_final)

ans = -1.1077 + 1.0901i -1.1077 - 1.0901i

Specify signal types used in MPC. Assume both reactor temperature and concentration are measurable.

plant_final.InputGroup.UnmeasuredDisturbances = [1 2]; plant_final.InputGroup.ManipulatedVariables = 3; plant_final.OutputGroup.Measured = [1 2]; plant_final.InputName = {'CAi','Ti','Tc'}; plant_final.OutputName = {'T','CA'};

Create MPC controller with default prediction and control horizons.

mpcobj_final = mpc(plant_final,Ts);

-->The "PredictionHorizon" property of "mpc" object is empty. Trying PredictionHorizon = 10. -->The "ControlHorizon" property of the "mpc" object is empty. Assuming 2. -->The "Weights.ManipulatedVariables" property of "mpc" object is empty. Assuming default 0.00000. -->The "Weights.ManipulatedVariablesRate" property of "mpc" object is empty. Assuming default 0.10000. -->The "Weights.OutputVariables" property of "mpc" object is empty. Assuming default 1.00000. for output(s) y1 and zero weight for output(s) y2

Set nominal values, scale factors and weights in the controller.

mpcobj_final.Model.Nominal = struct('X',x0,'U',[0;0;u0(3)],'Y',y0,'DX',[0 0]); Uscale = [10;30;50]; Yscale = [50;10]; mpcobj_final.DV(1).ScaleFactor = Uscale(1); mpcobj_final.DV(2).ScaleFactor = Uscale(2); mpcobj_final.MV.ScaleFactor = Uscale(3); mpcobj_final.OV(1).ScaleFactor = Yscale(1); mpcobj_final.OV(2).ScaleFactor = Yscale(2); mpcobj_final.Weights.OV = [0 1]; Dist = ss(getindist(mpcobj_final)); Dist.B = eye(2)*10; setindist(mpcobj_final,'model',Dist);

-->Converting model to discrete time. -->The "Model.Disturbance" property of "mpc" object is empty: Assuming unmeasured input disturbance #1 is integrated white noise. Assuming unmeasured input disturbance #2 is integrated white noise. Assuming no disturbance added to measured output channel #2. Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

The following model uses the Multiple MPC Controllers block to implement three MPC controllers across the operating range.

```
mmpc_mdl = 'mpc_cstr_multiple';
open_system(mmpc_mdl)
```

Note that it has been configured to use the three controllers in a sequence: `mpcobj`

, `mpcobj_intermediate`

, and `mpcobj_final`

.

```
open_system([mmpc_mdl '/Multiple MPC Controllers'])
```

Note also that the two switches specify when to switch from one controller to another. The rules are:

If CSTR concentration >= 8, use

`mpcobj`

.If 3 <= CSTR concentration < 8, use

`mpcobj_intermediate`

.If CSTR concentration < 3, use

`mpcobj_final`

.

Simulate with the Multiple MPC Controllers block

open_system([mmpc_mdl '/Measurements']) open_system([mmpc_mdl '/MV']) sim(mmpc_mdl)

-->Converting model to discrete time. Assuming no disturbance added to measured output channel #2. Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel. -->Converting model to discrete time. Assuming no disturbance added to measured output channel #2. Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

The transition is now well controlled. The major improvement is in the transition through the open-loop unstable region. The plot of the switching signal shows when controller transitions occur. The MV character changes at these times because of the change in dynamic characteristics introduced by the new prediction model.

bdclose(plant_mdl) bdclose(mpc_mdl) bdclose(mmpc_mdl)

MPC Controller | Multiple MPC Controllers