Accelerating the pace of engineering and science

# Documentation

## Using Custom Constraints in Blending Process

A continuous blending process combines three feeds in a well-mixed container to produce a blend having desired properties. The dimensionless governing equations are:

$\begin{array}{l}\frac{dv}{d\tau }=\sum _{i=1}^{3}{\varphi }_{i}-\varphi \\ V\frac{d{\gamma }_{j}}{d\tau }=\sum _{i=1}^{3}\left({\gamma }_{ij}-{\gamma }_{j}\right){\varphi }_{i}\end{array}$

where V is the mixture inventory (in the container), ϕi is the flow rate of the ith feed, ϕ is the demand, i.e., the rate at which the blend is being removed from inventory, γij is the concentration of constituent j in feed i, γj is the concentration of j in the blend, and τ is time. In this example, there are two important constituents, j = 1 and 2.

The control objectives are targets for the blend's two constituent concentrations and the mixture inventory. The challenge is that the demand ϕ and feed compositions γij vary. The inventory, blend compositions, and demand are measured, but the feed compositions are unmeasured.

At the nominal operating condition:

• Feed 1 ϕ1 (mostly constituent 1) is 80% of the total inflow ϕ.

• Feed 2 ϕ2 (mostly constituent 2) is 20%.

• Feed 3 ϕ3 (pure constituent 1) is not used.

The process design allows manipulation of the total feed entering the mixing chamber and the individual rates of feeds 2 and 3. In other words, the rate of feed 1 is:

ϕ1 = ϕT – ϕ2 – ϕ3

Each feed has limited availability:

0 ≤ ϕi ≤ ϕi,max

The equations are normalized such that—at the nominal steady state—the mean residence time in the mixing container is τ = 1. The target inventory is V = 1, and the target blend composition is γ1 = γ2 = 1.

The constraints ϕi,max = 0.8 is imposed by an upstream process and ϕ2,max = ϕ3,max = 0.6 is imposed by physical limits.

### MPC Controller with Custom Input/Output Constraints

1. Open the Simulink® model mpc_blendingprocess:

In the model, an MPC controller controls the blending process. The block labeled Blending incorporates the previously described model equations and includes unmeasured (step) disturbances in the feed compositions.

Signal u represents controller adjustments, i.e.,

u = [ϕT ϕ2 ϕ3].

Signal v represents the demand, ϕ which is a measured disturbance. The operator can vary the demand, and the resulting signal goes to the process and controller.

Consider the following scenario:

• At τ = 0, the process is operating at its nominal steady state.

• At τ = 0.5, the demand decreases from ϕ = 1 to ϕ = 0.9.

• At τ = 1, there is a large increase in the concentration of constituent 1 in feed 1 γ11 from 1.17 to 2.17.

The plant is mildly nonlinear. You can derive a linear model at the nominal steady state. This approach is quite accurate unless the (unmeasured) feed compositions change. If the change is sufficiently large, the steady-state gains of the nonlinear process change sign and the closed-loop system can become unstable.

2. Define a linear model.

```% Create a linear approximation -- a state-space model based on the nominal
% operating point:
ni = 3;                     % number of feed streams
nc = 2;                     % number of components
Fin_nom = [1.6, 0.4, 0];    % Nominal flow rate for the ith feed stream
F_nom  = sum(Fin_nom);      % Nominal flow rate for the exit stream (demand)
cin_nom = [0.7 0.2 0.8      % Nominal composition for jth constituent in the ith feed flow
0.3 0.8   0];
cout_nom = cin_nom*Fin_nom'/F_nom;  % Nominal product composition

% Normalize the linear model such that the target demand is 1 and the
% product composition is 1:
fin_nom = Fin_nom/F_nom;
gij = [cin_nom(1,:)/cout_nom(1)
cin_nom(2,:)/cout_nom(2)];

% Create the state-space model with feed flows |[F1, F2, F3]| as MVs:
A = [ zeros(1,nc+1)
zeros(nc,1) -eye(nc)];
Bu = [ones(1,ni)
gij-1];
% Change MV definition to [FT, F2, F3] where F1 = FT - F2 - F3
Bu = [Bu(:,1), Bu(:,2)-Bu(:,1), Bu(:,3)-Bu(:,1)];
% Add the blend demand as the 4th model input, a measured disturbance
Bv = [-1
zeros(nc,1)];
B = [Bu Bv];
% All the states (inventory and compositions) are measurable
C = eye(nc+1);
% No direct feed-through term
D = zeros(nc+1,ni+1);
% Construct the plant model
Model = ss(A, B, C, D);
Model.InputName = {'F_T','F_2','F_3','F'};
Model.InputGroup.MV = 1:3;
Model.InputGroup.MD = 4;
Model.OutputName = {'V','c_1','c_2'};```
3. Specify an MPC controller.

