%   Script ZIEGLER.M
%
%   Reglereinstellung nach Ziegler-Nichols
%   Beispiel: Temperaturregelung eines Reaktors
%
% J. Lunze
% 13.11.1995
% für 2. Auflage: 5.11.1998
% für 3. Auflage: 3.4.2001
% für 4. Auflage: 15.9.2003
% für 5. Auflage: 18.4.2005
% für 6. Auflage: 4.12.2006
% für 7. Auflage: 28.4.2008
% für 9. Auflage: 10.2.2012
% für 10. Auflage: 27.1.2014
% für 12. Auflage: 11.10.2019
%
%   
clear
close all
Bspreak;     %  "Reaktor" = Modell des Temperaturverhaltens

%%   Übergangsfunktion des Reaktors
%   
figure(1);
subplot(2, 1, 1)
T=[0:0.01:7];
[Y, X]=step(Reaktor, T);
%   Wendetangente
t1=T(80);
h1=Y(80);
t2=T(100);
h2=Y(100);
Anstieg=(h2-h1)/(t2-t1);
t0=t1-h1/Anstieg;
ks=dcgain(Reaktor);
t3=t1+(ks-h1)/Anstieg;
plot(T, Y, 'b', [t0, t3], [0, ks],'b--', [0, 7],[ks, ks],'k:' );
latextitle('\"Ubergangsfunktion der Temperaturregelstrecke');
latexxlabel('$$t$$ in min');
latexylabel('$$h$$ in K');
rotateY;
axis([0 6 0 1.6]);
xticks([0 2 4 6 8]);
yticks([0 1]);



