clc
clg
echo on
% demo943 program to demonstrate constrained single objective mpc of
% the Shell column
% Copyright 1992 by Babu Joseph
%  Washington University
%
% First define some parameters
% to change these parameters you must edit this file
pause % press any key
clc;
ts=5;  %  This is the sample time
p=30;pp=p  %  this is the output horizon
n=2;   % this is the control move horizon
yr1=0.5 ;%  this is the set point for y1 : the top end point
yr2=0;%  this is the set point for y2 : the side end point
% y3 is the bottom reflux temperature this is a constrained variable
nt=40; %  this is the no of time steps for simulation
l1d=0;% the disturbance in l1   input
l2d=0;% the disturbance in l2 input
beta=4;% move suppression factor
u3set=-.5; % set point for u3 this is taken as the desired target value of u3
y3min=-0.5;% the minimum value allowed for y3 this is an output constraint
pause  % press any key to continue
clc
% the model is generated using another file called demo73 this must be
% in your directory for this program to run
% first generate a state space model of the column from transfer functions
% this program generates a transfer function and state space model
% consisting of five inputs and three outputs
% the inputs are
%         u1: top draw rate
%         u2  intermediate draw rate
%         u3  bottom duty
%         u4 or l1 the top reflux duty, a disturbance variable
%         u5 or l2 the intermediate reflux duty a disturbance variable
%  the outputs are
%          y1   top purity
%          y2   intermdiate purity
%          y3   bottom reflux temperature
%
% this may take a few seconds
pause  % press any key
clc
ch18ex2a
p=pp;
% convert to discrete model
[ad,bd,cd,dd]=c2dm(ac,bc,cc,dc,ts);
% generate dynamic matrices using another  program  demo76
ch18ex2b
echo on
pause % press any key
clc
pack
% initialize vectors
time=0;
% predicted values of y
yp1=zeros(p,1);yp2=zeros(p,1);yp3=zeros(p,1); yp=zeros(3*p,1);
% measured values of y
ym1=0;ym2=0;ym3=0;
% disturbance vectors
d1=zeros(p,1);d2=zeros(p,1);d3=zeros(p,1);d=zeros(3*p,1);
dyp1=zeros(p,1);dyp2=zeros(p,1);dyp3=zeros(p,1);dyp=zeros(3*p,1);
% set point vectors
yr1=ones(p,1) .*yr1;yr2=ones(p,1) .*yr2 ; yr=[yr1 ;yr2];
[na,ma]=size(ad);
% state vector for simulation
x0=zeros(na,1);
% input vectors
% these are two long for the simulation denotes current and value at next time
u1=[ 0 0 ]';u2=[0  0]';u3=[0 0]';u4=[l1d l1d]';u5=[l2d l2d]';
pu1=zeros(1,nt);pu2=pu1;pu3=pu1;py1=pu1;py2=pu1;py3=pu1;tpl=pu1;
% this vector stores values of y3 for constraint calculation
y3min=ones(p,1) .*(y3min);
B4=eye(n,n); for i=1:n ;for j=1:n;if j<=i B4(i,j) =1 ;end;end;end;
In1=ones(n,1);
pause  % press any key
clc;
echo off;
disp('u1      u2       u3       y1       y2       y3 ')
% start iterations
for ii=1:nt;
time=time+ts;
% compute disturbance
dd1=ym1-yp1(1); dd2=ym2-yp2(1);dd3=ym3-yp3(1);
for j=1:p d1(j)=dd1;d2(j)=dd2;d3(j)=dd3; end;
d=[d1;d2;d3];
% compute control vector   using quadratic program
% set up matrices first
B1=[A11 A12 A13];B2=[A21 A22 A23];
k1=yp1+d1-yr1;k2=yp2+d2-yr2;
I3n=eye(3*n,3*n); I3nc=ones(3*n,1);
k4=In1 .*(u3(1)-u3set);
znn=zeros(n,n);B44=[znn znn znn;znn znn znn;znn znn B4'*B4] .*5;
znn1=zeros(n,1);k44=[znn1;znn1;B4'*k4]' .*5;
H=B1'*B1+B2'*B2+I3n .*beta+B44  ;
C=(k1'*B1+k2'*B2+k44) ;
B3=-[A31 A32 A33];
k3=-y3min+yp3 +d3;
A=[I3n;B3];
b=[I3nc .*0.2 ;k3];
% call qp to solve quadratic program
du=qp(H,C,A,b);
% recover control action at current time
du1=du(1);du2=du(n+1);du3=du(2*n+1);
u1=u1+[du1 du1]';u2=u2+[du2 du2]';u3=u3+[du3 du3]';
u=[u1 u2 u3 u4 u5];
% do simulation for one step to compute the next measurement
[y,x]=dlsim(ad,bd,cd,dd,u,x0);
ym1=y(2,1);ym2=y(2,2);ym3=y(2,3);
% save state for next time step
x0=x(2,:);
% update predictions
%
dyp1=[0 a11'] .*du1+[0 a12'] .* du2 + [0 a13'] .* du3;
dyp2=[0 a21'] .*du1+[0 a22'] .* du2 + [0 a23'] .* du3;
dyp3=[0 a31'] .*du1+[0 a32'] .* du2 + [0 a33'] .* du3;
dyp1=dyp1(:,1:p)';dyp2=dyp2(:,1:p)';dyp3=dyp3(:,1:p)';
% update
yp1=yp1+dyp1+d1;yp2=yp2+dyp2+d2;yp3=yp3+dyp3+d3;
yp1=yp1(2:p);yp2=yp2(2:p);yp3=yp3(2:p);
yp1(p)=yp1(p-1);yp2(p)=yp2(p-1);yp3(p)=yp3(p-1);
% save variables for plotting
tpl(ii+1)=time;pu1(ii+1)=u1(1);pu2(ii+1)=u2(1);pu3(ii+1)=u3(1);
py1(ii+1)=ym1;py2(ii+1)=ym2;py3(ii+1)=ym3;
% display intermediate  results on screen
fprintf( '% 7.2f % 7.3f % 7.3f  ',u1(1),u2(1),u3(1));
fprintf ('% 7.3f % 7.3f % 7.3f\n ',ym1,ym2,ym3)
end % end of time iterations
echo on
clg;clc
% the results will now be plotted
% the upper graph shows the behavior of y1(---),y2(+++) and y3( ***)
% the lower graph show the behavior of u1(---) u2(+++) and u3 (***)
% check graph for constraint violations
pause % press any key
subplot(211),plot(tpl,py1,'-',tpl,py2,'+',tpl,py3,'*'); ...
title('output response  y1:---- y2:++++ y3: ****');...
subplot(212),plot(tpl,pu1,'-',tpl,pu2,'+',tpl,pu3,'*');...
title('input variables  u1:---- u2: ++++ u3: ****');pause