MATLAB Examples

Nonlinear MPC Control of an Ethylene Oxidation Plant

This example shows how to use a nonlinear model predictive controller (nonlinear MPC) to control an ethylene oxidation plant as it transitions from one operating point to another.

Contents

Product Requirement

This example requires Optimization Toolbox™ software to solve a nonlinear programming problem at each control interval.

if ~mpcchecktoolboxinstalled('optim')
    disp('Optimization Toolbox is required to run this example.')
    return
end

This example also 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

Add example file folder to MATLAB® path.

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

Ethylene Oxidation Plant

Oxidation 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

A mixture of air and ethylene is fed continuously. The plant is described by a first-principle nonlinear dynamic model, implemented as a set of ordinary differential equations (ODEs) in the oxidationCT function. For more information, see oxidationCT.m.

The plant contains four states:

  • Gas density in the reactor (x1)
  • C2H4 concentration in the reactor (x2)
  • C2H4O concentration in the reactor (x3)
  • Temperature in the reactor (x4)

The plant has two inputs:

  • Total volumetric feed flow rate (u1)
  • C2H4 concentration of the feed (u2)

The plant has one output:

  • C2H4O concentration in the effluent flow (y, equivalent to x3)

For convenience, all variables in the model are pre-scaled to be dimensionless. All the states are measurable as well. The plant equations and parameters are obtained from [1].

Control Objectives

In this example, the total volumetric feed flow rate (u1) is the manipulated variable (MV) and C2H4O concentration in the effluent flow (y) is the controlled variable (OV). Good tracking performance of y is required within an operating range from 0.03 to 0.05. The corresponding u1 values are 0.38 and 0.15.

The C2H4 concentration in the feed flow (u2) is a measured disturbance. Its nominal value is 0.5, and a typical disturbance has a size of 0.1. In addition to good reference tracking, the controller is required to reject such a disturbance.

An unmeasured step-like disturbance might occur in the C2H4O concentration in the effluent flow with a typical size of 0.01. The controller is required to reject such a disturbance as well.

The manipulated variable u1 has a range from 0.0704 to 0.7042 due to actuator limitations.

Plant Nonlinearity Assessment

A single linear MPC is often capable of controlling a nonlinear plant near its nominal operating point as long as the plant nonlinearity is not strong. One way to assess plant nonlinearity is to compare linear plant characteristics obtained from different operating points.

In this example, two linear plant models are obtained at operating points with low (y = 0.03) and high (y = 0.05) C2H4O concentration in the reactor, respectively.

Find the equilibrium operating point for low C2H4O concentration.

options = optimoptions('fsolve','Display','none');
uLow = [0.38;0.5];
xLow = fsolve(@(x) oxidationCT(x,uLow),[1 0.3 0.03 1],options);
[~,yLow] = oxidationCT(xLow,uLow);

Linearize the low concentration model by numerically perturbing the plant states and inputs.

plantLow = linearizeODEs(@oxidationCT,xLow,uLow);

Find the equilibrium operating point for high C2H4O concentration.

uHigh = [0.15;0.5];
xHigh = fsolve(@(x) oxidationCT(x,uHigh),[1 0.3 0.03 1],options);
[~,yHigh] = oxidationCT(xHigh,uHigh);

Linearize the high concentration model by numerically perturbing the plant states and inputs.

plantHigh = linearizeODEs(@oxidationCT,xHigh,uHigh);

Compare the frequency responses of the plants.

figure
bode(plantLow(:,1),plantHigh(:,1))
legend('low C2H4O','high C2H4O')

Some nonlinearity can be observed from the plot. For example, the DC gain of the plant at high C2H4O concentration is more than twice the DC gain of the plant at low C2H4O concentration. Therefore, if a controller is designed at one operating point and used at the other, there may be some performance and robustness issues.

Linear MPC Design

Design a linear MPC at the low C2H4O concentration.

Specify signal names and define MV and MD channels.

plantLow.InputName = {'Qin','CEin'};
plantLow.StateName = {'Den','CE','CEO','Tc'};
plantLow.OutputName = {'CEOout'};
plantLow = setmpcsignals(plantLow,'MV',1,'MD',2);

Create the controller, specifying sample time and horizons.

Ts = 5;
PredictionHorizon = 10;
ControlHorizon = 3;
mpcobj = mpc(plantLow,Ts,PredictionHorizon,ControlHorizon);
-->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 nominal plant conditions for the low concentration operating point.

mpcobj.Model.Nominal.X = xLow;
mpcobj.Model.Nominal.U = uLow;
mpcobj.Model.Nominal.Y = yLow;

Specify input and output scale factors using the plant nominal conditions.

mpcobj.OV.ScaleFactor = yLow;
mpcobj.MV.ScaleFactor = uLow(1);
mpcobj.DV.ScaleFactor = uLow(2);

Specify MV constraints based on the controller actuator limitations.

u1LowLimit = 0.0704;
u1HighLimit = 0.7042;
mpcobj.MV.Min = u1LowLimit;
mpcobj.MV.Max = u1HighLimit;

To achieve good robustness, specify an MV rate weight.

mpcobj.Weight.MVRate = 0.2;

