% ov678_pi.m
% Goal: Simulation of system h_p(s)=exp(-2s)/(s^2+3s+2) with PI-controller
% h_c(s)=Kp(1+Tis)/(Tis) with different settings for Kp og Ti. 
% The PI-controller parameters are obtained using
% the "Skogestad/SIMC method" in Exercise 6, and parameters from Exercises 7 and 8.
% Functions used: dread.m
% Written by: David Di Ruscio, 5. november 2002.

clear all
num=[0,0,1];                        % Definerer teller og nevner i h_p(s)
den=[1,3,2];                        % uten transportforsinkelsesleddet.
[A,B,D,E]=tf2ss(num,den);           % Tilstandsrommodell uten transportforsinkelse.

t0=0; t1=20; N=200;                 % Lager en passende tidshorisont.
t=linspace(t0,t1,N);
dt=t(2)-t(1);                       % Samplingsintervall.

nt=floor(2/dt); yt=zeros(nt,1);     % Array for implementering av transport
                                    % forsinkelse.
r=1;                                % Referansesignal, sprang fra r=0 til r=1 ved t=t0.
% PI-regulatorparametre
tekst=['PI-controller parameters from Exercise 6, 1b : (choice 6)'
       'PI-controller parameters from Exercise 7.    : (choice 7)' 
       'Ziegler-Nichols PI-controller parameters.    : (choice 8)'];
disp(tekst);
pi_reg=6;
pi_reg=dread('PI-controller choice no:',pi_reg);
if pi_reg==6
   Kp=0.56; Ti=1.25;                   % PI-regulatorparametre fra oving 6, 1b.
elseif pi_reg==7 
   Kp=0.6631; Ti=1;                    % PI-regulatorparametre fra oving 7.
elseif pi_reg==8 
   Kp=1.4009; Ti=5.4366;               % Ziegler-Nichols PI-reg.parametre.
else
   disp('Warning: hvilke PI-reg. param. oensker du ?');
   Kp=0.56; Ti=1.25;
   disp('Skogestad parametre?');
end
x0=[0;0];                           % Initialverdi for tilstandsvektoren, x.
z0=0;                               % Initialverdi for regulatortilstanden, z.

%%% Simuleringsloekke for det lukkede systemet.
x=x0; z=z0;
for i=1:N
   y=D*x;                           % Måling før transportforsinkelseselementet.
   
   ym=yt(nt);                       % Implementering av transportforsinkelse.
   for k=nt:-1:2
      yt(k)=yt(k-1);
   end
   yt(1)=y;
   
   e=r-ym;                          % Beregner reguleringsavviket.
   
   u=Kp*e+z;                        % Implementering av PI-regulator.
   z=z+dt*Kp*e/Ti;                  % (integrerer med eksplisitt Euler.)
   
   Y(i,1:2)=[ym y];                 % Lagrer systemvariable.
   R(i,1)=r;
   U(i,1)=u;  
   
   fx=A*x+B*u;                      % Paatrykker paadraget paa prosessen.
   x=x+dt*fx;                       % (integrerer med eksplisitt Euler.)
end

%%% Plotting results
subplot(211), plot(t,U), grid
title('Simulation af closed loop system'), ylabel('u')
subplot(212), plot(t,[R Y(:,1)]), grid
xlabel('Time [s]'), ylabel('y and r')

print -deps ov6_fig1
