MATLAB Examples

Review Model Predictive Controller for Stability and Robustness Issues

This example shows how to use the review command to detect potential issues with a model predictive controller design.


Fuel Gas Blending Process

The example application is a fuel gas blending process. The objective is to blend six gases to obtain a fuel gas, which is then burned to provide process heating. The fuel gas must satisfy three quality standards in order for it to burn reliably and with the expected heat output. The fuel gas header pressure must also be controlled. Thus, there are four controlled output variables. The manipulated variables are the six feed gas flow rates.


1. Natural Gas (NG)
2. Reformed Gas (RG)
3. Hydrogen (H2)
4. Nitrogen (N2)
5. Tail Gas 1 (T1)
6. Tail Gas 2 (T2)


1. High Heating Value (HHV)
2. Wobbe Index (WI)
3. Flame Speed Index (FSI)
4. Header Pressure (P)

The fuel gas blending process was studied by Muller et al.: "Modeling, validation, and control of an industrial fuel gas blending system", C.J. Muller, I.K. Craig, N.L. Ricker, J. of Process Control, in press, 2011.

Linear Plant Model

Use the following linear plant model as the prediction model for the controller. This state-space model, applicable at a typical steady-state operating point, uses the time unit of hours.

a = diag([-28.6120, -28.6822, -28.5134,  -0.0281, -23.2191, -23.4266, ...
          -22.9377, - 0.0101, -26.4877, -26.7950, -27.2210,  -0.0083, ...
          -23.0890, -23.0062, -22.9349,  -0.0115, -25.8581, -25.6939, ...
          -27.0793,  -0.0117, -22.8975, -22.8233, -21.1142,  -0.0065]);
b = zeros(24,6);
b( 1: 4,1) = [4, 4, 8, 32]';
b( 5: 8,2) = [2, 2, 4, 32]';
b( 9:12,3) = [2, 2, 4, 32]';
b(13:16,4) = [4, 4, 8, 32]';
b(17:20,5) = [2, 2, 4, 32]';
b(21:24,6) = [1, 2, 1, 32]';
c = [diag([ 6.1510,  7.6785, -5.9312, 34.2689]), ...
     diag([-2.2158, -3.1204,  2.6220, 35.3561]), ...
     diag([-2.5223,  1.1480,  7.8136, 35.0376]), ...
     diag([-3.3187, -7.6067, -6.2755, 34.8720]), ...
     diag([-1.6583, -2.0249,  2.5584, 34.7881]), ...
     diag([-1.6807, -1.2217,  1.0492, 35.0297])];
d = zeros(4,6);
Plant = ss(a, b, c, d);

By default, all the plant inputs are manipulated variables.

Plant.InputName = {'NG', 'RG', 'H2', 'N2', 'T1', 'T2'};

By default, all the plant outputs are measured outputs.

Plant.OutputName = {'HHV', 'WI', 'FSI', 'P'};

Transport delay is added to plant outputs to reflect the delay in the sensors.

Plant.OutputDelay = [0.00556  0.0167  0.00556  0];

Initial Controller Design

Construct an initial model predictive controller based on design requirements.

Specify sampling time, horizons and steady-state values.

The sampling time is that of the sensors (20 seconds). The prediction horizon is approximately equal to the plant settling time (39 intervals). The control horizon uses four blocked moves that have lengths of 2, 6, 12 and 19 intervals respectively. The nominal operating conditions are non-zero. The output measurement noise is white noise with magnitude of 0.001.

MPC_verbosity = mpcverbosity('off'); % Disable MPC message displaying at command line
Ts = 20/3600;   % Time units are hours.
Obj = mpc(Plant, Ts, 39, [2, 6, 12, 19]);
Obj.Model.Noise = ss(0.001*eye(4));
Obj.Model.Nominal.Y = [16.5, 25, 43.8, 2100];
Obj.Model.Nominal.U = [1.4170, 0, 2, 0, 0, 26.5829];

Specify lower and upper bounds on manipulated variables.

Since all the manipulated variables are flow rates of gas streams, their lower bounds are zero. All the MV constraints are hard (MinECR and MaxECR = 0) by default.

MVmin = zeros(1,6);
MVmax = [15, 20, 5, 5, 30, 30];
for i = 1:6
    Obj.MV(i).Min = MVmin(i);
    Obj.MV(i).Max = MVmax(i);

Specify lower and upper bounds on manipulated variable increments.

The bounds are set large enough to allow full range of movement in one interval. All the MV rate constraints are hard (RateMinECR and RateMaxECR = 0) by default.

for i = 1:6
    Obj.MV(i).RateMin = -MVmax(i);
    Obj.MV(i).RateMax =  MVmax(i);

Specify lower and upper bounds on plant outputs.

All the OV constraints are soft (MinECR and MaxECR = 0) by default.

OVmin = [16.5, 25, 39, 2000];
OVmax = [18.0, 27, 46, 2200];
for i = 1:4
    Obj.OV(i).Min = OVmin(i);
    Obj.OV(i).Max = OVmax(i);

Specify weights on manipulated variables.

MV weights are specified based on the known costs of each feed stream. This tells MPC controller how to move the six manipulated variables in order to minimize the cost of the blended fuel gas. The weights are normalized so the maximum weight is approximately 1.0.

Obj.Weights.MV = [54.9, 20.5, 0, 5.73, 0, 0]/55;

Specify weights on manipulated variable increments.

They are small relative to the maximum MV weight so the MVs are free to vary.

Obj.Weights.MVrate = 0.1*ones(1,6);

Specify weights on plant outputs.

The OV weights penalize deviations from specified setpoints and would normally be "large" relative to the other weights. Let us first consider the default values, which equal the maximum MV weight specified above.

Obj.Weights.OV = [1, 1, 1, 1];

