%% Nonlinear Blending Process with Custom Constraints
% This example shows how to use constraints on a combination of inputs and
% outputs to control a nonlinear blending process.

% Copyright 1990-2014 The MathWorks, Inc.

%% About the Blending Process 
% A continuous blending process combines three streams in a well-mixed
% vessel to produce a blend with desired properties. These properties
% depend on concentrations of two key constituents:
%
% * Feed rates
%
% * Inventory of materials in the blending vessel.
% 
% Let the feed rates be |F1|, |F2|, and |F3|, and the total |FT = F1 + F2 + F3|.
% The plant configuration allows |F2|, |F3|, and |FT| to be controlled
% independently.  Therefore, F1 = FT - F2 - F3 is determined by difference.
%
% The control objective is to adjust the feed rates in order to keep the
% concentrations in the vessel near target values. The MPC controller must
% also control the total inventory of material in the blending vessel.

%%
% The process operates under the following constraints:
%
% * Physical limitations require |0 <= FT <= 2| and |0 <= F2, F3 <= 0.6|.
%
% * Each feed has limited availability. For |F2| and |F3|, maximum
% availability equals the physical limit 0.6. For F1, a maximum
% availability of 0.8 is imposed by an upstream process.

%%
% During operation, it is important to prevent the controller from
% requesting |FT < F2 + F3|, which implies |F1 < 0| which is impossible.
% Similarly, |FT > F2 + F3 + 0.8|, implying |F1 > 0.8| should also be
% avoided. This can be accomplished using MPC controller with a mixed
% constraint (constraints on input and output combinations).

%%
% System disturbances include the feed stream compositions and the demand
% for blended material, |F|, defined as the rate at which material leaves
% the blending vessel.  This is set by a downstream process.  From the
% point of view of the blending process, |F| is a measured disturbance.
% Feed composition disturbances are unmeasured.
% 
% Measurements available for feedback control include the constituent
% concentrations in the vessel and the blend inventory.

%% Define Linear Model for the Blending Process 
% 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'};

%% Design 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

%% Specify Mixed Constraints
% Impose the following constraints on the control system:
%
% * Physical constraint |F1 = FT - F2 - F3 >= 0|.
% 
% * Availability constraint, |F1 <= 0.8|.  
%
% Putting these in the standard form, you obtain:
%
% |-FT + F2 + F3 <= 0|
% |FT - F2 - F3 <= 0.8|
%
% This defines the input constraint coefficient matrix:
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);

%% Simulate Using Simulink(R)
% To run this example, Simulink(R) is required.
if ~mpcchecktoolboxinstalled('simulink')
    disp('Simulink(R) is required to run this example.')
    return
end
%%
% Open and simulate the Simulink(R) model. The model includes a nonlinear
% model of the blending process. The time scale is normalized such
% that one time unit is the vessel's nominal mean residence time. 
mdl = 'mpc_blendingprocess';
open_system(mdl);
sim(mdl);
%%
% The simulation shows:
%
% * At time 0, the plant is operating steadily at the nominal conditions.
%
% * At time 1, the demand decreases 10% (from 1 to 0.9).
%
% * At time 2, there is a large (unmeasured) increase in the concentration
% of constituent 1 contained in feed 1.  This is a nonlinear effect but the
% linear MPC compensates well (target values are 1 for the three controlled
% outputs).


%% Validate No Mixed Constraint Violations
% Plot the input and output signals to check if the mixed constraint is 
% violated during simulation:
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:
%
% * When demand drops, the composition disturbance requires |F1| to be
% decreased.  During the transient, |F1| becomes zero. If the mixed
% constraint had not been included, it would have gone negative, i.e., the
% controller's requests for |FT|, |F2|, and |F3| would have been impossible
% to satisfy, causing a performance degradation.  With the constraint
% included, the controller does its best given the physical limits.
% 
% * The availability constraint |F1 <= 0.8| is maintained at all times.

%%
bdclose(mdl)
