
% main_kolreg_lq
% Simulation program for discrete LQ optimal PI control of a distillation column.
% DDIR november 2001. Comments translated November 7 2011
% Dependent functions: dlqdu_pi.m
clear
%path(path,'s:\tex3\home\dlqdu_paper\matlab')

ireg=2;
ireg=dread('LQ-optimal PI (2) or PI (1)',ireg);

u=[2;2.5];                                       % Steady state control values, u_1=R og u_2=V.
v=[1;0.5];                                       % v_1=F og v_2=x_F
r=[0.04;0.96];                                   % References for bottom and top products.

h=5.0;                                           % Sampling time.
t=0:h:250;                                       % Simulation horizon.
N=length(t);
X=zeros(N,8); Y=zeros(N,2); U=zeros(N,2); T=zeros(N,1);

x0=(0.2:0.1:0.9)'
x0=ssprof2(u,v,2);                               % Compute steady state composition profile, solve f(x,u,v)=0 for x

D=zeros(2,8); D(1,1)=1; D(2,8)=1;                % D matrix in measurement equation.

if ireg==2                                       % LQ-optimal PI controller.
  q=2500; q=dread('q=',q);
  [Ac,Bc,Cc,Dc]=fx2abcd(x0,u,v);                 % Matrices in linearized continuous state space model
  [A,B,D,E]=c2dm(Ac,Bc,Dc,zeros(2,2),h,'zoh')    % Linear discrete time model
  P=eye(2); Q=eye(2)*q;                          % Chosing weighting matrices P and Q
  [G1,G2]=dlqdu_pi(A,B,D,Q,P);                   % Feedback matrices in the MIMO PI controller
  
  %nx=size(A,1); nu=size(B,2); ny=size(D,1);
  %G1s=zeros(nu,nx); G1s(1,end)=G1(1,end); G1s(end,1)=G1(end,1); G1=G1s;
  %G2s=zeros(nu,ny); G2s(1,end)=G2(1,end); G2s(end,1)=G2(end,1); G2=G2s;
 
elseif ireg==1                                  % Conventional  PI controller parameters
  Kp1=-20.3; Ti1=6.12; 
  Kp2=239*0.2; Ti2=22*0.5;
  G1=[0,Kp1;Kp2,0]; G2=[0,Kp1/Ti1;Kp2/Ti2,0];
  z=-inv(G2)*u;
end


Tall=[]; Xall=[]; Uall=[]; Rall=[];
x=x0; x_old=x; y_old=D*x;                     % Initialising
for k=1:N-1
   y=D*x;                                     % Take measurements first..
   
   if t(k)>50; r(2)=0.97; end                 % Perform some changes in the references
   if t(k)>150; r(1)=0.05; end
   
   if ireg==1                                 % Conventional PI-controller
     % Diskret PI regulator
     %u=-G1*(r-D*x)+G2*z
     %z=z+h*(y-r);
   
     u=-G1*(r-y)-G2*z
     z=z+h*(r-y);
   
   elseif ireg==2                             % LQ-optimal MIMO PI-controller
     u=u+G1*(x-x_old)+G2*(y_old-r);
     y_old=y; x_old=x;
   end
   
   % Integrating the process model over the sampling intervall, i.e.. from time t(k) to time t(k)+h.
   X(k,:)=x'; Y(k,:)=y'; T(k,1)=t(k); U(k,:)=u';
   [Tk,Xk]=ode15s('fcol2',t(k):0.1:t(k+1),x,[],u,v);
   x=Xk(end,:)';
   Tall=[Tall;Tk]; Xall=[Xall;Xk];            % Lagring...
   nu=length(0:0.1:h);
   Uall=[Uall;[u(1)*ones(nu,1) u(2)*ones(nu,1)]];
   Rall=[Rall;[r(1)*ones(nu,1) r(2)*ones(nu,1)]];
end 
y=D*x; 
X(N,:)=x'; Y(N,:)=y'; T(N,1)=t(N); U(N,:)=u';

%%%%%%%%%%%%%% Plotting the results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
subplot(211), plot(Tall,[Xall(:,1) Rall(:,1)]), ylabel('x_1=x_B'), grid
title('Top and bottom compositions')
subplot(212), plot(Tall,[Xall(:,8) Rall(:,2)]), ylabel('x_2=x_D'), grid

figure(2)
subplot(211), plot(Tall,Uall(:,1)), ylabel('u_1=R'), grid
title('Column control inputs')
subplot(212), plot(Tall,Uall(:,2)), ylabel('u_2=V'), grid