Improve the Initial Design

Review the initial controller design.


The summary table shown above lists three warnings and one error. Consider these in turn. Click QP Hessian Matrix Validity and scroll down to display the warning It indicates that the plant signal magnitudes differ significantly. Specifically, the pressure response is much larger than the others.

Examination of the specified OV bounds shows that the spans are quite different, and the pressure span is two orders of magnitude larger than the others. It is good practice to specify MPC scale factors to account for the expected differences in signal magnitudes. We are already weighting MVs based on relative cost, so we focus on the OVs only.

Calculate OV spans.

OVspan = OVmax - OVmin;

Use these as the specified scale factors.

for i = 1:4
    Obj.OV(i).ScaleFactor = OVspan(i);

Verify that the scale factor warning has disappeared.


The next warning indicates that the controller does not drive the OVs to their targets at steady state. Click Closed-Loop Steady-State Gains to see a list of the non-zero gains.

The first entry in the list shows that adding a sustained disturbance of unit magnitude to the HHV output would cause the HHV to deviate 0.0860 units from its steady-state target, assuming no constraints are active. The second entry shows that a unit disturbance in WI would cause a steady-state deviation ("offset") of -0.0345 in HHV, etc.

Since there are six MVs and only four OVs, excess degrees of freedom are available and you might be surprised to see nonzero steady-state offsets. The nonzero MV weights we have specified in order to drive the plant toward the most economical operating condition are causing this.

Nonzero steady-state offsets are often undesirable but are acceptable in this application because:

  1. The primary objective is to minimize the blend cost. The gas quality (HHV, etc.) can vary freely within the specified OV limits.
  2. The small offset gain magnitudes indicate that the impact of disturbances would be small.
  3. The OV limits are soft constraints. Small, short-term violations are acceptable.

View the second warning by clicking Hard MV Constraints, which indicates a potential hard-constraint conflict.

If an external event causes the NG to go far below its specified minimum, the constraint on its rate of increase might make it impossible to return the NG within bounds in one interval. In other words, when you specify both MV.Min and MV.RateMax, the controller would not be able to find an optimal solution if the most recent MV value is less than (MV.Min - MV.RateMax). Similarly, there is a potential conflict when you specify both MV.Max and MV.RateMin.

An MV constraint conflict would be extremely unlikely in the gas blending application, but it is good practice to eliminate the possibility by softening one of the two constraints. Since the MV minimum and maximum values are physical limits and the increment bounds are not, we soften the increment bounds as follows:

for i = 1:6
    Obj.MV(i).RateMinECR = 0.1;
    Obj.MV(i).RateMaxECR = 0.1;

Review the new controller design.


The MV constraint conflict warning has disappeared.

Click Soft Constraints to view the error message.

We see that the delay in the WI output makes it impossible to satisfy bounds on that variable until at least three control intervals have elapsed. The WI bounds are soft but it is poor practice to include unattainable constraints in a design. We therefore modify the WI bound specifications such that it is unconstained until the 4th prediction horizon step.

Obj.OV(2).Min = [-Inf(1,3), OVmin(2)];
Obj.OV(2).Max = [ Inf(1,3), OVmax(2)];

Rerunning the review command verifies that this eliminates the error message (see the next step).

Diagnose Impact of Zero Output Weights

Given that the design requirements allow the OVs to vary freely within their limits, consider zeroing their penalty weights:

Obj.Weights.OV = zeros(1,4);

Review the impact of this design change.


A new warning regarding QP Hessian Matrix Validity has appeared.

Click QP Hessian Matrix Validity warning to see the details.

The review has flagged the zero weights on all four output variables. Since the zero weights are consistent with the design requirement and the other Hessian tests indicate that the quadratic programming problem has a unique solution, this warning can be ignored.

Click Closed-Loop Steady-State Gains to see the second warning. It shows another consequence of setting the four OV weights to zero. When an OV is not penalized by a weight, any output disturbance added to it will be ignored, passing through with no attenuation.

Since it is a design requirement, non-zero steady-state offsets are acceptable provided that MPC is able to hold all the OVs within their specified bounds. It is therefore a good idea to examine how easily the soft OV constraints can be violated when disturbances are present.

Review Soft Constraints

Click Soft Constraints to see a list of soft constraints - in this case an upper and lower bound on each OV.

The Impact Factor column shows that using the default MinECR and MaxECR values give the pressure (P) a much higher priority than the other OVs. If we want the priorities to be more comparable, we should increase the pressure constraint ECR values and adjust the others too. For example, consider:

Obj.OV(1).MinECR = 0.5;
Obj.OV(1).MaxECR = 0.5;
Obj.OV(3).MinECR = 3;
Obj.OV(3).MaxECR = 3;
Obj.OV(4).MinECR = 80;
Obj.OV(4).MaxECR = 80;

Review the impact of this design change.


Notice from the Sensitivity Ratio column that all the sensitivity ratios are now less than unity. This means that the soft constraints will receive less attention than other terms in the MPC objective function, such as deviations of the MVs from their target values. Thus, it is likely that an output constraint violation would occur.

In order to give the output constraints higher priority than other MPC objectives, increase the Weights.ECR parameter from the default, 1e5, to a higher value to harden all the soft OV constraints.

Obj.Weights.ECR = 1e8;

Review the impact of this design change.


The controller is now a factor of 100 more sensitive to output constraint violations than to errors in target tracking.

Review Data Memory Size

Click Memory Size for MPC Data to see the estimated memory size needed to store the MPC data matrices used on the hardware.

In this example, if the controller is running using single precision, it requires 250 KB of memory to store its matrices. If the controller memory size exceeds the memory available on the target system, you must redesign the controller to reduce its memory requirements. Alternatively, increase the memory available on the target system.

[~, hWebBrowser] = web;