%   Script DNAFU.M zur Herstellung von Bildern
%
%    Dezentrale FÜ-Regelung (DNA-Entwurf)
%
% verwendet BSPFUREG.M, EPSBILD.M
% Funktion lambdaP, gerkreis
%
% J. Lunze
% 16.10.1996
% für 2. Auflage: 15.12.2001
% für 3. Auflage: 15.3.2004
% für 5. Auflage: 18.10.2007
% für 7. Auflage 25.5.2012
% für 10. Auflage: 26.3.20120
%
clear
close all
Bspfureg;
CE1=C1;
CE2=C2;
CE3=C3;

%%
%
%   Regler für das Netz 1
%
%
figure(1);
%    Ausblenden des Teilsystems 1 (A, Bs1, Es1, Cs1, Ds1, Fs1)
%    aus Gesamtmodell  (vgl. BSPFU.M)
Bs1=B(1:n, 1);
Es1=E(1:n, 1);
Cs1=CE1;
Ds1=0;
Fs1=FE1(1,1);
Teilsystem1=ss(A, Bs1, Cs1, Ds1);
[Betrag, Phase, Omega]=bode(Teilsystem1);
subplot(4,1,[1,2]);;
semilogx(Omega, 20*log10(Betrag(:)), 'b');
grid on;
hold on;
plot([0.03, 30], [20, -40 ],'b--');
xleer;
latexylabel('$$|G_{11}|$$ in dB');
rotateY;
axis([0.01 10 -40 20]);
yticks([-40 -20 0 20 40]);
latextitle('Bodediagramm von Teilnetz 1');
hold off;
%
subplot(4,1,[3,4]);;
semilogx(Omega, Phase(:), 'b');
latexxlabel('$$\omega$$ in $${\rm rad/s}$$'); %  in $$\frac{\rm rad}{\rm s}$$
latexylabel('$${\rm arg}(G_{11})$$');
rotateY;
axis([0.01 10 -90 0]);
xticks([0.01 0.1 1 10 100]);
yticks([-90 -45 0]);
grid on


%%   Offene Kette des Regelkreises 1
%
%
figure(2);
KI1=3.33;
KP1=1;
A01=[A  zeros(n,1)
     Cs1 0        ];
B01=[Bs1
      0  ];
C01=[KP1*Cs1 KI1];
D01=0;
Regelkreis1=ss(A01, B01, C01, D01);
[Betrag, Phase, Omega]=bode(Regelkreis1);
subplot(4,1,[1,2]);;
semilogx(Omega, 20*log10(Betrag(:)), 'b');
grid on;
xleer;
latexylabel('$$|G_0|$$ in dB');
rotateY;
latextitle('Offene Kette des Teilreglers 1');
axis([0.01 10 -40 50]);
yticks([-40 -20 0 20 40]);
%
subplot(4,1,[3,4]);;
semilogx(Omega, Phase(:), 'b');
latexylabel('$${\rm arg}(G_0)$$');
rotateY;
latexxlabel('$$\omega$$ in $${\rm rad/s}$$'); %  in $$\frac{\rm rad}{\rm s}$$
axis([0.01 10 -180 0]);
xticks([0.01 0.1 1 10 100]);
yticks([-180 -135 -90 -45 0]);
grid on;



%%    Verhalten des ersten Regelkreises bei Lastabfall von 100 MW
%
%
figure(3);
KI1=KI1/20;
KP1=KP1/20;
Ag1=[A-Bs1*KP1*Cs1  -Bs1*KI1
     Cs1             0      ];
Eg1=[Es1-Bs1*KP1*Fs1
     Fs1            ];
Cg1=[C(1, 1:n) 0;
     C(2, 1:n) 0];
Fg1=[F(1, 1)
     F(2, 1)];
