% DEMO_REACT.M
% Purpose: Illustrate conventional PID control of a chemical reaktor.
% DDIR 27.03.2001

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

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

%% Compute concentration of A, x1, at the inlet.
x10= (k1*x(1)+k3*x(1)^2)/u+x(1)   % Compute nominal inlet concentration, x10=C_Ai.

%% Check RHS of dot(x)=f(x,u) if f(x,u)=0.
f=zeros(2,1);
f(1)=-k1*x(1)-k3*x(1)^2+(x10-x(1))*u;
f(2)= k1*x(1)-k2*x(2)-x(2)*u;

h=0.0011; t1=0.1; t=0:h:t1; N=length(t); % Time axis data.

%% Get some parameters from the user
Text=['*** CONTROL OF CHEMICAL REACTOR ***'
      '1 : PI-control                     '
      '2 : Open loop step response in u   '
      '3 : PID control                    '];
disp(Text);   
ireg=1;
ireg=dread('Regulert (1) uregulert (2) ',ireg);
if (ireg==1) ~= (ireg==3)
   Kp=25; Ti=10;
   z=25;
   r=1.05;
   r=dread('Spesifiser referanse: r=',r);
   Kp=dread('Kp=',Kp);
   Ti=dread('Ti=',Ti);
   if ireg==3; Td=5; Td=dread('Td=',Td); end
elseif ireg == 2
   u=30;
   u=dread('Sprang i foedepaadrag: u=',u);
end

%% Main simulation loop
X=zeros(N,2); Y=zeros(N,1); U=zeros(N,1);
D=[0,1]; eold=r-D*x;
for i=1:N
   y=D*x;
   
   if ireg == 1                          % Discrete PI-controller
      e=r-y;
      u=Kp*e+z;
      z=z+Kp*e/Ti;
   elseif ireg == 2                      % Open loop step response in u.
      u=u;
   elseif ireg ==3
      e=r-y;
      u=Kp*e+z+Kp*Td*(e-eold);
      z=z+Kp*e/Ti;
      eold=e;
   end
   
   U(i,1)=u; Y(i,1)=y; X(i,:)=x';        % 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

subplot(211), plot(t,X(:,1)), grid, ylabel('x_1')
title('Chemical reactor responses')
subplot(212), plot(t,X(:,2)), grid, ylabel('x_2')
xlabel('Continuous time: 0 \leq t \leq 0.1')