MATLAB Examples

Economic MPC Control of Ethylene Oxide Production

This example shows how to maximize the production of an ethylene oxide plant for profit using an economic MPC controller with custom cost and constraint functions.

Contents

This example requires Simulink® software to simulate nonlinear MPC control of the ethylene oxidation plant in Simulink.

if ~mpcchecktoolboxinstalled('simulink')
    disp('Simulink(TR) is required to run this example.')
    return
end

The example also uses the fsolve and fmincon commands from the Optimization Toolbox.

if ~mpcchecktoolboxinstalled('optim')
    disp('Optimization Toolbox must be installed to run this example.')
    return
end

Add example file folder to MATLAB® path.

addpath(fullfile(matlabroot,'examples','mpc_featured','main'));

Nonlinear Ethylene Oxidation Plant

Conversion of ethylene (C2H4) to ethylene oxide (C2H4O) occurs in a cooled, gas-phase catalytic reactor. Three reactions occur simultaneously in the well-mixed gas phase within the reactor:

C2H4 + 0.5*O2 -> C2H4O
C2H4 + 3*O2 -> 2*CO2 + 2*H2O
C2H4O + 2.5*O2 -> 2*CO2 + 2*H2O

The first reaction is wanted and the other two reduce C2H4O production. A mixture of air and ethylene is continuously fed into the tank. The first-principle nonlinear dynamic model of the reactor is implemented as a set of ordinary differential equations (ODEs) in the oxidationPlantCT function. For more information, see oxidationPlantCT.m.

The plant has six states:

  • Gas density in the reactor ($x_1$)
  • C2H4 concentration in the reactor ($x_2$)
  • C2H4O concentration in the reactor ($x_3$)
  • Temperature in the reactor ($x_4$)
  • Total volumetric flow rate exiting the reactor ($x_5$)
  • Integrator state of a PI reactor temperature controller ($x_6$)

The plant has three inputs:

  • Total volumetric feed flow rate ($u_1$)
  • C2H4 concentration in the feed ($u_2$)
  • Reactor temperature setpoint ($u_3$)

The plant has two outputs:

  • C2H4O concentration in the reactor ($y_1$, equivalent to $x_3$)
  • Total volumetric flow rate exiting the reactor ($y_2$, equivalent to $x_5$)

All variables in the model are scaled to be dimensionless and of unity order. The basic plant equations and parameters are obtained from [1].

The reactions are exothermic. Therefore, continuous feedback control of the reactor temperature has been added, as would normally be done for safe and stable operation. To hold the reactor temperature at a specified setpoint, the PI control law manipulates the coolant temperature.

The plant is asymptotically open-loop stable.

Control Objectives

The primary control objective is to maximize the time-averaged ethylene oxide (C2H4O) production rate, which in turn maximizes profit. The instantaneous C2H4O production rate is ${y_1}{y_2}$.

Also, during the production operation, the ethylene feed rate, ${u_1}{u_2}$, must always equal the available ethylene coming from an upstream process. This equality constraint is referred to as C2H4 availability in this example.

Therefore, the MPC controller must maximize the C2H4O production rate from the three chemical reactions under the C2H4 availability constraint by manipulating $u_1$ and $u_2$.

Nominal Condition

At the nominal condition, the C2H4 availability (${u_1}{u_2}$) is 0.175, and the reactor is operating at a steady state with plant inputs uNom.

uNom = [0.5; 0.35; 1.15];

To obtain the nominal operating point, compute the nominal plant states based on the given plant inputs usi fsolve from the Optimization Toolbox.

options = optimoptions('fsolve','Display','none');
x0 = [0.9,0.4,0.1,1.1,0.3,0];  % Initial guess
xNom = fsolve(@(x) oxidationPlantCT(x,uNom),x0,options);
[dxNom,yNom] = oxidationPlantCT(xNom,uNom);

Verify that all time derivatives are close to zero at the nominal operating point.

norm(dxNom)
ans =

   8.8674e-09

Linear Plant Model