Regelkreis1=ss(Ag1, Eg1, Cg1, Fg1);
[Y, T]=step(Regelkreis1);
subplot(4,1,[1,2]);;
plot(T, 100*Y(:, 1), 'b', [0 50], [0 0],'k:');
xleer;
latexylabel('$$f$$ in mHz');
rotateY;
latextitle('St\"orverhalten bei Regelung von Netz 1');
axis([0 25 -40 20]);
yticks([-40 -20 0 20]);
%
subplot(4,1,[3,4]);;
plot(T, 100*Y(:, 2), 'b', [0 50], [0 0],'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$p_{\rm t1}$$ in MW');
rotateY;
axis([0 25 -80 20]);
xticks([0 10 20 30]);
yticks([-80 -40 0]);


%%  Regler für Teilnetz 2
%
%
figure(4)
%    Ausblenden des Teilsystems 2 (A, Bs2, Es2, Cs2, Ds2, Fs2)
%    aus Gesamtmodell  (vgl. BSPFU.M)
Bs2=B(1:n, 2);
Es2=E(1:n, 2);
Cs2=CE2;
Ds2=0;
Fs2=FE(2,2);
Teilsystem2=ss(A, Bs2, Cs2, Ds2);
[Betrag, Phase, Omega]=bode(Teilsystem2);
subplot(4,1,[1,2]);;
semilogx(Omega, 20*log10(Betrag(:)), 'b');
hold on;
plot([0.015, 15], [20, -40 ],'b--');
xleer;
grid on;
latexylabel('$$|G_{22}|$$ in dB');
rotateY;
axis([0.01 10 -40 20]);
yticks([-40 -20 0 20]);
latextitle('Bodediagramm von Teilnetz 2');
hold off;
subplot(4,1,[3,4]);;
semilogx(Omega, Phase(:), 'b');
latexxlabel('$$\omega$$ in $${\rm rad/s}$$'); %  in $$\frac{\rm rad}{\rm s}$$
latexylabel('$${\rm arg}(G_{22})$$');
rotateY;
axis([0.01 10 -100 20]);
xticks([0.01 0.1 1 10 100]);
yticks([-90 -45 0]);
grid on;


%%   Offene Kette des Regelkreises 2
%
%
figure(5);
KI2=6.667/60;
KP2=1/60;
A02=[A  zeros(n,1)
     Cs2 0        ];
B02=[Bs2
      0  ];
C02=[KP2*Cs2 KI2];
D02=0;
Regelkreis2=ss(A02, B02, C02, D02);
[Betrag, Phase, Omega]=bode(Regelkreis2);
subplot(4,1,[1,2]);;
semilogx(Omega, 20*log10(Betrag(:)), 'b');
grid on;
xleer;
latexylabel('$$|G_0|$$ in dB');
rotateY;
latextitle('Offene Kette des Teilreglers 2');
yticks([-100 -50 0 50]);
%
subplot(4,1,[3,4]);;
semilogx(Omega, Phase(:), 'b');
grid on;
latexylabel('$${\rm arg}(G_0)$$');
rotateY;
latexxlabel('$$\omega$$ in $${\rm rad/s}$$'); %  in $$\frac{\rm rad}{\rm s}$$
axis([0.01 1000 -190 -80]);
xticks([0.01 0.1 1 10 100]);
yticks([-180 -135 -90 -45 0]);


%%   Verhalten des zweiten Regelkreises bei Lastabfall von 100 MW
%
%
figure(6);
Ag2=[A-Bs2*KP2*Cs2  -Bs2*KI2
     Cs2             0      ];
Eg2=[Es2-Bs2*KP2*Fs2
     Fs2            ];
Cg2=[C(1, 1:n) 0;
     C(3, 1:n) 0];
Fg2=[F(1, 2)
     F(3, 2)];
Regelkreis2=ss(Ag2, Eg2, Cg2, Fg2);
[Y, T]=step(Regelkreis2);
subplot(4,1,[1,2]);;
plot(T, 100*Y(:, 1), 'b', [0 50], [0 0],'k:');
xleer;
latexylabel('$$f$$ in mHz');
rotateY;
latextitle('St\"orverhalten bei Regelung von Netz 2');
axis([0 25 -40 20]);
yticks([-40 -20 0 20]);
%
subplot(4,1,[3,4]);;
plot(T, 100*Y(:, 2), 'b', [0 50], [0 0],'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$p_{\rm t2}$$ in MW');
rotateY;
axis([0 25 -80 20]);
xticks([0 10 20 30]);
yticks([-80 -40 0]);


