%   Script TUNAHL.M zur Herstellung von Bildern
%
%   Einstellung von PI-Reglern am Beispiel der AHL-Anlage
%
% verwendet EPSBILD.M
% Daten aus BSPAHL.M
%
% J. Lunze
% 12.1.1997             Daten wie 21.3.1996
% für 2. Auflage: 16.12.2001
% für 3. Auflage: 19.4.2004
% für 5. Auflage: 18.10.2007
% für 7. Auflage 31.3.2012
% für 10. Auflage: 11.9.2019
%
%
clear
close all
%           Daten aus BSPAHL.M
%   Daten aus (Lunze, v.Kurnatowski msr 1987 (zeitdiskret)
%               --> Lunze 1988 (kontinuierlich) )
%
%        Eingang: u1 = Zuflussmenge der Harnstofflösung
%                 u2 = Temperatur des Eingangsgemisches
%                      (Sollwert für unterlagerten Regelkreis)
%        Ausgang: y1 = Harnstoffkonzentration
%                 y2 = Stickstoffkonzentration
%
%        Zeit in Minuten
%
%   [A B C D]  Zustandsraummodell
%
%
%             Modell für Zeiteinheit Sekunden
%             Daten aus o.g. Literaturstelle; System ist nicht vollständig
%             steuerbar
%
A = [-0.00606    0      0        0       0        0        0        0
      0.00166 -0.00166  0        0       0        0        0        0
          0      0   -0.00641    0       0        0        0        0
          0      0    0.00084 -0.00084   0        0        0        0
          0      0      0        0    -0.00641    0        0        0
          0      0      0        0     0.00166 -0.00166    0        0
          0      0      0        0       0        0     -0.00084    0
          0      0      0        0       0        0      0.0004  -0.0004 ];
B = [0.00606     0;
       0         0;
       0       0.00641;
       0         0;
     0.00641     0;
       0         0;
       0        0.00084;
       0         0];
C = [ 0  2.2  0  0.033  0      0     0      0
      0  0    0    0    0   -0.26    0    0.09 ];
D = [0 0
     0 0];
n=8;
m=2;
r=2;
%       Umrechnung auf Zeiteinheit Minuten
%
k=60;
As=k*A;
Bs=k*B;
Cs=C;
Ds=D;
AHLAnlage=ss(As, Bs, Cs, Ds);
%
%     2 zusätzliche Messglieder, Zeitkonstante 2 Minuten
%
Am=[-0.5 0
     0   -0.5];
Bm=[0.5  0
     0  0.5];
Cm=[1  0
    0  1];
Dm=[0  0
    0  0];
Messglieder=ss(Am, Bm, Cm, Dm);
AHLAnlageM=series(AHLAnlage, Messglieder);
n=n+2;

%%   Statischer Verstärkungsfaktor
%
%
Ks=dcgain(AHLAnlageM)
Ksinv=inv(Ks);

%%  Übergangsfunktionsmatrix der Regelstrecke
%
figure(1);
Te=250;
T=[0:Te];
Y=step(AHLAnlageM, T);
subplot(2, 2, 1);
plot(T, Y(:,1,1),'b', [0 Te], [Ks(1,1), Ks(1,1)], 'k:');
latexylabel('$$h_{11}$$');
rotateY;
xleer;
latextitle('\"Ubergangsfunktionsmatrix');
axis([0 Te 0 2.5]);
yticks([0 2]);
subplot(2, 2, 2);
plot(T, Y(:,1,2),'b', [0 Te], [Ks(1,2), Ks(1,2)], 'k:');
yticks([0 0.04]);
latexylabel('$$h_{12}$$');
rotateY;
xleer;
latextitle(' AHL-Anlage');
axis([0 Te 0 0.04]);
subplot(2, 2, 3);
plot(T, Y(:,2,1),'b', [0 Te], [Ks(2,1), Ks(2,1)], 'k:');
yticks([-0.2 0]);
latexylabel('$$h_{21}$$');
rotateY;
latexxlabel('$$t$$ in min');
latextitle('');
axis([0 Te -0.3 0]);
xticks([0 100 200]);
%
subplot(2, 2, 4);
plot(T, Y(:,2,2),'b', [0 Te], [Ks(2,2), Ks(2,2)], 'k:');
latexylabel('$$h_{22}$$');
rotateY;
yticks([0 0.1]);
latexxlabel('$$t$$ in min');
latextitle('');
axis([0 Te 0 0.1]);
xticks([0 100 200]);



