clear;
% True plant & true initial state
% Initial state
e1=0;e2=0;e3=0;e4=0;e5=0; %  Fractional error in gains

% Uses kij  to represent gains, tij to represent time constants and dij to represent delays
% u1=top draw, u2=side draw, u3=botttoms duty, u4= reflux duty, u5= upper reflux duty
% y1=top end point, y2=side end point, y3= bottoms reflux temp
% gij represents transfer function relating input i and output j.
%


% Defining the transfer function of the process

delt=0;

num=4.05+2.11*e1;den=[50 1];delay=27;  
sys11=ss(tf(num,den,'InputDelay',delay));
num=5.39+3.29*e1;den=[50 1];delay=18;  
sys12=ss(tf(num,den,'InputDelay',delay));

num=3.66+2.29*e1;den=[33 1];delay=20;  
sys13=ss(tf(num,den,'InputDelay',delay));

num=1.77+.39*e2;den=[60 1];delay=28;  
sys21=ss(tf(num,den,'InputDelay',delay));

num=5.72+.57*e2;den=[60 1];delay=14;  
sys22=ss(tf(num,den,'InputDelay',delay));

num=4.42+3.11*e2;den=[44 1];delay=22;  
sys23=ss(tf(num,den,'InputDelay',delay));

num=5.88+.59*e3;den=[50 1];delay=27;  
sys31=ss(tf(num,den,'InputDelay',delay));

num=6.90+.89*e3;den=[40 1];delay=15;  
sys32=ss(tf(num,den,'InputDelay',delay));

num=7.2+1.33*e3;den=[19 1];delay=0;  
sys33=ss(tf(num,den,'InputDelay',delay));

num=1.2+.12*e4;den=[45 1];delay=27;  
sys41=ss(tf(num,den,'InputDelay',delay));

num=1.52+.13*e4;den=[25 1];delay=15;  
sys42=ss(tf(num,den,'InputDelay',delay));

num=1.14+.18*e4;den=[27 1];delay=0;  
sys43=ss(tf(num,den,'InputDelay',delay));

num=1.44+.16*e5;den=[40 1];delay=27;  
sys51=ss(tf(num,den,'InputDelay',delay));

num=1.83+.13*e5;den=[20 1];delay=15;  
sys52=ss(tf(num,den,'InputDelay',delay));

num=1.26+.18*e5;den=[32 1];delay=0;  
sys53=ss(tf(num,den,'InputDelay',delay));

sys=[sys11 sys21 sys31 ;sys12 sys22 sys32;sys13 sys23 sys33];

Ts=5;     %Sample time
Tstop=300;  %Simulation time
model=c2d(sys,Ts); % prediction model
p=30; % output horizon
moves=4; % Control move horizon

umin=-[.5 .5 .5]; umax=[.5 .5 .5 ]; dumin=-[.25 .25 .25]; 
dumax=[.25 .25 .25];
ymin=[-.5 -Inf -.5]; ymax=[.5  Inf Inf]; % output limits

limits={umin,umax,dumin,dumax,ymin,ymax};
 uweights=[0 0 0];duweights=[1 1 1];yweights=[1 1 1];  
 weights={uweights ,duweights,yweights};
clear MPCadd;MPCadd=[];
open_system('ch18ex2.mdl')
sim('ch18ex2',Tstop);
subplot(2,1,1);plot(time,[y1 y2 y3]);
subplot(2,1,2);plot(time,[u1 u2 u3]);