%%    Regler für Netz 3
%
%
figure(7)
%    Ausblenden des Teilsystems 3 (A, Bs3, Es3, Cs3, Ds3, Fs3)
%    aus Gesamtmodell  (vgl. BSPFU.M)
Bs3=B(1:n, 3);
Es3=E(1:n, 3);
Cs3=CE3;
Ds3=0;
Fs3=FE(3,3);
Teilsystem3=ss(A, Bs3, Cs3, Ds3);
[Betrag, Phase, Omega]=bode(Teilsystem3);
subplot(4,1,[1,2]);;
semilogx(Omega, 20*log10(Betrag(:)), 'b');
grid on;
hold on;
plot([0.015, 15], [20, -40 ],'b--');
xleer;
latexylabel('$$|G_{33}|$$ in dB');
rotateY;
axis([0.01 10 -40 20]);
yticks([-40 -20 0 20]);
latextitle('Bodediagramm von Teilnetz 3');
hold off;
subplot(4,1,[3,4]);;
semilogx(Omega, Phase(:), 'b');
grid on;
latexxlabel('$$\omega$$ in $${\rm rad/s}$$'); %  in $$\frac{\rm rad}{\rm s}$$
latexylabel('$${\rm arg}(G_{33})$$');
rotateY;
axis([0.01 10 -100 0]);
xticks([0.01 0.1 1 10 100 1000]);
yticks([-180 -135 -90 -45 0]);


%%  Offene Kette des Regelkreises 3
%
%
figure(8);
KI3=6.667/60;
KP3=1/60;
A03=[A  zeros(n,1)
     Cs3 0        ];
B03=[Bs3
      0  ];
C03=[KP3*Cs3 KI3];
D03=0;
Regelkreis3=ss(A03, B03, C03, D03);
[Betrag, Phase, Omega]=bode(Regelkreis3);
subplot(4,1,[1,2]);;
semilogx(Omega, 20*log10(Betrag(:)), 'b');
grid on;
xleer;
latexylabel('$$|G_0|$$ in dB');
rotateY;
latextitle('Offene Kette des Teilreglers 3');
yticks([-100 -50 0 50]);
subplot(4,1,[3,4]);;
semilogx(Omega, Phase(:), 'b');
grid on;
latexylabel('$${\rm arg}(G_0)$$');
rotateY;
latexxlabel('$$\omega$$ in $${\rm rad/s}$$'); %  in $$\frac{\rm rad}{\rm s}$$
axis([0.01 1000 -180 -90]);
xticks([0.01 0.1 1 10 100 1000]);
yticks([-180 -135 -90 -45 0]);


%%   Verhalten des dritten Regelkreises bei Lastabfall von 100 MW
%
%
figure(9);
Ag3=[A-Bs3*KP3*Cs3  -Bs3*KI3
     Cs3             0      ];
Bg3=[Es3-Bs3*KP3*Fs3
     Fs3            ];
Cg3=[C(1, 1:n) 0;
     C(4, 1:n) 0];
Dg3=[F(1, 3)
     F(4, 3)];
Lastabfall3=ss(Ag3, Bg3, Cg3, Dg3);
[Y, T]=step(Lastabfall3);
subplot(4,1,[1,2]);;
plot(T, 100*Y(:, 1), 'b', [0 50], [0 0],'k:');
xleer;
latexylabel('$$f$$ in mHz');
rotateY;
latextitle('St\"orverhalten bei Regelung von Netz 3');
axis([0 25 -40 20]);
yticks([-40 -20 0 20]);
%
subplot(4,1,[3,4]);;
plot(T, 100*Y(:, 2), 'b', [0 50], [0 0],'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$p_{\rm t3}$$ in MW');
rotateY;
axis([0 25 -80 20]);
xticks([0 10 20 30]);
yticks([-80 -40 0]);