%%
%
%   Entkopplungseigenschaft:
%     Übergangsfunktion für Regelstrecke mit Eingang Ks^-1 * u
%
%
figure(2);
Entkopplungsglied=ss(zeros(m,m), zeros(m,m), zeros(m,m), Ksinv);
ModifizAHLAnlageM=series(Entkopplungsglied, AHLAnlageM);
Y=step(ModifizAHLAnlageM, T);
subplot(2, 2, 1);
plot(T, Y(:,1,1), 'b', [0 Te], [0 0], 'k:');
latexylabel('$$y_1$$');
rotateY;
yticks([0 1]);
xleer;
axis([0 Te -0.2 1.2]);
%latextitle('\"Ubergangsfunktionsmatrix');
subplot(2, 2, 2);
plot(T, Y(:,1,2),'b', [0 Te], [0 0], 'k:');
yleer;
yticks([0 1]);
%latexylabel('$$y_1$$');
rotateY;
axis([0 Te -0.2 1.2]);
xleer;
latextitle('AHL-Anlage mit Eingang $$K_s^{-1} *u$$');
%
subplot(2, 2, 3);
plot(T, Y(:,2,1),'b', [0 Te], [0 0], 'k:');
latexylabel('$$y_2$$');
rotateY;
yticks([0 1]);
latexxlabel('$$t$$ in min');
latextitle('');
axis([0 Te -0.2 1.2]);
xticks([0 100 200]);
%
subplot(2, 2, 4);
plot(T, Y(:,2,2),'b', [0 Te], [0 0], 'k:');
yleer;
yticks([0 1]);
%latexylabel('$$y_2$$');
%rotateY;
latexxlabel('$$t$$ in min');
latextitle('');
axis([0 Te -0.2 1.2]);
xticks([0 100 200]);


%%     Wurzelortskurve
%
%
figure(3);
[A, B, C, D]=ssdata(AHLAnlageM);
Wurzelortskurve=ss([A zeros(n,m); C zeros(m,m)], [B; D], [C zeros(m,m)], zeros(m, m));
%  Pole und Nullstellen der offenen Kette
[Pole, Nullstellen]=pzmap(Wurzelortskurve);
for k=1:length(Pole)
    plot(real(Pole(k)), imag(Pole(k)), 'bx');
    hold on;
end
%for k=1:length(Nullstellen)
%    plot(real(Nullstellen(k)),imag(Nullstellen(k)),'bo');
%end
k=1;
R=[];
for a=[0.00001:0.0001:0.06, 0.06:0.001:0.15]
    KI=a*Ksinv;
    Ag=[A -B*KI;
        C zeros(m,m)];
    d=roots(poly(Ag));
    R(k,:)=d.';
    k=k+1;
end;
plot(R, 'b.');
plot([-1 1],[0 0],'k:', [0 0],[-1 1],'k:');
axis([-0.1 0.1 -0.1 0.1]);
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Wurzelortskurve der AHL-Anlage bez\"uglich des I-Tuningfaktors');
hold off
xticks([-0.1 0 0.1]);
yticks([-0.1 0 0.1]);



%%  Detailuntersuchung der dominierenden Pole
%  Markieren der Pole für die später verwendeten Einstellungen
%
figure(4)
plot(R,'b.');
hold on
plot([-1 1],[0 0],'k:', [0 0],[-1 1],'k:');
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Wurzelorte f\"ur $$a=0.005, 0.01, 0.02$$');
a=0.005;
KI=a*Ksinv;
Ag=[A -B*KI; C zeros(m,m)];
Pole=roots(poly(Ag));
for k=1:length(Pole)
    plot(real(Pole(k)), imag(Pole(k)), 'bd');
