% main_reacsim_pid.m
% Function to simulate PI control of a non-linear
% chemical reactor model using the explicit Euler method.
% Purpose: Supplement to Exercise 3

clear all
global u
% Nominal variables
u=25;
x=[2.5;1];

dt=0.0005;     % Sampling interval
tf=0.2;        % Final time instant
t=0:dt:tf;     % Time horizon
N=length(t);   % Number of discrete time instants

% PI controller parameters
Kp=62.5; Ti=1/125;
Kp=46.9; Ti=1/83.3; % SIMC tuning rules
Td=Ti/3;% Test

z=u;   % Initial controller state
r=1.0; % Step change in r from 1.0 to 1.02 at t=0

% Vectors to store the variables and allocate space
U=zeros(N,1); Y=zeros(N,1); X=zeros(N,2); R=zeros(N,1);

% The main simulation loop using explicit Euler method
y_old=x(2);
for k=1:N
    y=x(2);    % The measurement
    
    if k >= 41; r=1.02; end
    
    % The PI controller
    e=r-y;
    u=Kp*e+z - Kp*Td*(y-y_old)/dt; y_old=y;
    z=z+dt*Kp*e/Ti;    
       
    U(k,1)=u;  % Save u
    Y(k,1)=y;  % Save output y
    X(k,:)=x'; % Save the states
    R(k,1)=r;  % The reference r
    
    fx=fx_chemreac(0,x); % Calculate the right hand side in dot(x)=f(x,u). u enters as global.
    x=x+dt*fx;           % Calculate x_(k+1) using explicit Euler
end

% Plotting results
figure(2)
subplot(311), plot(t,U), grid, title('Control input u')
subplot(312), plot(t,X(:,1)), grid, title('State x_1')
subplot(313), plot(t,Y,'-r',t,R,'b'), grid, title('Output y=x_2 and reference r')
xlabel('Continuous time 0 \leq t \leq t_f')
% END main_reacsim_pid.m