%   Script PITUNING.M
%
%   I-Reglereinstellung
%   Beispiel: Temperaturregelung eines Reaktors wie in ZIEGLER.M
%
% J. Lunze
% 19.5.2002
% für 2. Auflage: 23.6.2007
% für 3. Auflage: 20.1.2012
% für 4. Auflage: 5.1.2016
% für 5. Auflage: 12.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)
%               Durchfluß   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

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]);
T=[0:0.01:7];
[Y, X]=step(Reaktor, T);
ks=dcgain(Reaktor)
plot(T, Y, '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]);
yticks([0 0.5 1 1.5]);
xticks([0 2 4 6]);

%%    Einstellung des I-Reglers
%
figure(2);
subplot(3,1,[1,2]);
T=[0:0.01:12];
nI = [1, 0];
for kI= [0.05 0.1 0.3 0.5];
    zI = kI;        % I-Regler [kI]
    IRegler=tf(zI, nI);
    %  geschlossene Kreise
    offeneKetteI = series(IRegler, Reaktor);
    RegelkreisI= feedback(offeneKetteI, 1);
    YI = step(RegelkreisI, T);
    plot(T, YI, 'b');
    hold on;
end
plot([0 20],[1 1],'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$h_{\rm w}$$');
rotateY;
%grid('on');
latextitle('F\"uhrungs\"ubergangsfunktion mit I-Temperaturregler');
axis([0 12 0 1.8]);
yticks([0 0.5 1 1.5]);
xticks([0 2 4 6 8 10 12 14]);
latextext(6, 0.2, '0,05');
latextext(6, 0.45, '0,1');
latextext(8, 1.2, '0,3');
latextext(4, 1.4, '$$k_{\rm I}=0,\!5$$');
hold off

%%  Lage der Pole bei unterschiedlicher Reglerverstärkung
%
figure(3)
Pole=[];
k1=0;
for kI=[0.001 0.1:0.0005:0.2 0.21:0.01:0.5]
    zI = kI;        % I-Regler [kI]
    IRegler=tf(zI, nI);
    %  geschlossene Kreise
    offeneKetteI = series(IRegler, Reaktor);
    RegelkreisI= feedback(offeneKetteI, 1);
    [num, den]=tfdata(RegelkreisI,'v');
    k1=k1+1;
    Pole(k1,:)=roots(den);
end
[Anz, Ord]=size(Pole);
for k1=1:Ord
    plot(real(Pole(:,k1)), imag(Pole(:,k1)),  'b')
    hold on
end
[num, den]=tfdata(Reaktor,'v');
Pole=roots(den);
plot(real(Pole), imag(Pole), 'bx');
plot(0, 0, 'bx');
kI=0.05;
IRegler=tf(kI, nI);
offeneKetteI = series(IRegler, Reaktor);
RegelkreisI= feedback(offeneKetteI, 1);
[num, den]=tfdata(RegelkreisI,'v');
Pole=roots(den);
plot(real(Pole), imag(Pole), 'bd');
kI=0.1;
IRegler=tf(kI, nI);
offeneKetteI = series(IRegler, Reaktor);
RegelkreisI= feedback(offeneKetteI, 1);
[num, den]=tfdata(RegelkreisI,'v');
Pole=roots(den);
plot(real(Pole), imag(Pole), 'bs');
kI=0.3;
IRegler=tf(kI, nI);
offeneKetteI = series(IRegler, Reaktor);
RegelkreisI= feedback(offeneKetteI, 1);
[num, den]=tfdata(RegelkreisI,'v');
Pole=roots(den);
plot(real(Pole), imag(Pole), 'bo');
kI=0.5;
IRegler=tf(kI, nI);
offeneKetteI = series(IRegler, Reaktor);
RegelkreisI= feedback(offeneKetteI, 1);
[num, den]=tfdata(RegelkreisI,'v');
Pole=roots(den);
plot(real(Pole), imag(Pole), 'b*');
plot([0 0],[-1 1],'k:', [-2 1],[0 0],'k:');
latexxlabel('Re');
latexylabel('Im');
rotateY;
axis([-1.2 0.4 -0.8 0.8]);
axis('square');
yticks([-0.5 0 0.5]);
xticks([-1 -0.5 0 0.5]);
latextitle('Dominante Pole bei der Reglereinstellung');
hold off

%%  Hinzunahme des P-Anteils
%
figure(4);
subplot(3,1,[1,2]);
kI=0.5;
nPI = [1 0];
for kP= [0  0.3  0.6 1 2 4 ];
    zPI = [kP kI];        % PI-Regler
    PIRegler=tf(zPI, nPI);
    %  geschlossene Kreise
    offeneKettePI = series(PIRegler, Reaktor);
    RegelkreisPI= feedback(offeneKettePI, 1);
    YPI = step(RegelkreisPI, T);
    plot(T, YPI, 'b');
    hold on;
end
plot([0 20],[1 1],'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$h_{\rm w}$$');  % in K
rotateY;
%grid('on');
latextitle('F\"uhrungs\"ubergangsfunktion mit PI-Temperaturregler');
axis([0 12 0 1.8]);
yticks([0 0.5 1 1.5]);
xticks([0 2 4 6 8 10 12 14]);
latextext(8.5, 1.7, '$$k_{\rm I}=0,\!5$$');
latextext(4.5, 1.25, '0,3');
latextext(4.5, 1.07, '0,6');
latextext(0.9, 0.6, '1');
latextext(0.8, 0.9, '2');
latextext(4, 1.5, '$$k_{\rm P}=0$$');
latextext(0.4, 1.2, '4');
hold off

%%  gemeinsame Einstellung von P- und I-Anteil
%
figure(5);
subplot(3,1,[1,2]);
kI=0.5;
kP=2;
zPI0=[kP, kI];
nPI = [1, 0];
for a=[0.4, 0.6, 1, 2]
    zPI=a*zPI0;
    PIRegler=tf(zPI, nPI);
    %  geschlossene Kreise
    offeneKettePI = series(PIRegler, Reaktor);
    RegelkreisPI= feedback(offeneKettePI, 1);
    YPI = step(RegelkreisPI, T);
    plot(T, YPI, 'b');
    hold on;
end
plot([0 20],[1 1],'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$h_{\rm w}$$');  % in K
rotateY;
%grid('on');
latextitle('F\"uhrungs\"ubergangsfunktion mit PI-Temperaturregler');
axis([0 12 0 1.8]);
yticks([0 0.5 1 1.5]);
xticks([0 2 4 6 8 10 12 14]);
latextext(8.5, 1.7, '$$k_{\rm I}=0,\!5$$');
latextext(8.5, 1.5, '$$k_{\rm P}=2$$');
latextext(0.5, 1.6, '$$a=$$');
latextext(0.5, 1.35, '$$2$$');
latextext(1.1, 1.2, '1');
%latextext(1.1, 0.8, '0.6');
latextext(3.5, 0.65, '0.4');
hold off

%%   Stellgröße bei der Einstellung des I-Reglers
%
figure(6);
subplot(3,1,[1,2]);
ks=dcgain(Reaktor);
T=[0:0.01:20];
nI = [1, 0];
for kI= [0.05 0.1 0.3 0.5];
    zI = kI;        % I-Regler [kI]
    IRegler=tf(zI, nI);
    %  geschlossene Kreise
    StellgrRegelkreisI= feedback(IRegler, Reaktor);
    YI = step(StellgrRegelkreisI, T);
    plot(T, YI, 'b');
    hold on;
end
plot([0 20],[1/ks 1/ks],'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$u$$');
rotateY;
latextitle('Stellgroesse bei I-Temperaturregler');
axis([0 20 0 1.4]);
yticks([0 0.5 1 1.5]);
xticks([0 5 10 15 20]);
latextext(6, 0.1, '0.05');
latextext(6, 0.4, '0.1');
latextext(7, 0.9, '0.3');
latextext(0.4, 1.2, '$$k_{\rm I}=0,\!5$$');
latextext(20.5, 1/ks, '$$\frac{1}{k_{\rm s}}$$');
hold off

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