end
a=0.01;
KI=a*Ksinv;
Ag=[A -B*KI; C zeros(m,m)];
Pole=roots(poly(Ag));
for k=1:length(Pole)
    plot(real(Pole(k)), imag(Pole(k)), 'bs');
end
a=0.02;
KI=a*Ksinv;
Ag=[A -B*KI; C zeros(m,m)];
Pole=roots(poly(Ag));
for k=1:length(Pole)
    plot(real(Pole(k)), imag(Pole(k)), 'b*');
end
hold off;
axis([-0.03 0.03 -0.03 0.03]);
xticks([-0.02 0 0.02]);
yticks([-0.02 0 0.02]);



%%    Einstellung des I-Anteils
%
%
figure(5);
clf;
T=[0:3:Te]';
for a=[0.005, 0.01, 0.02]
    KI=a*Ksinv;
    Ag=[A -B*KI;
        C zeros(m,m)];
    Bg=[zeros(n,m);
        -eye(m,m)];
    Cg=[C zeros(m,m)];
    Dg=zeros(m,m);
    IerweiterteAHLAnlageM=ss(Ag, Bg, Cg, Dg);
    Y=step(IerweiterteAHLAnlageM, T);
    subplot(2,2,1)
    plot(T, Y(:,1,1), 'b');
    hold on;
    subplot(2,2,2);
    plot(T, Y(:,1,2), 'b');
    hold on;
    subplot(2,2,3);
    plot(T, Y(:,2,1), 'b');
    hold on;
    subplot(2,2,4);
    plot(T, Y(:,2,2), 'b');
    hold on;
end;
subplot(2, 2, 1);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w11}$$');
rotateY;
yticks([0 1]);
xleer;
axis([0 Te -0.2 1.4]);
latextitle('F\"uhrungsverhalten ');
subplot(2, 2, 2);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w12}$$');
rotateY;
yticks([0 1]);
xleer;
axis([0 Te -0.2 1.4]);
latextitle('(a=0.005, 0.01, 0.02)');
subplot(2, 2, 3);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w21}$$');
rotateY;
yticks([0 1]);
xticks([0 100 200 300]);
latexxlabel('$$t$$ in min');
axis([0 Te -0.2 1.4]);
latextitle('');
subplot(2, 2, 4);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w22}$$');
rotateY;
yticks([0 1]);
xticks([0 100 200 300]);
latexxlabel('$$t$$ in min');
axis([0 Te -0.2 1.4]);
latextitle('');


%
%          Festlegung des I-Anteils:
KI=0.015*Ksinv;

%%     Wurzelortskurve
%
%
figure(6);
[A, B, C, D]=ssdata(AHLAnlageM);
Wurzelortskurve=ss([A -B*KI; C zeros(m,m)], [B; D], [C zeros(m,m)], zeros(m, m));
[Pole, Nullstellen]=pzmap(Wurzelortskurve);
for k=1:length(Pole)
    plot(real(Pole(k)), imag(Pole(k)), 'bx');
    hold on;
end
for k=1:length(Nullstellen)
    plot(real(Nullstellen(k)),imag(Nullstellen(k)),'bo');
end
k=1;
R=[];
for b=[0.0001:0.005:1, 1:0.01:8]
    KP=b*Ksinv;
    Ag=[A-B*KP*C -B*KI;
        C zeros(m,m)];
    d=roots(poly(Ag));
    R(k,:)=d.';
    k=k+1;
end;
plot(R, 'b.');
plot([-1 1],[0 0],'k:', [0 0],[-1 1],'k:');
axis([-0.15 0.02 -0.085 0.085]);
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('AHL-Anlage f\"ur P-Tuningfaktor');
hold off
xticks([-0.1 -0.05 0]);
yticks([-0.05 0 0.05]);



%  Detailuntersuchung der dominierenden Pole
%  Markieren der Pole für die später verwendeten Einstellungen
figure(7)
plot(R,'b.');
hold on
plot([-1 1],[0 0],'k:', [0 0],[-1 1],'k:');
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('b=0,5, 1, 2');
b=0.5;
KP=b*Ksinv;
Ag=[A-B*KP*C -B*KI;  C zeros(m,m)];
Pole=roots(poly(Ag));
for k=1:length(Pole)
    plot(real(Pole(k)), imag(Pole(k)), 'bd');