To obtain the nominal LTI plant model used by the MPC controller, linearize the ODEs by numerically perturbing the nominal states and inputs.

plant = linearizeOxidationPlant(@oxidationPlantCT,xNom,uNom);

Then, add C2H4 availability as an additional input ($u_4$). This input has no effect on the plant states or outputs. Instead, it is a measured disturbance used to define the equality constraint.

plant = ss(plant.A, [plant.B zeros(6,1)], plant.C, [plant.D zeros(2,1)]);

The nominal value of $u_4$ is ${u_1}{u_2}$.

uNom = [uNom; uNom(1)*uNom(2)];

Therefore, the MPC controller uses four plant inputs. The first two are manipulated variables, and the last two are measured disturbances.

plant.InputName = {'Qin','CEin','Tref','Availability'};
plant = setmpcsignals(plant,'MV',[1 2],'MD',[3 4]);

Specify names for the six states and two measured outputs.

plant.StateName = {'Den','C2H4','C2H4O','Tc','Qout','PIx'};
plant.OutputName = {'C2H4O','Qout'};

Traditional MPC Design

To create an economic MPC controller, first design a traditional MPC controller. The relatively large sample time of 5 seconds used here is appropriate when the plant is stable and the primary objective is economic optimization.

Ts = 5;     % Sample time
P = 20;     % Prediction horizon
M = 5;      % Control horizon
mpcobj = mpc(plant,Ts,P,M);
-->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.

Specify the nominal operating conditions.

mpcobj.Model.Nominal.X = xNom;
mpcobj.Model.Nominal.U = uNom;
mpcobj.Model.Nominal.Y = yNom;

Plant inputs $u_1$ and $u_2$ must stay within saturation limits:

$$\begin{array}{l}
0.05 \le {u_1} \le 0.7\\
0.1 \le {u_2} \le 3
\end{array}$$

mpcobj.MV(1).Min = 0.05;
mpcobj.MV(1).Max = 0.7;
mpcobj.MV(2).Min = 0.1;
mpcobj.MV(2).Max = 3;

The rates of change of $u_1$ and $u_2$ are also limited:

$$\begin{array}{l}
- 0.05 \le {u_1}\left( k \right) - {u_1}\left( {k - 1} \right) \le 0.0.5\\
- 0.2 \le {u_2}\left( k \right) - {u_2}\left( {k - 1} \right) \le 0.2
\end{array}$$

mpcobj.MV(1).RateMin = -0.05;
mpcobj.MV(1).RateMax =  0.05;
mpcobj.MV(2).RateMin = -0.2;
mpcobj.MV(2).RateMax =  0.2;

Economic MPC Design with Custom Cost and Constraints

Instead of using the standard quadratic objective function, use a custom cost function that can represent any arbitrary nonlinear cost. The custom cost function must be named mpcCustomCostFcn.m. The template mpcCustomCostFcn.txt provides a basis for customization. For details, see the comments in mpcCustomCostFcn.m.

In this example, the custom cost function is the C2H4O production rate time-averaged over the prediction horizon:

f = -sum(y(:,1).*y(:,2))/P;

where P is the MPC prediction horizon. The negative sign in f is used to maximize production, since the controller minimizes f during optimization. The custom cost function replaces the standard quadratic cost used in traditional MPC.

To force ethylene consumption ${u_1}{u_2}$ to equal C2H4 availability, use a custom constraint function. C2H4 availability is the second measured disturbance. The controller passes the measured disturbances to the custom constraint function in argument v. For this example, the second measured disturbance, v(:,2), corresponds to plant input $u_4$. Apply this constraint at each step of the prediction horizon, producing P nonlinear equality constraints.

Ceq = u(:,1).*u(:,2) - v(1:end-1,2);

The custom constraint function must be named mpcCustomConstraintFcn.m. The template mpcCustomConstraint.txt provides a basis for customization. For details, see the comments in mpcCustomConstraint.m. The constraints in the custom constraint function are used in addition to linear constraints defined in the MPC controller.

