% DLQ_EX3_DU.M
% Purpose: Illustrate mpc of a chemical reaktor.
% DDIR 29.10.01

Text=['*** CONTROL OF CHEMICAL REACTOR ***        '
      '1 : LQ optimal control with integral action'
      '2 : PI-control                             '];
disp(Text);   
ireg=1;
ireg=dread('LQ-control (1) PI-control (2) ',ireg);

%% Reaction coefficients.
k1=50; k2=100; k3=10;

%% Nominal process variables.
us=25;                                  % Steady state feed flowrate.
xs=zeros(2,1);
xs(1)=2.5;                              % Steady state concentration of A.
xs(2)=1;                                % Steady state concentration of B.
x10=10;
r=1;                                    % reference for x_2.

%% Linearized continuous and discrete model
h=0.002; t1=0.3; t=0:h:t1; N=length(t); % Time axis data.

Ac=[-125,0;50,-125]; Bc=[7.5;-1]; Dc=[0,1]; Ec=0;   % Continuous model.
[A,B,D,E]=c2dm(Ac,Bc,Dc,Ec,h,'zoh');                % Discrete model.

%% Specified parameters.
r=1;                                                % Nominal reference for y.
p=1; q=500;                                         % Weighting parameters.
q=dread('q=',q);

[G1i,G2i]=dlqdu_pi(A,B,D,q,p);

%% Make augmented state space model for LQ-design.
nx=size(A,1); nu=size(B,2); ny=size(D,1);
At=[A,zeros(nx,ny);D,eye(ny,ny)];Bt=[B;zeros(ny,nu)];Dt=[D,eye(ny,ny)];
P=p; 
Qt=q*Dt'*Dt;

%% Solve Riccati-equation and compute feedback matrix.
[K,Rr]=dlqr(At,Bt,Qt,P);
G=-K; G1=G(:,1:nx); G2=G(:,nx+1:nx+ny);

%% PI-controller parameters
Kp=46.9; Ti=1/83; Ti=Ti/h;

%% Main simulation loop
X=zeros(N,2); Y=zeros(N,1); U=zeros(N,1); R=zeros(N,1); f=zeros(2,1);
x=xs; u=us; y_old=D*xs; x_old=xs;
for k=1:N
   y=D*x;
   
   if k>14; r=0.7; end
   if k>84; r=1.0; end
   
   if ireg==1
    u=u+G1*(x-x_old)+G2*(y_old-r);
    y_old=y; x_old=x;
   elseif ireg==2
    u=u+(-Kp)*(y-y_old)+(-Kp/Ti)*(y_old-r);
    y_old=y;
   end 
    
   U(k,1)=u; Y(k,1)=y; X(k,:)=x'; R(k)=r;    % Save variables.
   
   f(1)=-k1*x(1)-k3*x(1)^2+(x10-x(1))*u;     % Update process model.
   f(2)= k1*x(1)-k2*x(2)-x(2)*u;   
   x=x+h*f;
end

if ireg==1
   subplot(311), plot(t,U), grid, ylabel('u')
   title('Chemical reactor responses: MPC-deviation control')
   subplot(312), plot(t,X(:,1)), grid, ylabel('x_1')
   subplot(313), plot(t,Y,'-',t,R,'-.'), grid, 
   legend('y_k','r_k'), ylabel('y=x_2')
   xlabel('Continuous time: 0 \leq t \leq 0.3')
elseif ireg==2
   subplot(311), plot(t,U), grid, ylabel('u')
   title('Chemical reactor responses: conventional PI-control')
   subplot(312), plot(t,X(:,1)), grid, ylabel('x_1')
   subplot(313), plot(t,Y,'-',t,R,'-.'), grid, 
   legend('y_k','r_k'), ylabel('y=x_2')
   xlabel('Continuous time: 0 \leq t \leq 0.3')

end