end
b=1;
KP=b*Ksinv;
Ag=[A-B*KP*C -B*KI;  C zeros(m,m)];
Pole=roots(poly(Ag));
for k=1:length(Pole)
    plot(real(Pole(k)), imag(Pole(k)), 'bs');
end
b=2;
KP=b*Ksinv;
Ag=[A-B*KP*C -B*KI;  C zeros(m,m)];
Pole=roots(poly(Ag));
for k=1:length(Pole)
    plot(real(Pole(k)), imag(Pole(k)), 'b*');
end
hold off;
axis([-0.08 0.02 -0.05 0.05]);
xticks([-0.05 0]);
yticks([-0.05 0 0.05]);





%%    Einstellung des P-Anteils
%
%
figure(8);
T=[0:0.7:Te]';
%  reiner P-Regler
b=0;
KP=b*Ksinv;
Ag=[A-B*KP*C -B*KI;
    C zeros(m,m)];
Bg=[B*KP;
    -eye(m,m)];
Cg=[C zeros(m,m)];
Dg=zeros(m,m);
PgeregelteAHLAnlageM=ss(Ag, Bg, Cg, Dg);
Y=step(PgeregelteAHLAnlageM, T);
subplot(2,2,1)
plot(T, Y(:,1,1), 'b--');
hold on;
subplot(2,2,2);
plot(T, Y(:,1,2), 'b--');
hold on;
subplot(2,2,3);
plot(T, Y(:,2,1), 'b--');
hold on;
subplot(2,2,4);
plot(T, Y(:,2,2), 'b--');
hold on;
for b=[0.5, 1, 2, 5]
    KP=b*Ksinv;
    Ag=[A-B*KP*C -B*KI;
        C zeros(m,m)];
    Bg=[B*KP;
        -eye(m,m)];
    Cg=[C zeros(m,m)];
    Dg=zeros(m,m);
    PgeregelteAHLAnlageM=ss(Ag, Bg, Cg, Dg);
    Y=step(PgeregelteAHLAnlageM, T);
    subplot(2,2,1)
    plot(T, Y(:,1,1),'b');
    subplot(2,2,2);
    plot(T, Y(:,1,2),'b');
    subplot(2,2,3);
    plot(T, Y(:,2,1),'b');
    subplot(2,2,4);
    plot(T, Y(:,2,2),'b');
end;
subplot(2, 2, 1);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w11}$$');
rotateY;
yticks([0 1]);
xleer;
axis([0 Te -0.2 1.2]);
latextitle('PI-geregelte AHL-Anlage');
subplot(2, 2, 2);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w12}$$');
rotateY;
yticks([0 1]);
xleer;
axis([0 Te -0.2 1.2]);
latextitle('(b=0, 0,5, 1, 2, 5)');
subplot(2, 2, 3);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w21}$$');
rotateY;
yticks([0 1]);
latexxlabel('$$t$$ in min');
axis([0 Te -0.2 1.2]);
latextitle('');
subplot(2, 2, 4);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w22}$$');
rotateY;
yticks([0 1]);
latexxlabel('$$t$$ in min');
axis([0 Te -0.2 1.2]);
latextitle('');



%%
%
%    Einstellung des P-Anteils (KP = diag(b1, 2) Ks^-1)
%
%
figure(9);
b2=2;
b1=2;
KP=[b1 0; 0 b2]*Ksinv;
Ag=[A-B*KP*C -B*KI;
    C zeros(m,m)];
Bg=[B*KP;
    -eye(m,m)];
