% Function: ex3b_half_pi_du.m
% Simulating integrator+delay process with PI-controller
% on velocity (deviation) form,
% and with SIMC (the Skogestad), and some alternative settings
% for the PI controller parameters
% Process model: y=h_p(s)(u+v), h_p(s)=K*exp(-tau*s)/s, K=1, tau=1

% Written: November 15th 2011, David Di Ruscio

clear all

v=0;             % Disturbance
r=1;             % reference
kg=1; tau=1;     % slope and time delay, hp=kg*exp(-tau*s)/s;

dt=0.1;          % Sampling time
t=0:dt:80; t=t'; % Simulation time interval
N=length(t);     % Number of discrete time instants

nt=tau/dt;       % Size of array to implement time delay
Delay=zeros(nt,1); % arrays to make time delay, hold (nt) samples

itun=2;
if itun == 1
   Kp=1/(2*kg*tau); Ti=8*tau;          % Skogestad setting
elseif itun == 3
   Kp=1/(2*kg*tau); Ti=5.83*tau;       % Settings in this note
elseif itun == 2
   Kp=1/(2*kg*tau); Ti=5*tau           % Buttherworth setting
end

x=0; z=-v;
Delay=Delay+x;
ireg=2;
y_old=0; g1=-Kp; g2=dt*Kp/Ti; u=-v;  % Initializations for the velocity formulation
for k=1:N
    yp=x;
    %%% Implementation of time delay, using array method.
    % programming of time delay for process model, 
    % hold (y) in array of length (nt) samples
    y=Delay(nt); % the delayed process output Y=y(i-nt)
    for id=nt:-1:2;
      Delay(id)=Delay(id-1);
    end
    Delay(1)=yp;
    %%% End implementation of time delay
    
    if k>floor(40/dt); v=1; end      % Disturbance step at time, t=40
    
    % PI controller
    if ireg==1         % Standard discrete implementation of PI controller
       e=r-y;
       u=Kp*e+z;
       z=z+dt*Kp*e/Ti;
    elseif ireg==2     % Velocity/deviation formulation of PI controller
        u=u+g1*(y-y_old)+g2*(r-y_old);
        y_old=y;       % save: y_old=y(k-1)
    end
    
    % Save for plotting puposes
    Y(k,1)=y;
    R(k,1)=r;
    
    % Process update
    x=x+dt*kg*(u+v);
end

plot(t,Y,'-',t,R,'-.'), grid
legend('Output, y','Reference, r')
xlabel('Continuous time:  0 <= t <= 80')