%% Use Custom Constraints in Blending Process
% This example shows how to design an MPC controller for a blending process 
% using custom mixed input/output constraints.

% Copyright 2015 The MathWorks, Inc.

%% 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\limits_{i = 1}^3 {{\phi _i}}  - \phi \\
% V\frac{{d{\gamma _j}}}{{d\tau }} = \sum\limits_{i = 1}^3 
% {\left( {{\gamma _{ij}} - {\gamma _j}} \right){\phi _i}} 
% \end{array}$$
%
% where
%
% * $V$ is the mixture inventory (in the container).
% * $\phi_i$ is the plow rate for feed $i$.
% * $\phi$ is the rate at which the blend is being removed from inventory,
% that is the demand.
% * $\gamma _{ij}$ is the concentration of constituent $j$ in feed $i$.
% * $\gamma _j$ is the concentration of constituent $j$ in the blend.
% * $\tau$ is time.
%
%
% In this example, there are two important constituents, $j$ = 1 and 2.
%
% The control objectives are targets for the two constituent concentrations
% in the blend and the mixture inventory. The challenge is that the demand,
% $\phi$, and feed compositions, $\gamma _{ij}$, vary. The inventory, blend
% compositions, and demand are measured, but the feed compositions are
% unmeasured.
%
% At the nominal operating condition:
%
% * Feed 1, $\phi _1$, (mostly constituent 1) is 80% of the total inflow.
% * Feed 2, $\phi _2$, (mostly constituent 2) is 20%.
% * Feed 3, $\phi _3$, (pure constituent 1) is not used.
%
% The process design allows manipulation of the total feed entering the 
% mixing chamber, $\phi _T$, and the individual rates of feeds 2 and 3. In 
% other words, the rate of feed 1 is:
%
% $$\phi _1=\phi _T-\phi _2-\phi _2$$
%
% Each feed has limited availability:
%
% $$0 \le {\phi _i} \le {\phi _{i,\max }}$$
%
% The equations are normalized such that, at the nominal steady state, the 
% mean residence time in the mixing container is $\tau=1$.
%
% The constraint $\phi _{1,\max }=0.8$ is imposed by an upstream process,
% and the constraints $\phi _{2,\max }=\phi _{3,\max }=0.6$ are imposed by
% physical limits.

%% Define Linear Plant Model
% The blending process is mildly nonlinear, however 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.

%%
% Specify the number of feeds, |ni|, and the number of constituents, |nc|.
ni = 3;
nc = 2;

%%
% Specify the nominal flow rates for the three input streams and the output
% stream, or demand. At the nominal operating condition, the output flow
% rate is equal to the sum of the input flow rates.
Fin_nom = [1.6,0.4,0];
F_nom  = sum(Fin_nom);

%%
% Define the nominal constituent compositions for the input feeds, where 
% |cin_nom(i,j)| represents the composition of constituent |i| in feed |j|.
cin_nom = [0.7 0.2 0.8;0.3 0.8 0];

%%
% Define the nominal constituent compositions in the output feed.
cout_nom = cin_nom*Fin_nom'/F_nom;

%%
% 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 a state-space model with feed flows |F1|, |F2|, and |F3| as MVs:
A = [zeros(1,nc+1); zeros(nc,1) -eye(nc)];
Bu = [ones(1,ni); gij-1];

%%
% Change the MV definition to [FT, F2, F3] where F1 = FT - F2 - F3
Bu = [Bu(:,1), Bu(:,2)-Bu(:,1), Bu(:,3)-Bu(:,1)];

%%
% Add the measured disturbance, blend demand, as the 4th model input.
Bv = [-1; zeros(nc,1)];
B = [Bu Bv];

%%
% Define all of the states as measurable. The states consist of the mixture
% inventory and the constituent concentrations.
C = eye(nc+1);

%%
% Specify that there is no direct feed-through from the inputs to the
% outputs.
D = zeros(nc+1,ni+1);

%%
% Construct the linear 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'};

%% Create MPC Controller
% Specify the sample time, prediction horizon, and control horizon for the
% controller.
Ts = 0.1;   
p = 10;
m = 3;

%%
% Create the controller.
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
% for all outputs.
mpcobj.Model.Nominal.Y = [1 1 1]; 

%%
% Specify the normalized nominal values for the manipulated variables, |u(1)|,
% |u(2)| and |u(3)|, and the measured disturbance, |u(4)|.
mpcobj.Model.Nominal.U = [1 fin_nom(2) fin_nom(3) 1];

%%
% Specify output tuning weights. To pay more attention to controlling the
% inventory and the composition of the first constituent, use larger
% weights for the first two outputs.
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.

%% Specify Mixed Constraints
% Given the specified upper bounds on the feed 2 and 3 rates (0.6), it is
% possible that their sum could be as much as 1.2. Since the nominal total
% feed rate is 1.0, the controller can request a physically impossible
% condition, where the sum of feeds 2 and 3 exceeds the total feed rate,
% which implies a negative feed 1 rate.
%
% The following constraint prevents the controller from requesting an 
% unrealistic $\phi _1$ value.
%
% $$0 \le {\phi _1} = {\phi _T} - {\phi _2} - {\phi _3} \le 0.8$$
%
% Specify this constraint in the form $Eu + Fy \le g$.
E = [-1 1 1; 1 -1 -1];
g = [0;0.8];

%%
% Since no outputs are specified in the mixed constraints, set their
% coefficients to zero.
F = zeros(2,3);

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

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

%%
% Set the custom constraints in the MPC controller.
setconstraint(mpcobj,E,F,g,v,h)

%% Simulate Model in Simulink
% The Simulink model contains a nonlinear model of the blending process and
% an unmeasured disturbance in the constituent 1 feed composition.
%
% The |Demand|, $\phi$, is modeled as a measured disturbance. The operator 
% can vary the demand value, and the resulting signal goes to both the 
% process and the controller.
%
% The model simulates the following scenario:
% 
% * At $\tau=0$, the process is operating at steady state.
% * At $\tau=1$, the |Total Demand| decreases from $\phi=1.0$ to $\phi=0.9$.
% * At $\tau=2$, there is a large step increase in the concentration of
% constituent 1 in feed 1, from 1.17 to 2.17.

%%
% Open and simulate the Simulink model.
mdl = 'mpc_blendingprocess';
open_system(mdl)
sim(mdl)

%%
% In the simulation:
%
% * At time 0, the plant operates steadily at the nominal
% conditions.
% * At time 1, the demand decreases by 10%, and the controller maintains
% the inventory close to its setpoint.
% * At time 2, there is a large unmeasured increase in the concentration of
% constituent 1 contained in feed 1. This disturbance causes a prediction
% error and a large disturbance in the blend composition.
%
% The disturbance is a nonlinear effect, but the linear MPC controller
% recovers well and drives the blend composition back to its setpoint

%% Verify Effect of Custom Constraints
% Plot the feed rate signals.
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 unmeasured disturbance occurs at time 2, which requires the
% controller to decrease F1. During the transient, F1 becomes zero. If the
% mixed input/output constraint had not been included, F1 would have been
% negative. The controller requests for FT, F2, and F3 would have been
% impossible to satisfy, which would lead to a performance degradation.
% With the constraint included, the controller does its best given the
% physical limits of the system.

bdclose(mdl)