To assess the robustness of the linear MPC at both low and high C2H4O concentration, first compute the open-loop transfer function, L = plant * controller, at each operating condition. Use an unconstrained MPC as the nominal controller for linear analysis.

L_Low = -ss(mpcobj)*c2d(plantLow(:,1),Ts);
L_High = -ss(mpcobj)*c2d(plantHigh(:,1),Ts);
-->Converting model to discrete time.
-->Assuming output disturbance added to measured output channel #1 is integrated white noise.
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Plot the open-loop frequency responses, displaying the minimum stability margins.

figure
h = bodeplot(L_Low,L_High);
showCharacteristic(h,'MinimumStabilityMargins')
legend('low C2H4O','high C2H4O','Location','SouthWest')

The phase margin of L_low is 64 degrees, while the phase margin of L_High is less than 14 degrees. The loss of robustness at high C2H4O concentration is expected because, as the plant gain becomes larger there, the open-loop bandwidth becomes higher with the same controller. Therefore, the phase margin decreases since phase compensation from the controller remains the same.

This decrease in robustness means that the linear MPC designed at low C2H4O concentration would most likely perform poorly at high C2H4O concentration.

Linear MPC Design Simulation

Validate the robustness assessment by simulating the system in Simulink.

Open the model.

mdlMPC = 'oxidationMPC';
open_system(mdlMPC)

Run the simulation and view the output.

sim(mdlMPC)
open_system([mdlMPC '/y'])
bdclose(mdlMPC)

In the simulation, a step setpoint change from 0.03 to 0.05 occurs at 100 seconds and a step measured disturbance of 0.1 occurs at 200 seconds. In both cases, the linear MPC exhibits a highly underdamped response with large overshoot, as expected. A step output disturbance of 0.01 occurs at 300 seconds and is poorly rejected as well.

Note that designing the linear MPC at a high C2H4O concentration would not help either. In that case, the controller would produce very sluggish tracking performance because the open-loop bandwidth would be too small.

Nonlinear MPC Design

When a single linear MPC does not work well, more advanced MPC techniques should be considered. One option is to use different linear MPC controllers at different operating points. Both adaptive MPC and gain scheduled MPC fall into this category. In this example, another option is explored, which is using nonlinear MPC with a nonlinear prediction model. This choice is reasonable in this example because:

  1. A nonlinear plant model is available.
  2. The controller sample time of 5 seconds is slow.

Because the plant model is nonlinear, the optimal control problem is converted into a nonlinear optimization problem with a nonlinear cost function and nonlinear constraints.

In this example, to compute the optimal control sequence, use the fmincon function from the Optimization Toolbox as the nonlinear programming solver. Select the sequential quadratic programming (SQP) algorithm, which is one of the most common approaches in nonlinear MPC applications.

The cost function used in this example is described in the oxidationFmincon function. It is similar to the standard cost function used by linear MPC, where both tracking error accumulated across the prediction horizon and aggressive manipulated variable moves are suppressed. For more information, see oxidationFmincon.m.

There are no nonlinear constraints in this case. Actuator limits are used as the limits of optimization variables in fmincon.

In this example, since all the states are measurable, full state feedback is used by the nonlinear MPC. In general, when there are unmeasurable states, you must design a nonlinear state estimator, such as an Extended Kalman Filter (EKF) or a Moving Horizon Estimator (MHE).

To reject step-like output disturbance, the measured plant output is also fed back to the controller. The difference between the measurement y and the output of prediction model at current time y(k) is added to each prediction from y(k+1) to y(k+P), such that the same cost function can be used to achieve zero steady state error. It effectively adds integral action to the controller.

To achieve a fair comparison, the same sample time, horizons, constraints, and weights are used in both the linear and nonlinear MPC designs. Reference previewing is also enabled in both cases.

Nonlinear MPC Simulation

To assess nonlinear MPC performance, simulate the system in Simulink.

Open the model.

mdlNMPC = 'oxidationNMPC';
open_system(mdlNMPC)

Run the simulation and view the output.

sim(mdlNMPC)
open_system([mdlNMPC '/y'])

The nonlinear MPC produces much better reference tracking and disturbance rejection performance, as expected. Since it uses a nonlinear prediction model, the nonlinear MPC can be used as a single controller for the whole operating range.

Note that although reference change occurs at t=100, MPC knows about the change as early as at t=50 because of reference previewing. Since the objective is to minimize tracking errors across whole horizon, MPC decides to move the plant in advance such that tracking error is the smallest across the prediction horizon. If previewing is disabled, MPC would start reacting at t=100 but it generates larger tracking error.

Note that fmincon does not support code generation.

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

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

Conclusion

This example illustrates a general workflow to design and simulate nonlinear MPC in MATLAB and Simulink. Depending on the specific nonlinear plant characteristics and control requirements, the implementation details can vary significantly. The key design challenges include:

  • Choosing an ODE solver
  • Designing a nonlinear state estimator
  • Choosing a nonlinear optimization solver
  • Designing the nonlinear cost function and constraint function

You can use the functions and Simulink model in this example as templates for other nonlinear MPC design and simulation tasks.

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.