```% Create the controller object with sampling period, prediction and control
% horizons:
Ts = 0.1;
p=10;
m=3;
MPCobj = mpc(Model, Ts, p, m);

% The outputs are the inventory |y(1)| and the constituent concentrations
% |y(2)| and |y(3)|.  Specify nominal values of unity after normalization:
MPCobj.Model.Nominal.Y = [1 1 1];

% The manipulated variables are |u1 = FT|, |u2 = F2|, |u3 = F3|.  Specify
% nominal values after normalization:
MPCobj.Model.Nominal.U = [1 fin_nom(2) fin_nom(3) 1];

% Specify output tuning weights. Larger weights are assigned to the first
% two outputs because we pay more attention to controlling the inventory
% and composition of the first blending material:
MPCobj.Weights.OV = [1 1 0.5];

% Specify the hard bounds (physical limits) on the manipulated variables:
umin = [0 0 0];
umax = [2 0.6 0.6];
for i = 1:3
MPCobj.MV(i).Min = umin(i);
MPCobj.MV(i).Max = umax(i);
MPCobj.MV(i).RateMin = -0.1;
MPCobj.MV(i).RateMax =  0.1;
end```

The total feed rate and the rates of feed 2 and feed 3 have upper bounds. Feed 1 also has an upper bound, determined by the upstream unit supplying it. Under normal conditions, the plant operates far from these bounds but for the scenario outlined previously, the controller must reduce the rate of feed 1 drastically, as it is bringing in excess constituent 1. To do this, the controller must increase the rates of feeds 2 and 3 (keeping the total feed rate close to the demand rate to maintain the target inventory.)

4. Specify constraints.

Given the specified bounds on the feed 2 and 3 rates (= 0.6), it is possible that their sum could be as much as 1.2. Because the total feed rate is of order 0.9 to 1.0, the controller can request a physically impossible condition in which the sum of feeds 2 and 3 exceeds the total feed rate. This implies a negative feed 1 rate.

The constraint

0 ≤ ϕ1 = ϕT – ϕ2 – ϕ3 ≤ 0.8

prevents the controller from requesting an unrealistic ϕ1 value.

To specify this constraint, enter:

```E = [-1 1 1; 1 -1 -1];

% No outputs are specified in the mixed constraints, so set their
% coefficients to zero:
F = zeros(2,3);

% Specify vector g in E*u + F*y <= g:
g = [0; 0.8];

% Specify that both constraints are hard (ECR = 0):
v = zeros(2,1);

% Specify zero coefficients for the measured disturbance:
h = zeros(2,1);

% Include the mixed constraints in the controller object:
setconstraint(MPCobj, E, F, g, v, h);
```

v = zeros(2,1) defines hard constraints, which are reasonable here because the constraints involve manipulated variables only. If the constraints involved a mixture of input and output variables, use soft constraints.

5. Simulate the model and plot the input and output signals.

```sim('mpc_blendingprocess')
figure
plot(MVs.time,[MVs.signals(1).values(:,2), ...
(MVs.signals(2).values + MVs.signals(3).values), ...
(MVs.signals(1).values(:,2)-MVs.signals(2).values-MVs.signals(3).values)])
grid
legend('FT','F2+F3','F1')```

The plot shows the evolution of the total feed rate (blue curve) and the sum of feeds 2 and 3 (green curve). They coincide between τ = 1.7 and τ = 2.2.

If the custom input constraints had not been included, the controller would have requested a negative feed 1 rate during this period, as shown in by the red curve.

The controller maintains the inventory very close to its setpoint, but the severe disturbance in the feed composition causes a prediction error and a large disturbance in the blend composition (especially for constituent 1). Despite this, the controller recovers and drives the blend composition back to its setpoint, as shown in the following output of the CVs scope.