%   Script ZNTUNING.M
%
%   Reglereinstellung nach Ziegler-Nichols
%   Beispiel: Temperaturregelung eines Reaktors aus BSPREAK
%    Reaktordaten wie in RT1
%    hier um Messglied für Reaktortemperatur ergänzt
%
% J. Lunze
% 4.12.2002
% für 2. Auflage: 2.1.2007
% für 3. Auflage: 20.1.2012
% für 4. Auflage: 4.11.2015
% für 5. Auflage: 11.01.2020
%
%
clear
close all
%    Heizungsregelkreis für einen chem. Reaktor
%
%   Beispiel: Zylindrischer Behälter mit folgenden Parametern
%               Volumen     V = 10 000 cm^3
%               Höhe        h = 30 cm
%               Durchmesser d = 20 cm
%               Manteldicke     0.4 cm
%               Masse       m = 10 kg
%               spez. Wärme c = 4.19 kJ/kg K  (Wasser)
%               Durchfluss   w = 10 kg/min
%               Zeitkonstante der Temperaturmessung   Tm=0,2 min
%
%        Eingang: Heizleistung in kW
%        Ausgang: gemessene Temperatur der Flüssigkeit in oC
%        Störung: Temperatur des Zulaufes in oC
%        Zeit in min
%
%   [zs, ns] = Übertragungsfunktion Heizleistung -> Temperatur
%   [zd, nd] = Übertragungsfunktion Zulauftemperatur -> Temperatur
%
w=10;         % kg/min
m=10;         % kg         (Masser der Flüssigkeit im Reaktor)
c=4.19;       % kJ/(kg K)  (Wasser)
mh=7;         % kg         (Masse des Eisenmantels)
ch=0.47;      % kJ/(kg K)  (Stahl)
alpha=0.5;    % kW/(m^2 K)
A=0.22;       % m^2
Tm=0.2;       % min
%
Tend=7;
T=[0:0.01:Tend];
%
zs = [1/w/c*60];
ns = [mh*ch*m/w/alpha/A/60 (mh*ch/w/c + mh*ch/alpha/A/60 + m/w) 1];
Reaktorbehaelter=tf(zs,ns);
zm = [1];
nm = [Tm 1];
Messglied=tf(zm, nm);
Reaktor=Reaktorbehaelter*Messglied;
zd = [mh*ch/alpha/A/60 1];
nd = ns;
ReaktorbzglStoer=tf(zd, nd);
ReaktorStoer=ReaktorbzglStoer*Messglied;

%%   Übergangsfunktion des Reaktors
%
figure(1);
%
subplot(4,1,[1,2]);
[Y0, X]=step(Reaktor, T);
%   Wendetangente
t1=T(80);
h1=Y0(80);
t2=T(100);
h2=Y0(100);
Anstieg=(h2-h1)/(t2-t1);
t0=t1-h1/Anstieg;
ks=dcgain(Reaktor);
t3=t1+(ks-h1)/Anstieg;
plot(T, Y0, 'b', [t0, t3], [0, ks],'k--', [0, 7],[ks, ks],'k:', [2.54 2.54],[ks 0],'k:' );
latextitle('\"Ubergangsfunktion der Temperaturregelstrecke');
latexxlabel('\hspace{3cm} $$t$$ in min');
latexylabel('$$h$$ in K');
rotateY;
latextext(0.3, -0.25, '$$T_{\rm t}$$');
latextext(2.54, -0.25, '$$T+T_{\rm t}$$');
axis([0 Tend 0 1.6]);
xticks([1 2 4 5 6]);
yticks([0 0.5 1 1.5]);

