%   Script SMITHPR.M
%
%    Smithprädiktor für einen Kontireaktor
%
% J. Lunze
% 11.12.1995
% für 2. Auflage: 6.11.1998
% für 3. Auflage: 2.4.2001
% für 4. Auflage: 18.9.2003
% für 5. Auflage: 18.4.2005
% für 7. Auflage: 28.4.2008
% für 8. Auflage: 2.2.2010
% für 9. Auflage: 15.2.2012
% für 10. Auflage: 27.1.2014
% für 12. Auflage: 5.08.2019
%
%   
clear
close all
%%                  
%
%   Beispiel: Zylindrischer Behälter mit folgenden Parametern
%               Volumen     V = 10 dm^3
%               Durchfluß   w = 10 dm^3/min
%               Zeitkonstante der Konzentrationsmessung  Tm = 0.2 min
%               Totzeit zwischen Zulaufventil und Reaktor Tt = 0.4 min
%
%        Eingang: Volumenanteil eines Stoffes im Zulauf (Abweichung vom APkt)
%        Ausgang: Volumenanteil des Stoffes im Reaktor
%        Störung: Konzentrationsänderung im Reaktor
%        Zeit in min
%
%   Entwurf eines Smith-Prädiktors
%   
%
%  Daten aus BSPREA2.M
%
V = 10;     % dm^3
w = 10;     % dm^3/min
Tm = 0.2;   % min
Tt = 0.4;   % min
%Tt=1.5;
n = 2;
m = 1;
r = 1;
A = [-w/V 0
     1/Tm -1/Tm];
b = [w/V
     0];
c = [0 1];
d = 0;
Reaktor=ss(A, b, c, d);
%  Eigenwerte bei -1, -5

%%  Übergangsfunktion der Regelstrecke mit Totzeit 0,4
%   
figure(1)
subplot(2,1,1);
dT=0.1;
T = [0:dT:7];
[Y, X]= step(Reaktor, T);
kT=ceil(Tt/dT);
Ysize=size(Y);
for k=1:kT
   Y1(k)=0;
end;
for k=kT+1:Ysize(1)
   Y1(k)=Y(k-kT);
