%   Script BENGLIF.M zur Herstellung von Bildern
%
%   Dezentrale FÜ-Regelung nach Bengtsson/Lindahl
%
%  verwendet BSPFUREG.M, EPSBILD.M
%  Normierung der Bilder wie in DNAFU.M
%
% J. Lunze
% 14.9.1996
% für 2. Auflage: 16.12.2001
% für 3. Auflage: 15.3.2004
% für 5. Auflage: 18.10.2007
% für 7. Auflage: 31.3.2012
% für 8. Auflage: 28.3.2014 
% für 10. Auflage: 7.9.2019
%
%
clear
close all
Bspfureg;   % definiert A, B, C, D, E, F, CE, EE, FE


% I-erweiterte Strecke
AI=[A  zeros(n, 3)
    CE zeros(3, 3)];
BI=[B
    zeros(3,3)];
CI=[CE         zeros(3,3)    % Ausgang = Vektor der drei Netzkennlinienfehler
    zeros(3,n) eye(3,3)];    % Integratorausgänge der drei Netze
DI=zeros(6,3);
EI=[E;
    FE];
FI=[FE
    zeros(3,3)];
%%
%
%   Optimalregler für das Netz
%
%
figure(1);
CQ = [C(1:4, 1:10)  zeros(4,3)  % bewertet Frequenz, 3 Übergabeleistungen
      zeros(3,10)   eye(3,3)];  % bewertete 3 I-Anteile
Q = CQ'*diag([k1+k2+k3, 1, 1, 1, 1, 1, 1])*CQ;
R=10*eye(3,3);
K = lqr(AI, BI, Q, R);
%   Regelkreis mit Störung durch Lastanstieg um 100MW
Ag=AI-BI*K;
Bg=EI;
Cg=[C(1, 1:10) zeros(1,3)];  % Ausgang = Frequenz
Dg=F(1, 1:m);
Regelkreis1=ss(Ag, Bg, Cg, Dg);
RegelkreisStellgr=ss(Ag, Bg, -K, zeros(3,3));
[F1, T]=step(Regelkreis1);
U=step(RegelkreisStellgr, T);
subplot(6,1,[1,2]);
plot(T, 100*F1(:,1,1),'b',  [0 30], [0 0],'k:');
xleer;
latexylabel('$$f$$ in mHz');
rotateY;
axis([0 25 -40 20]);
latextitle('Energienetz mit PI-Zustandsr\"uckf\"uhrung');
yticks([-40 -20 0 20]);
%
subplot(6,1,[3,4]);
plot(T, 100*U(:, 1, 1),'b', T, 100*U(:, 2, 1),'b', T, 100*U(:, 3, 1),'b', [0 25], [100 100], 'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$u$$ in MW');
rotateY;
axis([0 25 -10 170]);
xticks([0 10 20]);
yticks([0 50 100 150]);




%%
%
%    Näherung der ZRF durch dezentrale ARF
%
echo  off
figure(2);
NetzZRF=ss(AI-BI*K, BI, zeros(1, n+3), zeros(1,3));
P=pzmap(NetzZRF);
plot(real(P(:)), imag(P(:)), 'bx');
hold on
plot([0 0],[-10 10],'k:', [-10 10],[0 0],'k:');
latexylabel('Im');
rotateY;
latexxlabel('Re');
latextitle('Netz mit Zustandsr\"uckf\"uhrung');
axis([-0.5 0.1 -0.3 0.3]);
axis('square');
xticks([-0.5 -0.25 0]);
yticks([-0.25 0 0.25]);
hold off


K1=K(1, :);
K2=K(2, :);
K3=K(3, :);
CS1=[C1          0 0 0
     zeros(1, n) 1 0 0];
CS2=[C2          0 0 0
     zeros(1, n) 0 1 0];
CS3=[C3          0 0 0
     zeros(1, n) 0 0 1];
BS1=BI(:, 1);
BS2=BI(:, 2);
BS3=BI(:, 3);
[V, Lambda]=eig(AI-BI*K);
W=eye(n+m, n+m);
CVW1=CS1*V*W;
CVWplus1=CVW1'*inv(CVW1*CVW1');
Ky1=real(K1*V*W*CVWplus1);  % real(...) um komplexe Darstellung zu verhindern
CVW2=CS2*V*W;
CVWplus2=CVW2'*inv(CVW2*CVW2');
Ky2=real(K2*V*W*CVWplus2);
CVW3=CS3*V*W;
CVWplus3=CVW3'*inv(CVW3*CVW3');
Ky3=real(K3*V*W*CVWplus3);
Ky=[Ky1(1,1) 0        0       Ky1(1,2)   0        0
    0       Ky2(1,1)  0        0        Ky2(1,2)  0
    0        0       Ky3(1,1)  0         0        Ky3(1,2)];