%%   Überpr\"ufung der Diagonaldominanz
%
%
figure(10);
% I-erweiterte Strecke
AI=[A  zeros(n, 3)
    CE zeros(3, 3)];
BI=[B
    zeros(3,3)];
CI=[CE zeros(3,3)
    zeros(3,n) eye(3,3)];
DI=eye(6,3);
%  Dezentraler Regler
K=[KP1 0   0  KI1  0   0;
    0  KP2 0   0  KI2  0;
    0  0  KP3  0   0  KI3];
KP=K(1:3, 1:3);
KI=K(1:3, 4:6);
A0=AI;
B0=BI;
C0=K*CI;
D0=eye(3,3);
W=logspace(-2, 1, 100);
Wsize=max(size(W));
Strecke=ss(A0, B0, C0, D0);
[Betrag, Phase]=bode(Strecke, W);
G=Betrag.*cos(Phase*pi/180) + Betrag.*sin(Phase*pi/180)*i;
for k=1:Wsize
   FF=[G(1,1,k) G(1,2,k) G(1,3,k);
       G(2,1,k) G(2,2,k) G(2,3,k);
       G(3,1,k) G(3,2,k) G(3,3,k)];
   detF(k) = det(FF);
   PMatrix = inv(diag([G(1,1,k) G(2,2,k) G(3,3,k)]))*(FF- diag(diag(FF)));
   lP(k) = lambdaP(abs(PMatrix));                                             % function lambdaP
end;
plot(real(detF), imag(detF), 'b', real(detF), -imag(detF), 'b--', [0 0], [-12 12], 'k:', [-7 7], [0 0], 'k:');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Ortskurve von $$\det F(j{\omega}))$$');
axis([-7 5 -6 6]);
xticks([-5 0 5]);
yticks([-5 0 5]);
axis('square');


figure(11);
subplot(4,1,[1,2]);
semilogx(W, lP, 'b');
latexxlabel('$$\omega$$ in $${\rm rad/s}$$'); %  in $$\frac{\rm rad}{\rm s}$$
latexylabel('$$\lambda_{\rm P}$$');
rotateY;
latextitle('Perronwurzel f\"ur Pr\"ufung der Diagonaldominanz');
xticks([0.01 0.1 1 10 100]);
yticks([0 0.1 0.2 0.3]);


figure(12);
Real=real(G(1,1,:)-ones(1,1,Wsize));
Imag=imag(G(1,1,:));
plot(Real(:), Imag(:), 'b', [0,0], [-5,2], 'k:', [-2, 10], [0,0], 'k:');
hold on;
R1(1,1,:)=Betrag(1,2,:)+Betrag(1,3,:);  %  Zeilendominanz
for k=1:2:Wsize
   gerkreis(G(1,1,k)-1, R1(1,1,k))   % function gerkreis
end;
plot([-1 -1 ],[-0.05 0.05], 'k');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Gershgorinkreise f\"ur F\"U-Regelung (TR 1)');
axis([-1.5 .5 -1.5 0.5]);
xticks([-1 0 1]);
yticks([-1 0 1]);
axis('square');
hold off;



figure(13);
Real=real(G(2,2,:)-ones(1,1,Wsize));
Imag=imag(G(2,2,:));
plot(Real(:), Imag(:), 'b', [0,0], [-5,2], 'k:', [-2, 10], [0,0], 'k:');
hold on;
R2(1,1,:)=Betrag(2,1,:)+Betrag(2,3,:);  %  Zeilendominanz
for k=1:2:Wsize
   gerkreis(G(2,2,k)-1, R2(1,1,k))   % function gerkreis
end;
plot([-1 -1 ],[-0.05 0.05], 'k');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Gershgorinkreise f\"ur F\"U-Regelung (TR 2)');
axis([-1.5 .5 -1.5 0.5]);
xticks([-1 0 1]);
yticks([-1 0 1]);
axis('square');
hold off;


