%NLH: Videregående modellering og regulering våren 2002
%Finn Haugen, 3.2.02
%Løsning til datamaskinøving 3

%---------------------------

clear

%Deloppgave 1a:

%Definerer prosessparametre:
K=1;
w0=0.5;
zeta=1;
tidsforsinkelse=1;
%Lager kontinuerlig prosessmodell i form av LTI-modell:
hpkont=tf([K*w0^2],[1, 2*zeta*w0,w0^2],'InputDelay',tidsforsinkelse)

%Lager diskret prosessmodell basert på diskretinsering med nullte ordens holdeelement
%(zero order hold = zoh) på prosessinngangen:
Ts=0.1; %Samplingsintervallet
hpdisk=c2d(hpkont,Ts,'zoh') 

%---------------------------

%Deloppgave 1b:

%Plotter sprangresponsen for både hpkont og hpdisk - de bør være (og er) nokså like:
figure(1)
step(hpkont,hpdisk)
%Kommantar til sprangresponsen vist i figur 1:
%Den er som forventet: 
% * Stasjonær respons lik 1
% * Intet oversving (egentlig på grensen til oversving da zeta = 1)
% * Responstid (63% stigetid) er ca 4 sek, som stemmer ok med anslaget 1,5/w0 = 1,5/0,5 = 3

%---------------------------

%Deloppgave 2:

%Finner først PI-regulatorens z-transferfunksjon ved å foreta
%substitusjonen s = (z-1)/(Ts*z) i den kontinuerlige regulatorens
%transferfunksjon: Kp*(1+Ti*s)/(Ti*s). Resultatet blir
%hrdisk(z) = (Kp*(Ts+Ti)*z-Kp)/(Ti*z-Ti)
%Vi kunne ha funnet denne transferfunksjonen også ved å
%diskretisere den kontinuerlige regulatorfunksjonen
%u = Kp*e + (Kp/Ti)*Integral_0^t e*dt
%vha. Eulers bakovermetode og så z-transformere den resulterende
%differenslikningen, som er 
%u(k)=u(k-1)+Kp*(1+Ts/Ti)*e(k)-Kp*e(k-1).

%Lager LTI-modell for hrdisk(z). Setter da Ti svært stor, slik at PI-regulatoren 
%blir lik en P-regulator (ihht. Ziegler-Nichols' metode).
Ti=100000000
%Setter også i utgangspunktet Kp=0 (ihht. Z-Ns metode):
Kp=4.5
%Lager LTI-modell for diskret PI-regulator:
hrdisk=tf([Kp*(Ts+Ti),-Kp*Ti],[Ti,-Ti],Ts)

%Trenger reguleringssystemets modell, følgeforholdet M(z):
htilbakekopling=tf(1,1);
M=feedback(series(hrdisk,hpdisk),htilbakekopling)

%Simulerer gjentatte ganger responsen i prosessutgangen y etter sprang i referansen r 
%for økende Kp inntil det oppstår stående svingninger i prosessutgangen:
figure(2)
tidslutt=50;
step(M,tidslutt)
%break 
%Finner fra simuleringene:
Kpkritisk = 4.5
Tkritisk = 6.7
%hvilket gir regulatorparametrene
KpZN=Kpkritisk*0.45 % = 2,0
TiZN=Tkritisk/1.2 % = 5,6

%---------------------------

%Deloppgave 3a:

%Sprangresponsen i y med ferdig innstilt regulator:
Kp=KpZN;
Ti=TiZN;
hrdisk=tf([Kp*(Ts+Ti),-Kp*Ti],[Ti,-Ti],Ts)
M=feedback(series(hrdisk,hpdisk),htilbakekopling)
figure(3)
tidslutt=50;
step(M,tidslutt)
%Reguleringssystemets sprangrespons er ok (se figur 3)!

%---------------------------

%Deloppgave 3b:

%Øker Kp:
Kp=KpZN*2;
Ti=TiZN;
hrdisk=tf([Kp*(Ts+Ti),-Kp*Ti],[Ti,-Ti],Ts)
M=feedback(series(hrdisk,hpdisk),htilbakekopling)
figure(4)
tidslutt=50;
step(M,tidslutt)
%Reguleringssystemet er ustabilt!

%---------------------------

%Deloppgave 4a:

%Velger Z-N-regulatorparameterverdier:
Kp=KpZN;
Ti=TiZN;
hrdisk=tf([Kp*(Ts+Ti),-Kp*Ti],[Ti,-Ti],Ts)
M=feedback(series(hrdisk,hpdisk),htilbakekopling)
%Beregner poler (og nullpunkter) for M(z):
[polerM,nullpunkterM]=pzmap(M)
%Beregner polenes avstander fra origo:
abs(polerM)
%Vi ser at alle polene ligger innenfor enhetssirkelen, 
%så reguleringssystemet er asymptotisk stabilt!

%Plotter polene:
figure(5)
pzmap(M)
%Joda, polene synes å ligge innenfor enhetssirkelen.

%---------------------------

%Deloppgave 4b:

%Som i 4a, men nå med fordoblet Kp:
Kp=KpZN*2;
Ti=TiZN;
hrdisk=tf([Kp*(Ts+Ti),-Kp*Ti],[Ti,-Ti],Ts)
M=feedback(series(hrdisk,hpdisk),htilbakekopling)
%Beregner poler (og nullpunkter) for M(z):
[polerM,nullpunkterM]=pzmap(M)
%Beregner polenes avstander fra origo:
abs(polerM)
%Vi ser to av polene har absoluttverdi 1,0029, dvs. at de ligger utenfor enhetssirkelen,
%og reguleringssystemet er derfor ustabilt (stemmer med simuleringen i 3b)!

%Plotter polene:
figure(6)
pzmap(M)
%Joda, to av polene synes å ligge utenfor enhetssirkelen, altså ustabilt system.

%---------------------------

%Deloppgave 5a:

%Velger Z-N-regulatorparameterverdier:
Kp=KpZN;
Ti=TiZN;
hrdisk=tf([Kp*(Ts+Ti),-Kp*Ti],[Ti,-Ti],Ts)
M=feedback(series(hrdisk,hpdisk),htilbakekopling)
%Bodeplott av M(z):
figure(7)
bode(M)
%Avlesning i Bodeplottet gir -3 dB-båndbredde wb lik ca 
wb=0.98 %rad/s.

%Reguleringssystemets statiske følgeegenskaper er gode siden
%amplitudekarakteristikken er på 0 dB = 1 for lave frekvenser.

%Anslag for responstiden:
Tr_bode=1.5/wb
%altså ca 1,5 sek.
%Avlesning av Tr fra sprangresponsen i figur 3 gir
Tr_sim=3
%Vel, anslaget er i hvertfall ikke helt vilt (feil med faktor på 2).

%---------------------------

%Deloppgave 5b:

%Som i 5a, men nå med økt Kp slik at systemet får dårlig stabilitet (men er stabilt):
Kp=KpZN*1.7;
Ti=TiZN;
hrdisk=tf([Kp*(Ts+Ti),-Kp*Ti],[Ti,-Ti],Ts)
M=feedback(series(hrdisk,hpdisk),htilbakekopling)

%Sjekker polene (for å bekrefte at reguleringssystemet fremdeles er stabilt):
[polerM,nullpunkterM]=pzmap(M)
abs(polerM)
%Ser at polene er innenfor enhetssirkelen (absoluttverdi mindre enn 1), altså stabilt system.

%Bodeplott av M(z):
figure(8)
bode(M)
%Vi ser at amplitudekurven får svært stor resonanstopp, typisk for systemer med svært dårlig stabilitet.