%
Ag=AI-BI*Ky*CI;
Bg=EI-BI*Ky*FI;
Cg=[C(1, 1:10) zeros(1,3)];
Dg=F(1, 1:m);
%
figure(3)
NetzARF=ss(Ag, Bg, Cg, Dg);
P=pzmap(NetzARF);
plot(real(P(:)), imag(P(:)), 'bx');
hold on
plot([0 0],[-10 10],'k:', [-10 10],[0 0],'k:');
latexylabel('Im');
rotateY;
latexxlabel('Re');
latextitle('Netz mit dezentraler Regelung');
axis([-0.5 0.1 -0.3 0.3]);
axis('square');
xticks([-0.5 -0.25 0]);
yticks([-0.25 0 0.25]);
hold off


%%    Stör\"ubergangsfunktion des dezentral geregelten Netzes
%
figure(4);
Regelkreis1=ss(Ag, Bg, Cg, Dg);
Regelkreis2=ss(Ag, Bg, -Ky*CI, -Ky*FI);
[F1, T]=step(Regelkreis1);
U1=step(Regelkreis2, T);
subplot(6,1,[1,2]);
plot(T, 100*F1(:, 1, 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*U1(:,1),'b', [0 25], [100 100], 'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$u_1$$ in MW');
rotateY;
axis([0 25 0 170]);
xticks([0 10 20]);
yticks([0 50 100 150]);


%%
%
%    Verbesserte Regelung
%
%
figure(5);
W(1,1)=0;
W(13,13)=0;
W(12,12)=0;
W(2,2)=80;
W(3,3)=80;
W(4,4)=50;
W(5,5)=50;
W(6,6)=50;
W(7,7)=50;
CVW1=CS1*V*W;
CVWplus1=CVW1'*inv(CVW1*CVW1');
Ky1=real(K1*V*W*CVWplus1);
CVW2=CS2*V*W;
CVWplus2=CVW2'*inv(CVW2*CVW2');
Ky2=real(K2*V*W*CVWplus2);
CVW3=CS3*V*W;
CVWplus3=CVW3'*inv(CVW3*CVW3');
Ky3=real(K3*V*W*CVWplus3);
Ky=[Ky1(1,1) 0        0       Ky1(1,2)   0        0
    0       Ky2(1,1)  0        0        Ky2(1,2)  0
    0        0       Ky3(1,1)  0         0        Ky3(1,2)];
%
Ag=AI-BI*Ky*CI;
Bg=EI-BI*Ky*FI;
Cg=[C(1, 1:10) zeros(1,3)];
Dg=F(1, 1:m);
NetzARF2=ss(Ag, Bg, Cg, Dg);
P=pzmap(NetzARF2);
plot(real(P(:)), imag(P(:)), 'bx');
hold on
plot([0 0],[-10 10],'k:', [-10 10],[0 0],'k:');
latexylabel('Im');
rotateY;
latexxlabel('Re');
latextitle('Netz mit verbesserter Regelung');
axis([-0.5 0.1 -0.3 0.3]);
axis('square');
xticks([-0.5 -0.25 0]);
yticks([-0.25 0 0.25]);
hold off


%%  Stör\"ubergangsfunktion bei verbesserter Regelung
%
%
figure(6);
Regelkreis1=ss(Ag, Bg, Cg, Dg);
[F1, T]=step(Regelkreis1);
Regelkreis2=ss(Ag, Bg, -Ky*CI, -Ky*FI);
U1=step(Regelkreis2, T);
subplot(6,1,[1,2]);
plot(T, 100*F1(:,1,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 verbesserter dezentraler Regelung');
subplot(6,1,[3,4]);
plot(T, 100*U1(:,1),'b', [0 25], [100 100], 'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$u_1$$ in MW');
rotateY;
axis([0 25 0 170]);
xticks([0 10 20]);
yticks([0 50 100 150]);

%




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