figure(14);
Real=real(G(3,3,:)-ones(1,1,Wsize));
Imag=imag(G(3,3,:));
plot(Real(:), Imag(:), 'b', [0,0], [-5,2], 'k:', [-2, 10], [0,0], 'k:');
hold on;
%xticks([]);
%yticks([]);
R3(1,1,:)=Betrag(3,1,:)+Betrag(3,2,:);  %  Zeilendominanz
for k=1:2:Wsize
   gerkreis(G(3,3,k)-1, R3(1,1,k))   % function gerkreis
end;
plot([-1 -1 ],[-0.05 0.05], 'b');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Gershgorinkreise f\"ur F\"U-Regelung (TR 3)');
axis([-1.5 .5 -1.5 0.5]);
xticks([-1 0 1]);
yticks([-1 0 1]);
axis('square');
hold off;



%%    Verhalten des dezentral geregelten Systems
%
%
figure(15)
CI=[CE         zeros(3,3)
    zeros(3,n) eye(3,3)];
FI=[   FE
    zeros(3,3)];
EI=[E;
    FE];
Ag=AI-BI*K*CI;
Eg=[E-B*KP*FE
   FE];
Cg=[C(1, 1:n) zeros(1,3)
    C(2, 1:n) zeros(1,3)];
 Fg=F(1:2, :);
Regelsys1=ss(Ag, Eg, Cg, Fg);
[F1, T]=step(Regelsys1, T);
Regelsys2=ss(Ag, Eg, -K*CI, -K*FI);
U=step(Regelsys2, T);
subplot(6,1,[1,2]);
plot(T, 100*F1(:, 1), 'b', [0 40], [0 0],'k:');
xleer;
latexylabel('$$f$$ in mHz');
rotateY;
axis([0 25 -40 20]);
yticks([-40 -20 0 20]);
latextitle('Energienetz mit dezentraler PI-Regelung');
subplot(6,1,[3,4]);
plot(T, 100*F1(:, 2), 'b', [0 40], [0 0],'k:');
xleer;
latexylabel('$$p_{\rm t1}$$ in MW');
rotateY;
axis([0 25 -80 20]);
yticks([-80 -40 0]);
subplot(6,1,[5,6]);
plot(T, 100*U(:, 1), 'b', T, 100*U(:, 2), 'b', T, 100*U(:, 3), 'b', [0 40], [0,0], 'k:', [0 40],[100,100],'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$u_i$$ in MW');
rotateY;
axis([0 25 -20 170]);
xticks([0 10 20 30]);
yticks([0 50 100 150]);


%%    Verlauf der Generatorleistungen
%
%
figure(16);
Cg=[C(5:8, 1:10) zeros(4,3)];
Fg=F(5:8, 1:m);
Regelsys1=ss(Ag, Eg, Cg, Fg);
[G, T]=step(Regelsys1,T);
subplot(4,1,[1,2]);;
plot(T, 100*G(:,1,1), 'b', T, 100*G(:,2,1), 'b', T, 100*G(:,3,1), 'b', T, 100*G(:,4,1), 'b', [0 40], [0,0], 'k:');
xleer;
latexylabel('$$p_{{\rm G}i}$$ in MW');
rotateY;
axis([0 25 -20 100]);
yticks([ 0 50 100]);
latextitle('Generatorleistungen bei dezentraler PI-Regelung');
subplot(4,1,[3,4]);;
plot(T, 100*U(:,1,1), 'b', T, 100*(G(:,1,1)+G(:,2,1)+G(:,3,1)+G(:,4,1)), 'b', [0 40], [0,0], 'k:', [0 40],[100,100],'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$u_1, p_{\rm G}$$ in MW');
rotateY;
latextext(5, 120, '$$u_1$$');
latextext(13, 80, '$$p_{\rm G}$$');
axis([0 25 0 170]);
xticks([0 10 20 30]);
yticks([0 50 100 150]);

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