end;
plot(T, Y1, 'b', [0 7], [1 1], 'k:');
latextitle('Totzeitbehaftete Regelstrecke einschl. Messglied');
latexylabel('$$h$$ in $$\frac{\rm g}{\rm m^3}$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([0 7 -0.1 1.2]);
xticks([0 Tt 2 4 6 8 10]);
yticks([0 0.5 1]);


%%  PI-Regler an Regelstrecke ohne Totzeit
%   
figure(2);
subplot(2,1,1)
kP=2;
kI=2;      %  Regler als Übertragungsfunktion: K(s) = k_P(1+ 1/(T_Is)
           %   mit  k_P = 2;  T_I = 1
Aquer=[A-b*kP*c b*kI;0 -1 0];
bquer=[b*kP; 1];
cquer=[c 0];
Regelkreis=ss(Aquer, bquer, cquer, d);
[Yquer, Xquer]=step(Regelkreis, T);
plot(T, Yquer, 'b',  [-1 0 0 7], [0 0 1 1], 'b--');
axis([-1 7 -0.2 1.2]);
xticks([0 2 4 6 8]);
yticks([0 0.5 1]);
latextitle('Regelkreis mit PI-Regler und totzeitfreier Regelstrecke');
latexxlabel('$$t$$ in min');
latexylabel('$$y$$ in $$\frac{\rm g}{\rm m^3}$$');
rotateY;


%%   PI-Regler an Regelstrecke mit Totzeit
%   
figure(3);
subplot(6,1,[1,2])
Ts=0.02;
%  zeitdiskretes Streckenmodell einschließlich Totzeit
[Ad, bd, cd, dd]=c2dt(A, b, c, Ts, Tt);
StreckeDiskr=ss(Ad, bd, cd, dd);
%  zeitdiskreter PI-Regler [Ar, br, kI, kP]
[Ar, br]=c2d(0, 1, Ts);
Regler=ss(Ar, br, kI, kP);
%  geschlossener Kreis
offeneKetteDiskr=series(Regler, StreckeDiskr);
RegelkreisDiskr=feedback(offeneKetteDiskr, 1);
[Adk, bdk, cdk, ddk]=ssdata(RegelkreisDiskr);
[Ydk, Xdk]=dstep(Adk, bdk, cdk, ddk);
kend=length(Ydk);
Td=[0:kend-1];
Td=Td*Ts;
plot(Td, Ydk,'b--', [-1 0],[0 0],'b--');
axis([-1 7 -0.2 1.9]);
yticks([0 1]);
xleer;
latextitle('Regelkreis mit totzeitbehafteter Regelstrecke');
latexylabel('$$y$$ in $$\frac{\rm g}{\rm m^3}$$');
rotateY;
hold on;
latextext(2.5, 1.5, 'PI-Regler');
%
subplot(6,1,[3,4])
RegelkreisStellgr=feedback(Regler, StreckeDiskr);
[AdS, bdS, cdS, ddS]=ssdata(RegelkreisStellgr);
[Udk, Xdk]=dstep(AdS, bdS, cdS, ddS, 1, kend);
Udsize=size(Udk);
Td=[0:Udsize(1)-1];
Td=Td*Ts;
plot(Td, Udk,'b--', [-1 0 0], [0 0 Udk(1,1)],'b--');
hold on;
latexxlabel('$$t$$ in min');
latexylabel('$$u$$ in $$\frac{\rm g}{\rm m^3}$$');
rotateY;
axis([-1 7 -1 3]);
yticks([0 2]);
xleer;

%%   Realisierung des Reglers als Smith-Prädiktor
%   
figure(4);
subplot(2,1,1);
%  Kpi(s) * G(s) * (1-exp(-sTt))
%  zeitdiskretes Streckenmodell ohne Totzeit [Ado, bdo, cdo, ddo]
[Ado, bdo]=c2d(A, b, Ts);
cdo=c;
ddo=d;
StreckeOhneTotzeit=ss(Ado, bdo, cdo, ddo);
%  Berechnung der Stellgröße
RegelkreisOhneTotzeit=feedback(Regler, StreckeOhneTotzeit);
%  geschlossener Regelkreis: Stellgröße ist Eingang der totzeitbeh. Strecke
offeneKetteMitTotzeit=series(RegelkreisOhneTotzeit, StreckeDiskr);
RegelkreisMitTotzeit=feedback(offeneKetteMitTotzeit, 1);
[Aon, bon, con, don]=ssdata(offeneKetteMitTotzeit);
[Ysp, Xsp]=dstep(Aon, bon, con, don, 1, 7/Ts);
Yspsize=size(Ysp);
Time=[0:Yspsize(1)-1];
Time=Time*Ts;
plot(Time, Ysp, 'b', [-1 0],[0 0], 'b', [-1 0 0 7], [0 0 1 1], 'k:');
axis([-1 7 -0.1 1.2]);
yticks([0 0.5 1]);
latextitle('Regelkreis mit Smith-Pr\"adiktor');
latexylabel('$$y$$ in $$\frac{\rm g}{\rm m^3}$$');
rotateY;
xleer;
%
% Eintragen in figure(3)
figure(3)
subplot(6,1,[1,2])
plot(Time, Ysp, 'b', [-1 0],[0 0], [-1 0 0 7], [0 0 1 1], 'b--');
latextext(1.1, 0.3, 'Smith-Pr\"adiktor');
hold off;
%
%   Berechnung der Stellgröße des Smith-Prädiktors
figure(4)
subplot(2,1,2)
[AS, bS, cS, dS]=ssdata(RegelkreisOhneTotzeit);
Usp=dstep(AS, bS, cS, dS, 1, 7/Ts);
Uspsize=size(Usp);
Time=[0:Uspsize(1)-1];
Time=Time*Ts;
plot(Time, Usp, 'b', [-1 0 0],[0 0 Usp(1,1)], 'b', [-1 7], [0 0], 'k:');
latexylabel('$$u$$ in $$\frac{\rm g}{\rm m^3}$$');
latexxlabel('$$t$$ in min');
axis([-1 7 -0.2 2.3]);
yticks([0 1 2]);
xticks([0 2 4 6 8]);
rotateY;
%
%  dasselbe in figure(3) eintragen
figure(3)
subplot(6,1,[3,4]);
plot(Time, Usp, 'b', [-1 0 0],[0 0 Usp(1,1)], [-1 7], [0 0], 'k:');
hold off;
%  Vergleich des Reglereingangs bei PI-Regler und Smithprädiktor
subplot(6,1,[5,6])
%      Regelabweichung beim PI-Regler
plot(Td, 1-Ydk,'b--', [-1 0 0],[0 0 1-Ydk(1,1)],'b--', [-1 7], [0 0], 'k:');
hold on
%
%      Regelabweichung beim Smithprädiktor
offeneKetteRegelabw=series(Regler, StreckeOhneTotzeit);
Regelabweichung=feedback(1, offeneKetteRegelabw, -1);
[Ae, be, ce, de]=ssdata(Regelabweichung);
E=dstep(Ae, be, ce, de, 1, 7/Ts);
Esize=size(E);
Te=[0:Esize(1)-1];
Te=Te*Ts;
plot(Te, E, 'b', [-1 0 0],[0 0 E(1,1)], 'b');
hold on
latexylabel('$$e$$ in $$\frac{\rm g}{\rm m^3}$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([-1 7 -0.7 1.2]);
xticks([0 2 4 6 8]);
yticks([0 1]);
hold off;


%%  Smith-Prädiktor für unterschiedliche Totzeiten
%   mit denselben PI-Reglerparametern
%   aber korrektem Streckenmodell
%   
%  Kpi(s) * G(s) * (1-exp(-sTt))
%  zeitdiskretes Streckenmodell ohne Totzeit
%
figure(5);
subplot(2,1,1)
plot(Time, Ysp, 'b');
latextitle('Regelkreis mit Smith-Pr\"adiktor bei $$T_t=0 ... 1.2$$');
latexylabel('$$y$$ in $$\frac{\rm g}{\rm m^3}$$');
rotateY;
latexxlabel('$$t$$ in min');
hold on;
for Tt1=0.01:0.3:1.2                  % MATLAB-Fehlermeldung bei Benutzung von c2dt mit Tt1=0
   [Ad1, bd1, cd1, dd1]=c2dt(A, b, c, Ts, Tt1);
   Strecke1=ss(Ad1, bd1, cd1, dd1);
   offeneKetteMitTotzeit2=series(RegelkreisOhneTotzeit, Strecke1);
   [Aon1, bon1, con1, don1]=ssdata(offeneKetteMitTotzeit2);
  Ysp1=dstep(Aon1, bon1, con1, don1, 1, length(Time));
  plot(Time, Ysp1, 'b');
  axis([0 7 0 1.2]);
end;
axis([0 7 0 1.3]);
yticks([0 0.5 1]);
xticks([0 2 4 6]);
hold off

%%  Smith-Prädiktor für unterschiedliche Totzeiten
%   mit denselben PI-Reglerparametern
%   aber fehlerbehaftetem Streckenmodell  (Totzeit falsch)
%   
%  Kpi(s) * G(s) * (1-exp(-sTt))
%  zeitdiskretes Streckenmodell ohne Totzeit
%
figure(6);
plot(Time, Ysp,'b--');
latextitle('Regelkreis mit Smith-Pr\"adiktor f\"ur fehlerbehaftete Totzeit $$T_t=0 ... 1.2$$');
latexylabel('$$y$$ in $$\frac{\rm g}{\rm m^3}$$');
rotateY;
latexxlabel('$$t$$ in min');
hold on;
for Tt2=[0.01, 0.3, 0.6, 0.9, 1.2]
   [Ad2, bd2, cd2, dd2]=c2dt(A, b, c, Ts, Tt2);
   StreckeMitTotzeit2=ss(Ad2, bd2, -cd2, -dd2);
  Parallelschaltung=parallel(StreckeDiskr, StreckeMitTotzeit2);
  Regelkreis3=feedback(RegelkreisOhneTotzeit, Parallelschaltung);
  RegelkreisMitTotzeit3=series(Regelkreis3, StreckeDiskr);
  [Aksp2, bksp2, cksp2, dksp2]=ssdata(RegelkreisMitTotzeit3);
  Ysp2=dstep(Aksp2, bksp2, cksp2, dksp2, 1, length(Time));
  plot(Time, Ysp2, 'b');
end;
latextext(0.3, 1.9, 'Totzeit');
latextext(0.3, 1.75, 'der Strecke');
latextext(0.3, 1.4, '0');
latextext(0.3, 1.15, '0,3');
latextext(0.3, 1.0, '0,4');
latextext(0.3, 0.85, '0,6');
latextext(0.3, 0.6, '0,9');
latextext(1.3, 0.4, '1,2');
axis([0 7 -0.1 2]);
yticks([0 0.5 1 1.5 2]);
xticks([0 2 4 6]);
hold off;


%%   Abschätzung des möglichen Totzeitfehlers
%   (Sicherung der robusten Stabilität)
%
%   
figure(7);
subplot(2,1,1);
%  Amplitudengang von
%  |(PI-Regler*G ohne Totzeit)(1+PI-Regler*G ohne Totzeit) * exp(-0.4s)|
%  Totzeitglied exp(-0.4s) kann weggelassen werden, weil sie den Betrag
%  nicht verändert
PIRegler=ss(0, 1, kI, kP);
offeneKette=series(PIRegler, Reaktor);
Regelkreis=feedback(offeneKette, 1);
O=logspace(-1, 1, 300);
[Betrag, Phase]=bode(Regelkreis, O);
semilogx(O, 20*log10(Betrag(:)), 'b');
hold on;
%  Amplitudengang von 1/|1-exp(-dTt s)|
Osize=size(O);
for dTt=[-0.55:0.05:-0.05, 0.05:0.05:0.55]
    for k=1:Osize(2)
    Rob(k, 1)=1/abs(1-exp(-i*dTt*O(k)));
    end;
    semilogx(O, 20*log10(Rob(:)), 'b--');
end;
latextitle('Amplitudengang zur Robustheitsanalyse');
latexylabel('$$|G|$$ in dB');
rotateY;
latexxlabel('$$\omega$$ in rad/s');
xticks([0.1 1 10]);
axis([min(O), max(O), -20, 20]);
xticks([0.1 1 10]);
yticks([-20 0 20]);
hold off;

%%     Bilder auf .EPS- und .PS-Dateien
%
Epsbildc15('smithpr');