Cg=[C zeros(m,m)];
Dg=zeros(m,m);
PgeregelteAHLAnlageM=ss(Ag, Bg, Cg, Dg);
Y=step(PgeregelteAHLAnlageM, T);
subplot(2,2,1)
plot(T, Y(:,1,1), 'b--');
hold on;
subplot(2,2,2);
plot(T, Y(:,1,2), 'b--');
hold on;
subplot(2,2,3);
plot(T, Y(:,2,1), 'b--');
hold on;
subplot(2,2,4);
plot(T, Y(:,2,2), 'b--');
hold on;
for b1=[0.5, 1, 2]
    KP=[b1 0; 0 b2]*Ksinv;
    Ag=[A-B*KP*C -B*KI;
        C zeros(m,m)];
    Bg=[B*KP;
        -eye(m,m)];
    Cg=[C zeros(m,m)];
    Dg=zeros(m,m);
    PgeregelteAHLAnlageM=ss(Ag, Bg, Cg, Dg);
    Y=step(PgeregelteAHLAnlageM, T);
    subplot(2,2,1)
    plot(T, Y(:,1,1),'b');
    subplot(2,2,2);
    plot(T, Y(:,1,2),'b');
    subplot(2,2,3);
    plot(T, Y(:,2,1),'b');
    subplot(2,2,4);
    plot(T, Y(:,2,2),'b');
end;
subplot(2, 2, 1);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w11}$$');
rotateY;
yticks([0 1]);
xleer;
axis([0 Te -0.2 1.2]);
latextitle('PI-geregelte AHL-Anlage');
subplot(2, 2, 2);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w12}$$');
rotateY;
yticks([0 1]);
xleer;
axis([0 Te -0.2 1.2]);
latextitle('(b1=0, 0,5, 1, 2)');
subplot(2, 2, 3);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w21}$$');
rotateY;
latexxlabel('$$t$$ in min');
yticks([0 1]);
axis([0 Te -0.2 1.2]);
latextitle('');
subplot(2, 2, 4);
plot([0 Te], [0 0], 'k:');
latexylabel('$$h_{\rm w22}$$');
rotateY;
yticks([0 1]);
latexxlabel('$$t$$ in min');
axis([0 Te -0.2 1.2]);
latextitle('');


%%     Übergangsfunktionsmatrix der Regelstrecke für die Identifikation
%      innerhalb der Projektaufgabe
%
%
figure(10);
Te=150;
T=[0:Te];
Y=step(AHLAnlageM, T);
F=rand(size(Y))-0.5*ones(size(Y));
Y(:,1,1)=Y(:,1,1)+0.3*F(:,1,1);
Y(:,1,2)=Y(:,1,2)+0.003*F(:,1,2);
Y(:,2,1)=Y(:,2,1)+0.03*F(:,2,1);
Y(:,2,2)=Y(:,2,2)+0.01*F(:,2,2);
subplot(2, 2, 1);
plot(T, Y(:,1,1), 'b');
latexylabel('$$\begin{array}{c}h_{11}\\ \mbox{in} \% \end{array}$$');
rotateY;
yticks([0 2]);
xleer;
axis([0 Te 0 2.5])
latextitle('\"Ubergangsfunktionsmatrix');
subplot(2, 2, 2);
plot(T, Y(:,1,2), 'b');
latexylabel('$$\begin{array}{c}h_{12}\\ \mbox{in} \% \end{array}$$');
rotateY;
yticks([0 0.04]);
xleer;
latextitle(' AHL-Anlage');
axis([0 Te 0 0.04])
subplot(2, 2, 3);
plot(T, Y(:,2,1), 'b');
latexylabel('$$\begin{array}{c}h_{21}\\ \mbox{in} \% \end{array}$$');
rotateY;
yticks([-0.3 0]);
latexxlabel('$$t$$ in min');
xticks([0 100]);
latextitle('');
axis([0 Te -0.3 0.05])
subplot(2, 2, 4);
plot(T, Y(:,2,2), 'b');
latexylabel('$$\begin{array}{c}h_{22}\\ \mbox{in} \% \end{array}$$');
rotateY;
yticks([0 0.1]);
latexxlabel('$$t$$ in min');
xticks([0 100]);
latextitle('');
axis([0 Te 0 0.1])



%%     Ausgabe der Bilder als .eps-Dateien
%
Epsbildc15('tunahl');
