Accelerating the pace of engineering and science

# Documentation Center

• Trial Software

## Review Model Predictive Controller Design for Potential Stability and Robustness Issues

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

The 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.

Inputs:

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

Outputs:

```1. High Heating Value (HHV)
2. Wobbe Index (WI)
3. Flame Speed Index (FSI)

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

The plant is nonlinear but the basis for the controller design is the following linear predictive model in the state space format, which applies at a typical steady state operating point. Time units are 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);
end
```

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);
end
```

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);
end
```

Specify weights on manipulated variables.

MV weights are specified based on the relative cost of the blended streams. This allows MPC to choose among the six in order to minimize the cost of the blended fuel gas.

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

Specify weights on manipulated variable increments.

They are relatively small to allow freedom of movement.

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

Specify weights on plant outputs.

The OV weights reflect the relative magnitudes of the output variable signals, considering that the predictive model is not scaled. Since the nominal value of the header pressure (the 4th output) is relatively large, set its weight smaller to compensate.

```Obj.Weights.OV = [82, 65, 57, 16];
```

Using the review Command to Improve the Initial Design

Review the initial controller design.

```review(Obj)
```

The summary table shown above lists two warnings. The first warning is that the controller won't drive the output variables 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.0156 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.0469 in HHV, etc.

Since there are six MVs and only four OVs, enough degrees of freedom are available in control system and you might be surprised to see non-zero steady state offset. The root cause is the non-zero MV weights we use to force the plant toward the most economical operating condition. The non-zero steady state offsets are often undesirable but they are acceptable in this case because the primary objective is to minimize the blend cost and the gas quality is allowed to float as long as it stays within the specified limits, even when a process disturbance is present Small offset gain magnitudes indicate that the impact of disturbance would be limited. Moreover, small violations of those limits would be acceptable as well because output constraints are soft.

Review 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's good practice to eliminate the possibility by softening one of the two constraints. Since the MV minimum and maximum values are respected in this example, we soften the increment bounds as follows:

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

Review the new controller design.

```review(Obj)
```

The MV constraint conflict warning has disappeared.

Diagnosing the Impact of Zero Output Weights

Assuming that a new design requirement allows the controlled outputs vary freely within their limits, consider eliminating the weights on them:

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

Review the impact of this design change.

```review(Obj)
```

A new warning regarding QP Hessian Matrix Validity appears, in addition to the steady-state gain warning. Click QP Hessian Matrix Validity warning to see the details.

We see that the review has flagged the zero weights on all four output variables. But since the zero weights come from design requirement and the other Hessian tests indicate that the quadratic programming problem should have a unique solution, this warning can be safely 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.

Reviewing 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 gives 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, use

```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.

```review(Obj)
```

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 default 1e5 to a higher value to harden all the soft OV constraints.

```Obj.Weights.ECR = 1e8;
```

Review the impact of this design change.

```review(Obj)
```

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

```mpcverbosity(MPC_verbosity);
[~, hWebBrowser] = web;
close(hWebBrowser);
```