To create an economic MPC controller, enable the traditional MPC controller to use custom cost and constraint functions.

mpcobj.Optimizer.CustomCostFcn = true;
mpcobj.Optimizer.CustomConstraintFcn = true;

Simulink Model with Economic MPC Controller

Open the Simulink model.

mdl = 'mpc_economicEO';
open_system(mdl)

C2H4 availability is initially 0.175 and remains constant for the first 500 seconds, which allows the controller to drive the plant to a new steady state with higher C2H4O production. C2H4 availability then steps down to 0.15 at t = 500 and remains constant for another 500 seconds. The economic MPC controller again maximizes C2H4O production at this condition.

The model includes constant (zero) references for the two plant outputs. The MPC Controller block requires these reference signals, but they are ignored in the custom cost function. The reactor temperature setpoint remains constant at its nominal value of 1.15.

The Plant subsystem calculates the six plant states by integrating the ODEs in oxidationPlantCT.m. Two measurements are fed back to the MPC Controller block as measured outputs. The C2H4O plant output is the instantaneous C2H4O production rate, which is used for display purposes.

Simulate Model and Analyze Results

Run the simulation.

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

In this example, you are able to calculate the steady-state $u_1$ and $u_2$ that maximize the C2H4O production rate (${y_1}{y_2}$) for a given C2H4 availability ($u_4$) and reactor temperature ($u_3$).

The true optimum steady states are:

  • When ${u_4} = 0.175$:

$$\begin{array}{l}
{u_1} = 0.205\\
{u_2} = 0.854\\
{y_1}{y_2} = 0.0137
\end{array}$$

  • When ${u_4} = 0.150$:

$$\begin{array}{l}
{u_1} = 0.245\\
{u_2} = 0.613\\
{y_1}{y_2} = 0.0116
\end{array}$$

Recall that at the initial condition:

$$\begin{array}{l}
{u_1} = 0.5\\
{u_2} = 0.469\\
{y_1}{y_2} = 0.0129
\end{array}$$

Therefore, the C2H4O plant operating at the nominal condition is not optimal, and its profit can be improved.

In the first 500 seconds, the economic MPC controller gradually moves the plant to a new steady state under the same C2H4 availability constraint:

$$\begin{array}{l}
{u_1} = 0.373\\
{u_2} = 0.469\\
{y_1}{y_2} = 0.0133
\end{array}$$

Therefore, the economic MPC controller improves C2H4O production rate by

$$\left( {0.0133 - 0.0129} \right)/0.0129 = 3.1\% $$

which could be worth millions of dollars per year in large-scale production.

The economic MPC controller does not drive the plant to the true optimum ${y_1}{y_2} = 0.0137$. This result is due to the contoller predicting the influence of $u_1$ and $u_2$ on $y_1$ and $y_2$ using the LTI model obtained at the nominal condition. The plant behavior described by the ODEs is nonlinear, and therefore, the controller predictions are not perfect.

You can improve controller performance by updating the LTI plant model at run time when conditions change. Doing so requires an adaptive MPC controller with the same custom cost and constraints. Even with adaptive control, it is unlikely that economic MPC would drive a physical plant to operate at the true optimum because MPC optimizes based on an internal model, which is always an imperfect representation of the real plant.

In the second 500 seconds, the C2H4 availability drops from 0.175 to 0.15. The economic MPC controller moves the plant smoothly to a new steady state:

$$\begin{array}{l}
{u_1} = 0.355\\
{u_2} = 0.423\\
{y_1}{y_2} = 0.01146
\end{array}$$

Again, the economic MPC controller brings the C2H4O production rate close to the true optimum 0.0116.

References

[1] H. Durand, M. Ellis, P. D. Christofides, "Economic model predictive control designs for input rate-of-change constraint handling and guaranteed economic performance", Computers and Chemical Engineering, Vol. 92, pp 18-36, 2016.

Remove example file folder from MATLAB path, and close Simulink model.

rmpath(fullfile(matlabroot,'examples','mpc_featured','main'));
bdclose(mdl)