%   Script DNAFUPROJ.M zur Herstellung von Bildern
%
%    Dezentrale FÜ-Regelung (DNA-Entwurf)
%    Lösung der Projektaufgabe
%
% verwendet BSPFUPROJ.M, EPSBILD.M
%
% J. Lunze
% für 3. Auflage: 25.3.2004
% für 5. Auflage: 18.10.2007
% für 7. Auflage 25.5.2012
% für 10. Auflage: 26.3.2020
%
%
clear
close all
%
%  Modell des Netzes aus Datei BSPFUPROJ.M
%
Bspfuproj;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Analyse des Netzes
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%   1. Netzverhalten bei Erhöhung der Leistungssollwerte der Kraftwerke
%
figure(1)
Te=50;                 % Endzeit der Simulation
Netz=ss(A, B, C, D);   % Netzverhalten bez\"uglich u_1, u_2
[Y,T]=step(Netz);
Y=50*Y;               % Sollwerterhöhung 50 MW
subplot(6,2,1);
plot(T, Y(:,1,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -30 60]);
yticks([0 50]);
latexylabel('$$f$$');  % in mHz
rotateY;
xleer;
latextitle('$$p_{\rm G_1 Soll}$$=50 MW');
subplot(6,2,2);
plot(T, Y(:,1,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -30 60]);
%xticks([]);
%yticks([]);
yleer;
xleer;
latextitle('$$p_{\rm G_2 Soll}$$=50 MW');
subplot(6,2,3);
plot(T, Y(:,2,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
yticks([0 50]);
latexylabel('$$p_{\rm G1}$$');    % in MW
rotateY;
xleer;
subplot(6,2,4);
plot(T, Y(:,2,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
%xticks([]);
%yticks([]);
yleer;
xleer;
subplot(6,2,5);
plot(T, Y(:,3,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
yticks([0 50]);
latexylabel('$$p_{\rm G2}$$');  % in MW
rotateY;
xleer;
subplot(6,2,6);
plot(T, Y(:,3,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
%xticks([]);
%yticks([]);
yleer;
xleer;
subplot(6,2,7);
plot(T, Y(:,4,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
yticks([0 50]);
latexylabel('$$p_{\rm V1}$$');  % in MW
rotateY;
xleer;
subplot(6,2,8);
plot(T, Y(:,4,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
%xticks([]);
%yticks([]);
yleer;
xleer;
subplot(6,2,9);
plot(T, Y(:,5,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
yticks([0 50]);
latexylabel('$$p_{\rm V2}$$');   % in MW
rotateY;
xleer;
subplot(6,2,10);
plot(T, Y(:,5,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
%xticks([]);
%yticks([]);
yleer;
xleer;
subplot(6,2,11);
plot(T, Y(:,6,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -40 40]);
yticks([0 50]);
xticks([0 20 40]);
latexylabel('$$p_{\rm t}$$');  % in MW
rotateY;
latexxlabel('$$t$$ in s');
subplot(6,2,12);
plot(T, Y(:,6,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -40 40]);
xticks([0 20 40]);
yleer;
latexxlabel('$$t$$ in s');



%%  2. Netzverhalten bei Lasterhöhung um 50 MW
%
figure(2)
Netzyd=ss(A, E, C, F);  % Netzverhalten Lasterhöhung -> Ausgänge
[Y,T]=step(Netzyd);
Y=50*Y;                % Laständerung 50MW
subplot(6,2,1);
plot(T, Y(:,1,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -60 30]);
yticks([-30 0 ]);
latexylabel('$$f$$');  % in mHz
rotateY;
xleer;
latextitle('$$p_{\rm L_1}$$=50 MW');
subplot(6,2,2);
plot(T, Y(:,1,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -60 30]);
%xticks([]);
%yticks([]);
yleer;
xleer;
latextitle('$$p_{\rm L_2}$$=50 MW');
subplot(6,2,3);
plot(T, Y(:,2,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
yticks([0 50]);
latexylabel('$$p_{\rm G1}$$');
rotateY;
xleer;
subplot(6,2,4);
plot(T, Y(:,2,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
%xticks([]);
%yticks([]);
yleer;
xleer;
subplot(6,2,5);
plot(T, Y(:,3,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
yticks([0 50]);
latexylabel('$$p_{\rm G2}$$');
rotateY;
xleer;
subplot(6,2,6);
plot(T, Y(:,3,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
%xticks([]);
%yticks([]);
yleer;
xleer;
subplot(6,2,7);
plot(T, Y(:,4,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -60 20]);
yticks([-50 0 50]);
latexylabel('$$p_{\rm V1}$$');
rotateY;
xleer;
subplot(6,2,8);
plot(T, Y(:,4,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -60 20]);
%xticks([]);
%yticks([]);
yleer;
xleer;
subplot(6,2,9);
plot(T, Y(:,5,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -60 20]);
yticks([-50 0 50]);
latexylabel('$$p_{\rm V2}$$');
rotateY;
xleer;
subplot(6,2,10);
plot(T, Y(:,5,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -60 20]);
%xticks([]);
%yticks([]);
yleer;
xleer;
subplot(6,2,11);
plot(T, Y(:,6,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -40 40]);
xticks([0 20 40]);
yticks([-40 0 50]);
latexylabel('$$p_{\rm t}$$');
rotateY;
latexxlabel('$$t$$ in s');
subplot(6,2,12);
plot(T, Y(:,6,2), 'b', [0 Te],[0 0],'k:');
axis([0 Te -40 40]);
xticks([0 20 40]);
yleer;
latexxlabel('$$t$$ in s');



%%   3. Netzkennlinienfehler bei Veränderung der Leistungssollwerte der Kraftwerke
%
figure(3)
Netzeu=ss(A, B, CE, DE);   % Netzverhalten Stellgrößen -> Netzkennlinienfehler
[Y,T]=step(Netzeu);
Y=50*Y;                % Sollwertänderung 50MW im Netz 1
subplot(6,1,[1,2]);
plot(T, Y(:,1,1), 'b', [0 Te],[0 0],'k:', [0 Te],[50 50],'k:');
axis([0 Te -10 70]);
yticks([0 50]);
latexylabel('$$e_1$$ in MW');
rotateY;
xleer;
%latexxlabel('$$t$$ in s');
latextitle('Netzkennlinienfehler bei $$p_{\rm G1 Soll}$$=50 MW');
subplot(6,1,[3,4]);
plot(T, Y(:,2,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
xticks([0 20 40 60]);
yticks([0 50]);
latexylabel('$$e_2$$ in MW');
rotateY;
%latextitle('Netzverhalten bei p_G_2_ _S_o_l_l=50 MW');
latexxlabel('$$t$$ in s');


%%  4. Netzkennlinienfehler bei Lasterhöhung um 50 MW
%
figure(4)
Netzed=ss(A, E, CE, FE);    % Netzverhalten Lasterhöhung -> Netzkennlinienfehler
[Y,T]=step(Netzed);
Y=50*Y;                % Laständerung 50MW im Netz 1
subplot(6,1,[1,2]);
plot(T, Y(:,1,1), 'b', [0 Te],[0 0],'k:', [0 Te],[-50 -50],'k:');
axis([0 Te -70 10]);
yticks([-50 0]);
latexylabel('$$e_1$$ in MW');
rotateY;
xleer;
%latexxlabel('$$t$$ in s');
latextitle('Netzkennlinienfehler bei $$p_{\rm L_1}$$=50 MW');
subplot(6,1,[3,4]);
plot(T, Y(:,2,1), 'b', [0 Te],[0 0],'k:');
axis([0 Te -20 60]);
xticks([0 20 40 60]);
yticks([0 50]);
latexylabel('$$e_2$$ in MW');
rotateY;
latexxlabel('$$t$$ in s');


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Entwurf des Reglers für das Netz 1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(5);
%    Ausblenden des Teilsystems 1 (A, Bs1, Es1, Cs1, Ds1, Fs1)
%    aus Gesamtmodell  (vgl. BSPFUPROJ.M)
%    Eingangsgröße: u_1 (Stellgröße des Kraftwerks 1)
%    Ausgangsgröße: e_1 (Netzkennlinienfehler des Netzes 1)
Bs1=B(1:n, 1);
Es1=E(1:n, 1);
Cs1=CE(1,:);
Ds1=0;
Fs1=FE(1,1);
Teilsystem1=ss(A, Bs1, Cs1, 0);
%  Bodediagramm der Strecke und des gewählten Reglers
[Betrag, Phase, Omega]=bode(Teilsystem1);
subplot(4,1,[1,2]);;
semilogx(Omega, 20*log10(Betrag(:)), 'b');
grid on;
hold on;
%plot([0.01, 10], [20, -40 ],'b--');   % Gerade mit Neigung -20 dB/Dekade
%
KI1=0.09;                            % PI-Regler
KP1=0.09;
omega0=KI1/KP1;
BetragK=20*log10(sqrt(2)*KP1);
%  Amplitudengang des PI-Reglers
plot([0.01*omega0 omega0 10*omega0],[40+BetragK BetragK BetragK],'b--');
xleer;
latexylabel('$$|G_{11}|$$ in dB');
rotateY;
yticks([-40 -20 0 20]);
axis([0.01 10 -40 20]);
%xticks([]);
%yticks([]);
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 -180 10]);
xticks([0.001 0.01 0.1 1 10 100]);
yticks([-90 0]);
grid on


%%
%
%   Bodediagramm der offenen Kette des Regelkreises 1
%
%
figure(6);
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;
hold on;
plot([0.01, 10], [20, -40 ],'b--');   % Gerade mit Neigung -20 dB/Dekade
xleer;
latexylabel('$$|G_0|$$ in dB');
rotateY;
latextitle('Offene Kette des Teilregelkreises 1');
axis([0.01 10 -40 20]);
yticks([-40 -20 0 20]);
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 10]);
xticks([0.001 0.01 0.1 1 10 100]);
yticks([-90 0]);
grid on



%%
%
%    Verhalten des ersten Regelkreises bei Lastanstieg um 50 MW
%
%
figure(7);
Ag1=[A-Bs1*KP1*Cs1  -Bs1*KI1
     Cs1             0      ];
Eg1=[Es1-Bs1*KP1*Fs1
     Fs1            ];
Cg1=[C(1, 1:n) 0;              %  y_1 Frequenz
     C(6, 1:n) 0;              %  y_2 Übergabeleistung p_t_1
     Cs1       0];             %  e_1 Netzkennlinienfehler
Fg1=[F(1, 1)
     Fs1
     Fs1];
Regelkreis1=ss(Ag1, Eg1, Cg1, Fg1);
[Y, T]=step(Regelkreis1);
Y=50*Y;
subplot(6,1,[1,2]);
plot(T, Y(:, 1), 'b', [0 50], [0 0],'k:');
xleer;
latexylabel('$$f$$ in mHz');
rotateY;
latextitle('St\"orverhalten bei Regelung von Netz 1 f\"ur $$p_{\rm L_1}$$=50 MW');
axis([0 Te -60 30]);
yticks([-60 0]);
subplot(6,1,[3,4]);
plot(T, Y(:, 2), 'b', [0 50], [0 0],'k:');
xleer;
latexylabel('$$p_{\rm t}$$ in MW');
rotateY;
yticks([-50 0]);
axis([0 Te -60 20]);
%xticks([]);
%yticks([]);
subplot(6,1,[5,6]);
plot(T, Y(:, 3), 'b', [0 50], [0 0],'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$e_1$$ in MW');
rotateY;
yticks([-50 0]);
xticks([0 20 40 60]);
axis([0 Te -60 20]);


%%  Entwurf des Reglers für das Teilnetz 2
%
figure(8)
%    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=CE(2,:);
Ds2=0;
Fs2=FE(2,2);
Teilsystem2=ss(A, Bs2, Cs2, Ds2);
% Bodediagramm der Regelstrecke und des gewählten Reglers
[Betrag, Phase, Omega]=bode(Teilsystem2);
subplot(4,1,[1,2]);;
semilogx(Omega, 20*log10(Betrag(:)), 'b');
hold on;
%plot([0.01, 10], [20, -40 ],'b--');   % Gerade mit Neigung -20 dB/Dekade
%
KI2=0.05;                            % PI-Regler
KP2=0.05;
omega0=KI2/KP2;
BetragK=20*log10(sqrt(2)*KP2);
% Amplitudengang des PI-Reglers
plot([0.01*omega0 omega0 10*omega0],[40+BetragK BetragK BetragK],'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 -180 10]);
xticks([0.001 0.01 0.1 1 10 100]);
yticks([-90 0]);
grid on



%%
%
%   Bodediagramm der offenen Kette des Regelkreises 2
%
%
figure(9);
A02=[A  zeros(n,1)
     Cs2 0        ];
B02=[Bs2
      0  ];
C02=[KP2*Cs2 KI2];
D02=0;
offeneKette2=ss(A02, B02, C02, D02);
[Betrag, Phase, Omega]=bode(offeneKette2);
subplot(4,1,[1,2]);;
semilogx(Omega, 20*log10(Betrag(:)), 'b');
grid on;
xleer;
latexylabel('$$|G_0|$$ in dB');
rotateY;
latextitle('Offene Kette des Teilregelkreises 2');
axis([0.01 10 -60 20]);
yticks([-40 -20 0 20]);
%
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 10 -270 10]);
xticks([0.001 0.01 0.1 1 10 100]);
yticks([-90 -45 0]);
grid on



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



%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Überpr\"ufung der Diagonaldominanz
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(11);
% Die Diagonaldominanz der R\"uckf\"uhrdifferenzmatrix der offenen Kette mit den
% beiden dezentralen Teilreglern wird gepr\"uft
%
% I-erweiterte Strecke
AI=[A  zeros(n, 2)
    CE zeros(2, 2)];
BI=[B
    zeros(2,2)];
CI=[CE zeros(2,2)
    zeros(2,n) eye(2,2)];
DI=eye(4,2);
%  Dezentraler Regler
K=[KP1 0  KI1  0;
    0  KP2 0  KI2];
KP=K(1:2, 1:2);
KI=K(1:2, 2:4);
A0=AI;
B0=BI;
C0=K*CI;
D0=eye(2,2);
W=logspace(-2, 1, 50);
Wsize=max(size(W));
%  Bildung der R\"uckf\"uhrdifferenzfunktion
% (diese R\"uckf\"uhrdifferenzfunktion ist die Übertragungsfunktionsmatrix
% des folgenden Systems)
Rueckfuehrdifferenz=ss(A0, B0, C0, D0+eye(2,2));
[FBetrag, FPhase, W]=bode(Rueckfuehrdifferenz, W);
%     Überpr\"ufung der Zeilendominanz, 1. Zeile
subplot(4,1,[1,2]);
F11=20*log10(FBetrag(1,1,:));
F12=20*log10(FBetrag(1,2,:));
semilogx(W, F11(:), 'b', W, F12(:),'b--');
axis([0.01 10 -80 30]);
%xticks([]);
%yticks([]);
grid on
latexylabel('$$|F_{11}|, |F_{12}|$$');
rotateY;
yticks([-40 0]);
xleer;
latextitle('Zeilendominanz der R\"uckf\"uhrdifferenzmatrix');
%     Überpr\"ufung der Zeilendominanz, 2. Zeile
subplot(4,1,[3,4]);
F21=20*log10(FBetrag(2,1,:));
F22=20*log10(FBetrag(2,2,:));
semilogx(W, F22(:), 'b', W, F21(:),'b--');
axis([0.01 10 -80 30]);
%xticks([]);
%yticks([]);
grid on
latexylabel('$$|F_{21}|, |F_{22}|$$');
rotateY;
yticks([-80 0]);
latexxlabel('$$\omega$$ in $${\rm rad/s}$$');  %  in $$\frac{\rm rad}{\rm s}$$
xticks([0.001 0.01 0.1 1 10 100]);
yticks([-40 0]);
grid on


%%    Verhalten des dezentral geregelten Systems bei Laständerung
%
%
figure(12)
% Modell des Regelkreises mit den beiden Teilreglern
FI=[   FE
    zeros(2,2)];
EI=[E;
    FE];
Ag=AI-BI*K*CI;
Eg=[E-B*KP*FE
   FE];
Cg=[C(1, 1:n) zeros(1,2)
    C(6, 1:n) zeros(1,2)];    % Ausgänge: f, p_t
Fg=[F(1, :);
    F(6, :)];
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, 50*F1(:, 1), 'b', [0 Te], [0 0],'k:');
xleer;
latexylabel('$$f$$ in mHz');
rotateY;
axis([0 Te -50 30]);
yticks([-50 0]);
latextitle('Energienetz mit dezentraler PI-Regelung');
subplot(6,1,[3,4]);
plot(T, 50*F1(:, 2), 'b', [0 Te], [0 0],'k:');
xleer;
latexylabel('$$p_{\rm t}$$ in MW');
rotateY;
yticks([-50 0]);
axis([0 Te -60 20]);
yticks([-50 0]);
subplot(6,1,[5,6]);
plot(T, 50*U(:, 1, 1), 'b', T, 50*U(:, 2,1), 'b', [0 Te], [0,0], 'k:', [0 Te],[50,50],'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$u_1, u_2$$ in MW');
rotateY;
yticks([0 50]);
axis([0 Te -20 60]);
xticks([0 20 40 60]);
yticks([0 50]);



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