%%    Reglereinstellung nach Ziegler-Nichols
%       für Tt= 0,4  [min]
%           T1 = 2,5  [min]
%
figure(2);
ks=1.4;
T1 = 2.5;
Tt= 0.4;
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/10 TD
zPID = kPID*[TPID*TD, (TPID+TD/10), (1+TPID*TD/10)];
nPID = [TPID*TD/10, 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);
%
subplot(4,1,[1,2]);
plot(T, YP,  'b',T, YPI,  'b',T, YPID, 'b', [0 7], [1 1], 'k:');
%latexxlabel('$$t$$ in min');
xleer;
xticksempty([0 2 4 6 8]);
latexylabel('$$h_{\rm w}$$ in K');
rotateY;
%grid('on');
latextitle('Temperaturregelkreis mit erster Reglereinstellung');
latextext(0.1, 1.4, 'PID');
latextext(1.6, 1.6, 'PI');
latextext(3.4, 0.7, 'P');
axis([0 Tend 0 1.8]);
yticks([0 0.5 1 1.5]);
%
%    Berechnung der Stellgröße u
StellgrP=feedback(PRegler, Reaktor);
StellgrPI=feedback(PIRegler,Reaktor);
StellgrPID=feedback(PIDRegler,Reaktor);
UP=step(StellgrP, T);
UPI=step(StellgrPI, T);
UPID=step(StellgrPID, T);
%
subplot(4,1,[3,4]);
plot(T, UP, 'b', T, UPID, 'b', [0 7], [1/ks 1/ks], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$u$$ in kW');
rotateY;
%grid('on');
latextext(0.1, -2, 'PID');
%latextext(2.5, 4, 'PI');
latextext(4, 2, 'P');
axis([0 Tend -3 6]);
yticks([-3 0 3 6]);
xticks([0 2 4 6 8]);



%% Störübergangsfunktion für den PID-Regler
%
figure(3);
subplot(4,1,[1,2]);
Kreis=feedback(1, offeneKettePID);
Stoerverhalten= series(ReaktorStoer, Kreis);
Ystor = step(Stoerverhalten, T);
plot(T, Ystor, 'b', [0 7],[0 0],'k:');
hold on
% Vergleich zum System ohne Regler
Yd=step(ReaktorStoer, T);
plot(T, Yd, 'b--');    
latexxlabel('$$t$$ in min');
latexylabel('$$h_{\rm d}$$ in K');
rotateY;
%grid('on');
axis([0, Tend, -0.2, 1.1]);
yticks([0 0.5 1]);
xticks([0 2 4 6 8]);
yticks([0 0.5 1 1.5]);
latextitle('St\"or\"ubergangsfunktion des PID-Regelkreises');
hold off


%%   Reglereinstellung nach dem zweiten Verfahren
%       Bestimmung der Schwingungsdauer
%
figure(4);
%
subplot(6,1,[3,4]);
kP= 9.12;
zP = [kP];                 % P-Regler [zP, nP]
nP = [1];
PRegler=tf(zP, nP);
%  geschlossene Kreise
offeneKetteP = series(PRegler, Reaktor);
RegelkreisP= feedback(offeneKetteP, 1);
YP = step(RegelkreisP, T);
plot(T, YP, 'b',  [0 7], [1 1], 'k:', [2.023 2.023], [1.3 0.6], 'k:', [3.53 3.53], [1.3 0.6], 'k:');
xleer;
xticksempty([0 2 4 6 8]);
latexylabel('$$h_{\rm w}$$ in K');
rotateY;
axis([0 Tend 0 2]);
yticks([0 1 2]);
latextext(2.1, 0.8, '$$\leftarrow$$');
latextext(3.1, 0.8, '$$\rightarrow$$');
latextext(2.8, 1.3, '$$T_{\rm krit}$$');
%
subplot(6,1,[1,2]);
kP= 10.5;
zP = [kP];                 % P-Regler [zP, nP]
PRegler=tf(zP, nP);
%  geschlossene Kreise
offeneKetteP = series(PRegler, Reaktor);
RegelkreisP= feedback(offeneKetteP, 1);
YP = step(RegelkreisP, T);
plot(T, YP, 'b',  [0 7], [1 1], 'k:');
latextitle('Schwingung des P-geregelten Reaktors');
xleer
xticksempty([0 2 4 6 8]);
latexylabel('$$h_{\rm w}$$ in K');
rotateY;
axis([0 Tend 0 2]);
yticks([0 1 2]);
%
subplot(6,1,[5,6]);
kP= 6.5;
zP = [kP];                 % P-Regler [zP, nP]
PRegler=tf(zP, nP);
%  geschlossene Kreise
offeneKetteP = series(PRegler, Reaktor);
RegelkreisP= feedback(offeneKetteP, 1);
YP = step(RegelkreisP, T);
plot(T, YP, 'b',  [0 7], [1 1], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$h_{\rm w}$$ in K');
rotateY;
axis([0 Tend 0 2]);
yticks([0 1 2]);
xticks([0 2 4 6 8]);


%%    Reglereinstellung nach Ziegler/Nichols
%      für  Tkrit = 1,5 min
%           kkrit = 9.1
%
%
figure(5);
kkrit=9.1;
Tkrit=1.5;
kP=0.5*kkrit;
kPI= 0.45*kkrit;
TPI=0.85*Tkrit;
kPID = 0.6*kkrit;
TPID = 0.5*Tkrit;
TD = 0.12*Tkrit;
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/10 TD
zPID = kPID*[TPID*TD, (TPID+TD/10), (1+TPID*TD/10)];
nPID = [TPID*TD/10, 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);
%
subplot(4,1,[1,2]);
plot(T, YP, 'b', T, YPI, 'b', T, YPID, 'b', [0 7], [1 1], 'k:');
%latexxlabel('$$t$$ in min');
xleer;
xticksempty([0 2 4 6 8]);
latexylabel('$$h_{\rm w}$$ in K');
rotateY;
%grid('on');
latextitle('Temperaturregelkreis mit zweiter Reglereinstellung');
latextext(0.1, 1.4, 'PID');
latextext(1.7, 1.6, 'PI');
latextext(3.4, 0.65, 'P');
axis([0 Tend 0 1.8]);
yticks([0 0.5 1 1.5]);
%    Berechnung der Stellgröße u
StellgrP=feedback(PRegler, Reaktor);
StellgrPI=feedback(PIRegler,Reaktor);
StellgrPID=feedback(PIDRegler,Reaktor);
UP=step(StellgrP, T);
UPI=step(StellgrPI, T);
UPID=step(StellgrPID, T);
%
subplot(4,1,[3,4]);
plot(T, UP, 'b', T, UPI, 'b', T, UPID, 'b', [0 7], [1 1], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$u$$ in kW');
rotateY;
%grid('on');
latextext(0.1, -2, 'PID');
latextext(2.5, 4, 'PI');
%latextext(3.6, 0.73, 'P');
axis([0 Tend -3 6]);
yticks([-3 0 3 6]);
xticks([0 2 4 6 8]);


%%     Bilder in .EPS-Dateien
%
Epsbildc15('zntuning');