%%  Reglereinstellung nach Ziegler-Nichols
%       für Tt= 0,18  [min]
%           T1 = 2,1  [min]
%   
figure(2);
subplot(3,1,[1,2])
T1 = 2.1;
Tt= 0.18;
kP=1/ks*T1/Tt;
kPI= 0.9*kP;
TPI=3.33*Tt;
kPID = 1.2*kP;
TPID = 2*Tt;
TD = 0.5*Tt;
zP = [kP];                 % P-Regler [zP, nP]
nP = [1];
PRegler=tf(zP, nP);
zPI = kPI*[TPI, 1];        % PI-Regler [zPI, nPI]
nPI = [TPI, 0];
PIRegler=tf(zPI, nPI);
%  PID-Regler [zPID, nPID] mit verzögertem D-Anteil
%  Verzögerungszeit 1/20 TD
zPID = kPID*[TPID*TD, TPID, 1];
nPID = [TPID*TD/20, TPID, 0];
PIDRegler=tf(zPID, nPID);
%  geschlossene Kreise
offeneKetteP = series(PRegler, Reaktor);
RegelkreisP= feedback(offeneKetteP, 1);
offeneKettePI = series(PIRegler, Reaktor);
RegelkreisPI = feedback(offeneKettePI, 1);
offeneKettePID= series(PIDRegler, Reaktor);
RegelkreisPID = feedback(offeneKettePID,1);
YP = step(RegelkreisP, T);
YPI = step(RegelkreisPI, T);
YPID = step(RegelkreisPID, T);
plot(T, YP, 'b', T, YPI, 'b', T, YPID, 'b', [0 5], [1 1], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$h_{\rm w}$$');  % in K
rotateY;
%grid('on');
latextitle('F\"uhrungs\"ubergangsfunktion des Temperaturregelkreises');
latextext(0.1, 1.4, 'PID');
latextext(0.9, 1.5, 'PI');
latextext(2.1, 0.85, 'P');
axis([0 4 0 1.8]);
xticks([0 1 2 3 4]);
yticks([0 0.5 1 1.5]);


%%   Vergrößerung der Integratorzeitkonstante
%    auf das Doppelte
%   
figure(3);
subplot(3,1,[1,2])
kPID = 1.2*kP;
TPID = 4*Tt;
TD = 0.5*Tt;
%  PID-Regler [zPID, nPID] mit verzögertem D-Anteil
%  Verzögerungszeit 1/20 TD
zPID = kPID*[TPID*TD, TPID, 1];
nPID = [TPID*TD/20, TPID, 0];
PIDReglerneu=tf(zPID, nPID);
offeneKettePIDneu = series(PIDReglerneu, Reaktor);
RegelkreisPIDneu= feedback(offeneKettePIDneu, 1);
YPIDneu = step(RegelkreisPIDneu, T);
plot(T, YPIDneu, 'b', T, YPID, 'b', [0 5], [1 1], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$h_{\rm w}$$');
rotateY;
%grid('on');
axis([0, 4, 0, 1.8]);
xticks([0 1 2 3 4]);
yticks([0 0.5 1 1.5]);
latextitle('F\"uhrungs\"ubergangsfunktion des PID-Regelkreises');
latextext(0.8, 1.45, 'alter Regler');
latextext(1.05, 1.25, 'neuer Regler');


%%  Störübergangsfunktion für den PID-Regler
%   
figure(4);
subplot(2,1,1)
Kreis=feedback(1, offeneKettePIDneu);
Stoerverhalten= series(ReaktorStoer, Kreis);
Ystor = step(Stoerverhalten, T);
plot(T, Ystor, 'b', [0 4],[0 0],'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$h_{\rm d}$$ in K');
rotateY;
%grid('on');
axis([0, 4, -0.05, 0.15]);
xticks([0 1 2 3 4]);
yticks([0 0.1]);
latextitle('St\"or\"ubergangsfunktion des PID-Regelkreises');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%    Vergleich von Einstellregeln mit WOK-Entwurf
%
%  Gütevorgaben: Überschwingen deltah < 5%
%                    deshalb Knickpunktabstand a > 2
%                Einschwingzeit T_m < 3 min
%                    deshalb omegas > 1 rad/min
%
%
%
figure(5)
a=2;
omegas=1;
%   Bodediagramm der Regelstrecke
W=logspace(-1, 2, 1000);
[BReaktor, PReaktor]=bode(Reaktor, W);
subplot(2, 1, 1);
semilogx(W, 20*log10(BReaktor(:)), 'b');
hold on;
%    Geradenapproximation der gewünschten offenen Kette
plot(omegas, 0, 'bx');
latextext(0.95*omegas, -20, '$$\omega_{\rm s}$$');
plot(a*omegas, 0, 'bx');
omega0=0.1;
omega2=100;
omega1=a*omegas;
G0=20*log10(omegas/omega0);    % Startpunkt der Geradenapproximation
G1=-20*log10(omega1/omegas);
G2=G1-40*log10(omega2/omega1);
plot([omega0 omegas omega1 omega2],[G0 0 G1 G2],'b--');
latextext(0.95*omegas*a, -20, '$$\omega_1$$');
axis([0.1 100 -60 20]);
yticks([-60 -40 -20 0 20]);
latextitle('Bodediagramm des Reaktors');
xleer;
grid('on');
latexylabel('$$|G|$$ in dB');
rotateY;
hold off;
%
subplot(2, 1, 2);
semilogx(W, PReaktor(:), 'b');
hold on
axis([0.1 100 -200 0]);
latexxlabel('$$\omega$$ in rad/min');
xticks([0.1 1 10 100]);
latexylabel('$$\phi$$ in Grad');
rotateY;
yticks([-270 -180 -90 0]);
grid('on');

%%   PI-Regler
%   Knickfrequenz bei omegas
%   Verstärkung 1
figure(6)
kP=1;
TI=omegas;
PIRegler=tf(kP*[TI 1],[TI 0]);
offeneKette=PIRegler*Reaktor;
%   Amplitudengang der offenen Kette
[B, P]=bode(offeneKette, W);
subplot(2, 1, 1);
semilogx(W, 20*log10(B(:)), 'b');
hold on;
%    Geradenapproximation der gewünschten offenen Kette
plot(omegas, 0, 'bx');
latextext(0.95*omegas, -20, '$$\omega_{\rm s}$$');
plot(a*omegas, 0, 'bx');
plot([omega0 omegas omega1 omega2],[G0 0 G1 G2],'b--');
latextext(0.95*omegas*a, -20, '$$\omega_1$$');
axis([0.1 100 -60 20]);
yticks([-60 -40 -20 0 20]);
latextitle('Bodediagramm der offenen Kette');
latexxlabel('$$\omega$$ in rad/min');
xticks([0.1 1 10 100]);
grid('on');
latexylabel('$$|G|$$ in dB');
rotateY;
hold off;


%%   Störverhalten des Regelkreises
%
figure(7)
subplot(3,1,[1,2])
RegelkreisEntw = feedback(offeneKette,1);
YEntw = step(RegelkreisEntw, T);
plot(T, YEntw, 'b', T, YPI, 'b--', [0 5], [1 1], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$h_{\rm w}$$ in K');
rotateY;
%grid('on');
latextitle('F\"uhrungs\"ubergangsfunktion des Temperaturregelkreises');
axis([0 4 0 1.8]);     % wie bei Tuning-Regelung
xticks([0 1 2 3 4]);
yticks([0 0.5 1 1.5]);


%%   Schnelleres Einschwingen
%
%
%  Gütevorgaben: Überschwingen deltah < 5%
%                    deshalb Knickpunktabstand a > 2 (wie bisher)
%                Einschwingzeit T_m < 1 min
%                    deshalb omegas > 3 rad/min
%
%
%
figure(8)
a=2;
omegas=3;
%   Bodediagramm der Regelstrecke
subplot(2, 1, 1);
semilogx(W, 20*log10(BReaktor(:)), 'b');
hold on;
%    Geradenapproximation der gewünschten offenen Kette
plot(omegas, 0, 'bx');
latextext(0.95*omegas, 20, '$$\omega_{\rm s}$$');
plot(a*omegas, 0, 'bx');
omega0=0.1;
omega2=100;
omega1=a*omegas;
G0=20*log10(omegas/omega0);    % Startpunkt der Geradenapproximation
G1=-20*log10(omega1/omegas);
G2=G1-40*log10(omega2/omega1);
plot([omega0 omegas omega1 omega2],[G0 0 G1 G2],'b--');
latextext(0.95*omegas*a, 20, '$$\omega_1$$');
axis([0.1 100 -60 40]);
yticks([-60 -40 -20 0 20 40]);
latextitle('Bodediagramm des Reaktors');
latexxlabel('$$\omega$$ in rad/min');
xticks([0.1 1 10 100]);
grid('on');
latexylabel('$$|G|$$ in dB');
rotateY;
hold off;


%%   PI-Regler wie bisher, aber mit größerer Verstärkung
%   Knickfrequenz wie bisher
%   Verstärkung 3
figure(9)
kP=3;
PIReglerNeu=tf(kP*[TI 1],[TI 0]);
offeneKetteNeu=PIReglerNeu*Reaktor;
%   Amplitudengang der offenen Kette
[B, P]=bode(offeneKetteNeu, W);
subplot(2, 1, 1);
semilogx(W, 20*log10(B(:)), 'b');
hold on;
%    Geradenapproximation der gewünschten offenen Kette
plot(omegas, 0, 'bx');
latextext(0.95*omegas, 20, '$$\omega_{\rm s}$$');
plot(a*omegas, 0, 'bx');
plot([omega0 omegas omega1 omega2],[G0 0 G1 G2],'b--');
latextext(0.95*omegas*a, 20, '$$\omega_1$$');
axis([0.1 100 -60 40]);
yticks([-60 -40 -20 0 20 40 60]);
latextitle('Bodediagramm der offenen Kette');
latexxlabel('$$\omega$$ in rad/min');
xticks([0.1 1 10 100]);
grid('on');
latexylabel('$$|G|$$ in dB');
rotateY;
hold off;

%%   zusätzliches Korrekturglied
%   phasenanhebend zwischen omegas und omega1
figure(10)
TD=1/omegas;
T1=1/omega1;
Korrglied=tf([TD 1], [T1 1]);
offeneKetteKorr=offeneKetteNeu*Korrglied;
%   Amplitudengang der offenen Kette
[B, P]=bode(offeneKetteKorr, W);
subplot(2, 1, 1);
semilogx(W, 20*log10(B(:)), 'b');
hold on;
%    Geradenapproximation der gewünschten offenen Kette
plot(omegas, 0, 'bx');
latextext(0.95*omegas, 20, '$$\omega_{\rm s}$$');
plot(a*omegas, 0, 'bx');
plot([omega0 omegas omega1 omega2],[G0 0 G1 G2],'b--');
latextext(0.95*omegas*a, 20, '$$\omega_1$$');
axis([0.1 100 -60 40]);
yticks([-60 -40 -20 0 20 40 60]);
latextitle('Bodediagramm der offenen Kette');
latexxlabel('$$\omega$$ in rad/min');
xticks([0.1 1 10 100]);
grid('on');
latexylabel('$$|G|$$ in dB');
rotateY;
hold off;



%%   Störverhalten des Regelkreises mit neuem Regler
%
figure(11)
subplot(3,1,[1,2])
RegelkreisKorr = feedback(offeneKetteKorr,1);
YKorr = step(RegelkreisKorr, T);
plot(T, YKorr, 'b', T, YPI, 'b--', [0 5], [1 1], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$h_{\rm w}$$ in K');
rotateY;
%grid('on');
latextitle('F\"uhrungs\"ubergangsfunktion des Temperaturregelkreises');
axis([0 4 0 1.8]);     % wie bei Tuning-Regelung
xticks([0 1 2 3 4]);
yticks([0 0.5 1 1.5]);


%%     Bilder auf .EPS-Dateien
%
Epsbildc15('ziegler');
