%SIMEVAP Simulate evaporator under predictive control, with MPC Toolbox
%
% This is a script file which simulates the Newell and Lee evaporator
% running under predictive control, using the MPC Toolbox.
% It was used to produce the results shown in Chapter 9.2 of
% `PREDICTIVE CONTROL WITH CONSTRAINTS' by J.M.Maciejowski.
%
% It requires the Simulink model 'evaporator.mdl' and the function
% 'scmpcnl2.m', both of which are available on the book's web site.
% It also requires the MPC Toolbox for Matlab.
%
% Edit this file to run different scenarios.

% J.M. Maciejowski, 9 April 2000
% Cambridge University Engineering Dept.
% Copyright (C) 2000. All Rights Reserved.

%%%%% WRITE INITIAL INFO TO SCREEN:

disp('Simulation of Newell and Lee Evaporator')
disp(' ')
disp('EVALSIM is a script file. Edit the file to change')
disp('        parameters and scenarios.')
disp(' ')

%%%%% DEFINE VARIABLES NEEDED BY EVAPORATOR MODEL:

lambda = 38.5; % kW/(kg/min)
UA2 = 6.84;    % kW/K
Cp = 0.07;     % kW/K(kg/min)
lambdas = 36.6;% kW/(kg/min)
rhoA = 20;     % kg/m
% Sample time:
Ts = 1;  % minute - units of time are minutes throughout.
disp(['Sample time = ',num2str(Ts),' minutes']);

%%%%% DEFINE INITIAL CONDITIONS FOR EVAPORATOR:

% The states are ordered as follows:
% (use '[sizes,x0,xstring]=evaporator' to get this information)
%    STATE              EQUILIBRIUM VALUE
%    Tank Level L2        1.0   m
%    Composition X2      25.0   percent
%    Pressure P2         50.5   kPa
%    P100 Lag           194.7   kPa
%    F3 Lag              50.0   kg/min
%    F200 Lag           208.0   kg/min
%    F2 Lag               2.0   kg/min

x0 = [1.0; 25.0; 50.5; 194.7; 50.0; 208.0; 2.0]; % Initial state vector

%%%%% OBTAIN INITIAL LINEARISED INTERNAL MODEL:

% The inputs are ordered as follows:
% (ordering corresponds to Inport numbering)
% INPUT                 EQUILIBRIUM VALUE
%   Control inputs:
% Product flowrate F2               2.0   kg/min
% Steam pressure P100 demand      194.7   kPa
% Cooling water rate F200 demand  208.0   kg/min
%   Disturbance inputs:
% Circulating rate F3 demand       50.0   kg/min
% Feed flowrate F1                 10.0   kg/min
% Feed composition X1               5.0   percent
% Feed temperature T1              40.0   deg C
% Cooling water inlet temp T200    25.0   deg C

% Define input values at operating point:
u0 = [2.0; 194.7; 208.0; 50.0; 10.0; 5.0; 40.0; 25.0];

% Linearise:
[A,B,C,D] = linmod2('evaporator',x0,u0); % Continuous-time model

%Retain only 3 first outputs ( = L2, X2, P2):
C = C(1:3,:);
D = D(1:3,:);

% Discretise model:
[Ad,Bd,Cd,Dd] = ssdata(c2d(ss(A,B,C,D),Ts));

% Put model into MPC 'MOD' format:
n = size(Ad,1); % Number of states
nu = 3;  % Number of control inputs
nd = 0;  % Number of measured disturances
nw = 5;  % Number of unmeasured disturbances
nym = 3; % Number of measured outputs
nyu = 0; % Number of unmeasured outputs
minfo = [Ts,n,nu,nd,nw,nym,nyu];
imod = ss2mod(Ad,Bd,Cd,Dd,minfo); % Linear internal model

%%%%% DEFINE CONSTRAINTS ON EVAPORATOR:

% Input amplitude constraints:
%   (Min: 0.  Max: 2*nominal operating value approx.)
%ulim = [0,0,0,0,4,400,400,100];  % in MPC Toolbox format
ulim = [0,0,0,4,400,400];

% Input rate constraints: Assume no constraints, since we already
% have lags to represent local controllers for F2, P100, F200, F3.
%ulim = [ulim,1e4,1e4,1e4,1e4]; % Set rate limits to very high values
ulim = [ulim,1e4,1e4,1e4];

%Output amplitude constraints:
%  L2: 0 to 2m,  X2: 0 to 100 percent,  P2: 0 to 100 kPa.
ylim = [0,0,0,2,100,100];  % in MPC Toolbox format

%%%%% DEFINE PREDICTIVE CONTROL PARAMETERS:

P = 30;  % Prediction horizon (steps)
M = 3;   % Control horizon (steps)

ywt = [1000,100,100];   % Weight on tracking error penalty
%uwt = [1,1,1,1];    % Weight on control move penalty
uwt = [1,1,1];

% Set-point trajectory: 
% Ramp X2 from 25 percent to 15 percent over 20 minutes,
% Ramp P2 from 50.5 kPa to 70 kPa over 20 minutes,
% try to hold L2 at 1.0:
X2ramp = [25*ones(1,20), 25:(15-25)/20:15];
P2ramp = [50.5*ones(1,20), 50:(70-50)/20:70];
nsp = length(X2ramp);
setpoints = [1.0*ones(nsp,1), X2ramp', P2ramp'];

Kest = [];  xm0 = []; % Default observer
% No noise or measured disturbances, or unmeasured disturbances on u:
z = []; d = []; wu = [];

% Unmeasured disturbance: +- 10% changes on F1:
F1dist = [zeros(10,1);2*ones(10,1);-ones(10,1);zeros(10,1);ones(10,1);-ones(10,1)];
nd = length(F1dist);
%w = [zeros(nd,1),F1dist,zeros(nd,1),zeros(nd,1),zeros(nd,1)];
w = [];

% Hard constraints:
suwt = [];  sywt = []; normtype = [];
% Soft constraints:
%suwt = 1e6; sywt = []; normtype = 'inf-norm';

disp(['Prediction horizon = ',int2str(P),' steps,   ',...
     'Control horizon = ',int2str(M),' steps'])
disp(['Output weights : ',num2str(ywt),...
     '.   Input weights : ',num2str(uwt)])
if isempty(normtype),
   disp('Hard constraints')
else
   disp(['Soft constraints, norm = ',normtype])
end

%%%%% DEFINE SIMULATION PARAMETERS:

% Start and stop times:
tstart = 0;   tend = 100;  % minutes
disp(['End time = ',num2str(tend)])

%%%%% RUN SIMULATION:

disp(' ')
disp('To avoid re-linearisation enter interval > end time (Inf ok)')
badinput = 1;
while badinput,
   relin = input('Input interval between re-linearisations: ');
   if relin == 1,
      disp('*** Software limitation: Cannot re-linearise correctly at every step.')
      disp('*** Specify interval > 1')
   else
      badinput = 0;
   end
end

if relin >= tend,  % No re-linearisation requested
   
   disp('Running simulation without re-linearsation')
   tvec = [tstart,tend];
   [y,u]=scmpcnl2('evaporator',imod,ywt,uwt,M,P,tvec, ...
                 setpoints,ulim,ylim,Kest,z,d,w,wu, ...
                 x0,u0(1:3),[],xm0,suwt,sywt,normtype);
              
else
              
   disp(['Re-linearising every ',num2str(relin),' minutes (simulated time)'])
   % Set up first segment:
   tvec = [tstart,relin-1]; % max() handles case relin=1
   x0seg = x0;  u0seg = u0(1:3);
   if tvec(2)+1>nsp,
        spseg = [setpoints(tvec(1)+1:nsp,:); setpoints(nsp,:)];
   else
        spseg = setpoints(tvec(1)+1:tvec(2)+1,:); 
   end
   y = []; u = [];
   
   while (tvec(1)+1<tend),
     [yseg,useg,ymseg,xseg]=scmpcnl2('evaporator',imod,ywt,uwt,M,P,tvec, ...
                 spseg,ulim,ylim,Kest,z,d,w,wu, ...
                 x0seg,u0seg,[],xm0,suwt,sywt,normtype);
   
     seglength = size(xseg,1);
     x0seg = xseg(seglength,:)';  % Initial state for next segment
     u0seg = useg(seglength,:)';  % Initial input for next segment
     y = [y;yseg];  u = [u;useg]; % Write results into y and u
     tvec = tvec + relin; % New start and stop times
     tvec(2) = min(tvec(2),tend); % Don't simulate after tend
     % Set-points for next segment:
     if tvec(2)+1>nsp,
        spseg = [setpoints(tvec(1)+1:nsp,:); setpoints(nsp,:)];
     else
        spseg = setpoints(tvec(1)+1:tvec(2)+1,:); 
     end
     % Re-linearise at current state:
     [A,B,C,D] = linmod2('evaporator',x0seg,[u0seg;u0(4:8)]); % Continuous-time model
     %Retain only 3 first outputs ( = L2, X2, P2):
     C = C(1:3,:);     D = D(1:3,:);
     [Ad,Bd,Cd,Dd] = ssdata(c2d(ss(A,B,C,D),Ts)); % Discrete-time
     imod = ss2mod(Ad,Bd,Cd,Dd,minfo); % MPC format
   end
   
end
              
%%%%% PROCESS RESULTS:

tvector = (0:(size(y,1)-1))'*Ts;
% Pad setpoints with final values or truncate as necessary:
nt = length(tvector);
if nt>nsp,
   lastsp = setpoints(nsp,:);
   setpoints = [setpoints; ones(nt-nsp,1)*lastsp];
elseif nt<nsp,
   setpoints = setpoints(1:nt,:);
end

fhandle=figure;

subplot(321)
set(gca,'FontSize',8);
plot(tvector,y(:,1),'-',tvector,setpoints(:,1),'--'); grid;
ylabel('L2 (m)','FontSize',8),title('Separator Level','FontSize',8);

subplot(322)
set(gca,'FontSize',8);
plot(tvector,y(:,2),'-',tvector,setpoints(:,2),'--'); grid;
ylabel('X2 (%)','FontSize',8),title('Product Concentration','FontSize',8);

subplot(323)
set(gca,'FontSize',8);
plot(tvector,y(:,3),'-',tvector,setpoints(:,3),'--'); grid;
ylabel('P2 (kPa)','FontSize',8),title('Operating Pressure','FontSize',8);

subplot(324)
set(gca,'FontSize',8);
stairs(tvector,u(:,1),'-'); grid;
ylabel('F2 (kg/min)','FontSize',8),title('Product Flow Rate SP','FontSize',8);

subplot(325)
set(gca,'FontSize',8);
stairs(tvector,u(:,2),'-'); grid;
xlabel('Time (minutes)','FontSize',8)
ylabel('P100 (kPa)','FontSize',8),title('Steam Pressure SP','FontSize',8);

subplot(326)
set(gca,'FontSize',8);
stairs(tvector,u(:,3),'-'); grid;
xlabel('Time (minutes)')
ylabel('F200 (kg/min)','FontSize',8),title('Water Flow Rate SP','FontSize',8);

%%%%%%%%%%%%%%%%%%%%  End of script file SIMEVAP.M  %%%%%%%%%%%%%%%%%%%%%%%%%
