%   Script DNAUQ.M zur Herstellung von Bildern
%
%    Dezentrale UQ-Regelung (DNA-Entwurf)
%
% verwendet BSPUQREG.M, EPSBILD.M
% Funktionen: lambdaP, gerkreis
%
% J. Lunze
% 6.1.1997
% für 2. Auflage: 16.12.2001
% für 3. Auflage: 19.4.2004
% für 4. Auflage: 11.1.2006
% für 5. Auflage: 19.10.2007
% für 7. Auflage 31.3.2012
% für 8. Auflage: 4.6.2014
% für 10. Auflage: 9.9.2019
%
%
clear
close all
Bspuqreg;    % Modell (A, B, C, D)

%%
%
%   Übergangsfunktionsmatrix der Regelstrecke
%
%
figure(1);
Regelstrecke=ss(A, B, C, D);
T=0:0.03:4;
Y=step(Regelstrecke, T);
subplot(3,3,1);
plot(T, Y(:,1, 1), 'b');
xleer;
latexylabel('$$h_{11}$$');
rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
latextitle('\"Ubergangsfunktionsmatrix');
subplot(3,3,2);
plot(T, Y(:, 1,2), 'b');
xleer;
yleer;
%ylabel('{\it h}_{12}');
%rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
subplot(3,3,3);
plot(T, Y(:, 1,3), 'b');
xleer;
yleer;
%ylabel('{\it h}_{13}');
%rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
latextitle('UQ-Regelstrecke');
subplot(3,3,4);
plot(T, Y(:, 2,1), 'b');
xleer;
latexylabel('$$h_{21}$$');
rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
subplot(3,3,5);
plot(T, Y(:, 2,2), 'b');
xleer;
%ylabel('{\it h}_{22}');
%rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
subplot(3,3,6);
plot(T, Y(:, 2,3), 'b');
xleer;
%ylabel('{\it h}_{23}');
%rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
subplot(3,3,7);
plot(T, Y(:, 3,1), 'b');
latexxlabel('$$t$$ in s');
latexylabel('$$h_{31}$$');
rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
subplot(3,3,8);
plot(T, Y(:, 3,2), 'b');
latexxlabel('$$t$$ in s');
%ylabel('{\it h}_{32}');
%rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
subplot(3,3,9);
plot(T, Y(:, 3,3), 'b');
latexxlabel('$$t$$ in s');
%ylabel('{\it h}_{33}');
%rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);



%%   Regler für das Netz 1 wie in Beispiel RT1-11.1
%
%
figure(2);
subplot(4,1,[1,2]);
%    Ausblenden des Teilsystems 1 (A, Bs1, Cs1, Ds1)
%    aus Gesamtmodell  (vgl. BSPUQREG.M)
Bs1=B(1:n, 1);
Cs1=C(1, 1:n);
Ds1=0;
%    Regler (s^2 + 3s + 2)/(s^2 +0,22s)
zr=[1 3 2];
nr=[0.22 1 0];
[Ar, Br, Cr, Dr]=tf2ss(zr, nr);
%           nur TR1 an Regelstrecke: (Ag1, Bg1, Cg1, Dg1)
Ag1=[A-Bs1*Dr*Cs1  Bs1*Cr
     -Br*Cs1        Ar   ];
Bg1=[Bs1*Dr
     Br];
Cg1=[Cs1 zeros(1,2)];
Dg1=0;
geregeltesNetz1=ss(Ag1, Bg1, Cg1, Dg1);
Y=step(geregeltesNetz1, T);
plot(T, Y, 'b');
latexxlabel('$$t$$ in s');
yticks([0 1]);
latexylabel('$$h_{\rm w11}$$');
rotateY;
latextitle('F\"uhrungs\"ubergangsfunktion des ersten UQ-Regelkreises');
axis([0 4 0 1.2]);
xticks([0 1 2 3 4 5]);
yticks([0 1]);


% dieselben Regler für alle drei Kraftwerke

%%
%
%   Ortskurve der R\"uckf\"uhrdifferenzmatrix
%
%
figure(3);
%             Drei dezentrale Regler: (Ardez, Brdez, Crdez, Drdez)
nri=max(size(Ar));
Ardez=[Ar zeros(nri, 2*nri)
       zeros(nri, nri) Ar zeros(nri, nri)
       zeros(nri, 2*nri) Ar];
Brdez=[Br zeros(nri, 2)
       zeros(nri, 1) Br zeros(nri, 1)
       zeros(nri, 2) Br];
Crdez=[Cr zeros(1, 2*nri)
       zeros(1, nri) Cr zeros(1, nri)
       zeros(1, 2*nri) Cr];
Drdez=diag([Dr Dr Dr]);
%              Offene Kette mit dezentralem Regler: (A0, B0, C0 D0)
offeneKette1=ss(A, B, C, D);
dezRegler=ss(Ardez, Brdez, Crdez, Drdez);
[A0, B0, C0, D0]=ssdata(series(offeneKette1, dezRegler));
%             Berechnung der R\"uckf\"uhrdifferenzmatrix FF
W=logspace(-3, 2, 150);
Wsize=max(size(W));
Rueckfuehrdifferenz=ss(A0, B0, C0, D0+eye(3,3));
[Betrag, Phase, W]=bode(Rueckfuehrdifferenz, 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;
subplot(1,2,1)
plot(real(detF), imag(detF), 'b', real(detF), -imag(detF),'b--', [0 0], [-100 100], 'k:', [-100 100], [0 0], 'k:');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Ortskurve von $$det(F(j\omega))$$');
axis([-65 15 -35 35]);
xticks([-40 0 40]);
yticks([-40 0 40]);
axis('square');
%
subplot(1,2,2)
plot(real(detF), imag(detF), 'b', real(detF), -imag(detF),'b--', [0 0], [-5 2], 'k:', [-4 10], [0 0], 'k:');
latexxlabel('Re');
%latexylabel('Im');
yleer;
latextitle('Ortskurve von $$det(F(j\omega))$$');
axis([-4 1.5 -4.5 1]);
xticks([-4 -2 0]);
yticks([-4 -2 0 2]);
axis('square');



%%    Überprüfung der Diagonaldominanz der R\"uckf\"uhrdifferenzmatrix
%
%
figure(4);
subplot(4,1,[1,2]);
semilogx(W, lP, 'b');
latexxlabel('$$\omega$$ in $$\frac{\rm rad}{\rm s}$$');
yticks([0 1]);
latexylabel('$$\lambda_{\rm P}$$');
rotateY;
latextitle('Perronwurzel f\"ur Pr\"ufung der Diagonaldominanz');
axis([0.001 100 0 1]);
xticks([0.0001 0.001 0.01 0.1 1 10 100]);
yticks([0 1]);
grid('on');



%%
%
%     Gershgorin-Bänder für die drei Regelkreise
%
%
figure(5);
%     Offene Kette mit TR1: G0(jw)=G(1,1,:)-1
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], 'b');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Gershgorinkreise f\"u UQ-Regelung (TR 1)');
axis([-1.25 1.25 -2 0.5]);
xticks([-1 0 1]);
yticks([-2 -1 0]);
axis('square');
hold off;




figure(6);
%     Offene Kette mit TR2: G0(jw)=G(2,2,:)-1
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], 'b');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Gershgorinkreise f\"ur UQ-Regelung (TR 2)');
axis([-1.25 1.25 -2 0.5]);
xticks([-1 0 1]);
yticks([-2 -1 0]);
axis('square');
hold off;


figure(7);
%     Offene Kette mit TR3: G0(jw)=G(3,3,:)-1
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;
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 UQ-Regelung (TR 3)');
axis([-1.25 1.25 -2 0.5]);
xticks([-1 0 1]);
yticks([-2 -1 0]);
axis('square');
hold off;



%%
%
%     Verbesserte Gershgorin-Bänder für den dritten Regelkreis
%
%
figure(8);
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;
for k=1:length(Betrag(3,3,:))
    R3v(1,1,k)=Betrag(3,3,k).*lP(k)';
end
for k=1:3:Wsize
   gerkreis(G(3,3,k)-1, R3v(1,1,k))
end;
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Verbessertes Gershgorinband (UQ-TR 3)');
plot([-1 -1 ],[-0.05 0.05], 'b');
axis([-1.25 1.25 -2 0.5]);
xticks([-1 0 1]);
yticks([-2 -1 0]);
axis('square');
hold off;


%%
%
%    Vergleich beider Abschätzungen
%
%
figure(9);
subplot(4,1,[1,2]);
%    Radien der Gershgorinkreise bezogen auf das Hauptdiagonalelement
V1(1,1,:)=R1./Betrag(1,1,:);
V2(1,1,:)=R2./Betrag(2,2,:);
V3(1,1,:)=R3./Betrag(3,3,:);
semilogx(W, V1(:), 'b--', W, V2(:), 'b--', W, V3(:), 'b--', W, lP(:), 'b');
axis([0.001, 100, 0, 1]);%xticks([]);
%yticks([]);

latexxlabel('$$\omega$$ in $$\frac{\rm rad}{\rm s}$$');
yticks([0 1]);
xticks([0.0001 0.001 0.01 0.1 1 10 100]);
grid('on');
rotateY;
latextitle('Relative Radien der Gershgorinkreise');
latextext(0.005, 0.3, '$$TR_2$$');
%text(0.005, 0.63, 'TR_1');
latextext(0.005, 0.8, '$$TR_3$$');
%text(0.005, 0.51, '\lambda_P');


%%
%
%    Verhalten des dezentral geregelten Systems
%
%
figure(10)
%    Geschlossener Regelkreis: (Ag, Bg, Cg, Dg)
Ag=[A-B*Drdez*C  B*Crdez
    -Brdez*C     Ardez ];
Bg=[B*Drdez
    Brdez];
Cg=[C zeros(3, 3*nri)];
Dg=zeros(3,3);
Regelkreis=ss(Ag, Bg, Cg, Dg);
Y=step(Regelkreis, T);
subplot(3,3,1);
plot(T, Y(:,1, 1), 'b');
xleer;
latexylabel('$$h_{11}$$');
rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
latextitle('F\"uhrungs\"ubergangsfunktionsmatrix');
subplot(3,3,2);
plot(T, Y(:,1,2), 'b');
xleer;
%ylabel('{\it h}_{12}');
%rotateY;
yleer
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
subplot(3,3,3);
plot(T, Y(:, 1,3), 'b');
xleer;
%ylabel('{\it h}_{13}');
%rotateY;
latextitle('dezentrale UQ-Regelung');
axis([0, 4, 0, 1]);
%xticks([]);
yleer;
subplot(3,3,4);
plot(T, Y(:, 2,1), 'b');
xleer;
latexylabel('$$h_{21}$$');
rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
subplot(3,3,5);
plot(T, Y(:, 2,2), 'b');
xleer;
%ylabel('{\it h}_{22}');
%rotateY;
yleer;
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
subplot(3,3,6);
plot(T, Y(:, 2,3), 'b');
xleer;
%ylabel('{\it h}_{23}');
%rotateY;
yleer;
axis([0, 4, 0, 1]);
%xticks([]);
%yticks([]);
subplot(3,3,7);
plot(T, Y(:, 3,1), 'b');
latexxlabel('$$t$$ in s');
latexylabel('$$h_{31}$$');
rotateY;
yticks([0 1]);
axis([0, 4, 0, 1]);
xticks([0 2 4]);
%
subplot(3,3,8);
plot(T, Y(:, 3,2), 'b');
latexxlabel('$$t$$ in s');
%ylabel('{\it h}_{32}');
%rotateY;
yleer;
axis([0, 4, 0, 1]);
xticks([0 2 4]);
%
subplot(3,3,9);
plot(T, Y(:, 3,3), 'b');
latexxlabel('$$t$$ in s');
%ylabel('{\it h}_{33}');
%rotateY;
yleer;
axis([0, 4, 0, 1]);
xticks([0 2 4]);





%%
%
%    Abschätzung des Amplitudenganges des ersten Regelkreises
%
%
figure(11);
subplot(4,1,[1,2]);
%   Strecke mit TR1
StreckevonTR1=ss(Ag1, Bg1, Cg1, Dg1);
[BetragGw1, PhaseGw1]=bode(StreckevonTR1, W);
%semilogx(W, 20*log10(BetragGw1(:)));
%hold on;
%    Kopplungen in Regelstrecke
Regelstrecke=ss(A, B, C, D);
[Betrag, Phase]=bode(Regelstrecke,W);
%    Berechnung von Guw
%    Regelkreise bez\"uglich Eingang f und Ausgang u
Bf1=-Bg1;
Cf1=[-Dr*Cs1 Cr];
Df1=-Dr;
Regelkreis1=ss(Ag1, Bf1, Cf1, Df1);
[Betraguw1, Phaseuw1, W]=bode(Regelkreis1,  W);
%    Ausblenden des Teilsystems 2 (A, Bs2, Cs2, Ds2)
Bs2=B(1:n, 2);
Cs2=C(2, 1:n);
Ds2=0;
%    Regelkreis bzgl. Eingang f und Ausgang u
Ag2=[A-Bs2*Dr*Cs2  Bs2*Cr
     -Br*Cs2        Ar   ];
Bf2=[-Bs2*Dr
     -Br];
Cf2=[-Dr*Cs2 Cr];
Df2=-Dr;
Regelkreis2=ss(Ag2, Bf2, Cf2, Df2);
[Betraguw2, Phaseuw2, W]=bode(Regelkreis2,  W);
%    Ausblenden des Teilsystems 3 (A, Bs3, Cs3, Ds3)
Bs3=B(1:n, 3);
Cs3=C(3, 1:n);
Ds3=0;
%    Regelkreis bzgl. Eingang f und Ausgang u
Ag3=[A-Bs3*Dr*Cs3  Bs3*Cr
     -Br*Cs3        Ar   ];
Bf3=[-Bs3*Dr
     -Br];
Cf3=[-Dr*Cs3 Cr];
Df3=-Dr;
Regelkreis3=ss(Ag3, Bf3, Cf3, Df3);
[Betraguw3, Phaseuw3, W]=bode(Regelkreis3,  W);
%   Berechnung des Fehlers
for k=1:Wsize
   V2 = [   0             Betrag(2,1,k)  Betrag(1,3,k);
         Betrag(2,1,k)      0             Betrag(3,2,k);
         Betrag(3,1,k)    Betrag(3,2,k)         0      ];
   Gyf = abs(diag([1/G(1,1,k), 1/G(2,2,k), 1/G(3,3,k)]));
   Guw = diag([Betraguw1(1,1,k), Betraguw2(1,1,k), Betraguw3(1,1,k)]);
   Guf = Guw;
   Vyw11(k,1) = [1 0 0]*Gyf*V2*inv(eye(3,3) - Guf*V2)*Guw*[1; 0; 0];
   Vyw12(k,1) = [1 0 0]*Gyf*V2*inv(eye(3,3) - Guf*V2)*Guw*[0; 1; 0];
end;
for k=1:length(BetragGw1)
    Gdrucko(k)=BetragGw1(k)+Vyw11(k);
    Gdrucku(k)=BetragGw1(k)-Vyw11(k);
end
semilogx(W, 20*log10(Gdrucko(:)), 'b', W, 20*log10(Gdrucku(:)), 'b');
latexxlabel('$$\omega$$ in $$\frac{\rm rad}{\rm s}$$');
yticks([-60 0]);
latexylabel('$$|G_{\rm w11}|$$ in dB');
rotateY;
latextitle('Amplitudengang des ersten Regelkreises');
axis([0.001 100 -60 20]);
xticks([0.0001 0.001 0.01 0.1 1 10 100]);
grid('on');
hold off;


%%
%
%     Fehlerschranken V11 und V12
%
%
figure(12);
subplot(4,1,[1,2]);
semilogx(W, 20*log10(Vyw11(:)), 'b--', W, 20*log10(Vyw12(:)), 'b');
latexxlabel('$$\omega$$ in $$\frac{\rm rad}{\rm s}$$');
latexylabel('$$|V|$$ in dB');
rotateY;
latextitle('Amplitudengang der Fehlerschranken $$V_{11}$$ und $$V_{12}$$');
latextext(10, -90, '$$V_{12}$$');
latextext(10, -20, '$$V_{11}$$');
axis([0.001 100 -100 20]);
grid('on');
xticks([0.0001 0.001 0.01 0.1 1 10 100]);
yticks([-100 -50 0]);

%%
%
%    Führungsübergangsfunktion des ersten Regelkreises
%
%
figure(13);
subplot(4,1,[1,2]);
Regelkreis1=ss(Ag1, Bf1, Cf1, Df1);
Y10=step(Regelkreis1, T);
plot(T, Y(:,1,1), 'b', T, Y10(:,1), 'b--');
latexxlabel('$$t$$ in s');
yticks([0 1]);
latexylabel('$$h_{\rm w11}$$');
rotateY;
latextitle('F\"uhrungs\"ubergangsfunktion bei TR-1 - - und dezentralem Regler --');
axis([0 4 0 1.2]);
xticks([0 1 2 3 4 5]);


%%
%
%    Wiederholung der Abschätzung des Regelkreises bei
%    Vergrößerung der Reglerverstärkungen
%
%
figure(14)
subplot(4,1,[1,2]);
kn=6;
Arn=Ar;
Brn=kn*Br;
Crn=Cr;
Drn=kn*Dr;
%           nur TR1 an Regelstrecke: (Ag1, Bg1, Cg1, Dg1)
Ag1n=[A-Bs1*Drn*Cs1  Bs1*Crn
     -Brn*Cs1        Arn   ];
Bg1n=[Bs1*Drn
     Brn];
Cg1n=[Cs1 zeros(1,2)];
Dg1n=0;
%   Strecke mit TR1
RegelTR1=ss(Ag1n, Bg1n, Cg1n, Dg1n);
[BetragGw1n, PhaseGw1n, W]=bode(RegelTR1,W);
%semilogx(W, 20*log10(BetragGw1n(:)));
%hold on;
%    Berechnung von Guw
%    Regelkreise bez\"uglich Eingang f und Ausgang u
Bf1n=-Bg1n;
Cf1n=[-Drn*Cs1 Cr];
Df1n=-Drn;
RegelkreisGuw11=ss(Ag1n, Bf1n, Cf1n, Df1n);
[Betraguw1n, Phaseuw1n]=bode(RegelkreisGuw11, W);
RegelkreisGuw12=ss(Ag1n, Bf1n, [Cs1 0 0], 1);
[Betragyf1n, Phaseyf1n]=bode(RegelkreisGuw12, W);
Ag2n=[A-Bs2*Drn*Cs2  Bs2*Crn
     -Brn*Cs2        Arn   ];
Bf2n=[-Bs2*Drn
     -Brn];
Cf2n=[-Drn*Cs2 Crn];
Df2n=-Drn;
RegelkreisGuw21=ss(Ag2n, Bf2n, Cf2n, Df2n);
[Betraguw2n, Phaseuw2n]=bode(RegelkreisGuw21, W);
RegelkreisGuw22=ss(Ag2n, Bf2n, [Cs2 0 0], 1);
[Betragyf2n, Phaseyf2n]=bode(RegelkreisGuw22,  W);
Ag3n=[A-Bs3*Drn*Cs3  Bs3*Crn
     -Brn*Cs3        Arn   ];
Bf3n=[-Bs3*Drn
     -Brn];
Cf3n=[-Drn*Cs3 Crn];
Df3n=-Drn;
RegelkreisGuw31=ss(Ag3n, Bf3n, Cf3n, Df3n);
[Betraguw3n, Phaseuw3n]=bode(RegelkreisGuw31,  W);
RegelkreisGuw32=ss(Ag3n, Bf3n, [Cs3 0 0], 1);
[Betragyf3n, Phaseyf3n]=bode(RegelkreisGuw32,  W);
%   Berechnung des Fehlers
for k=1:Wsize
   V2 = [   0             Betrag(2,1,k)  Betrag(1,3,k);
         Betrag(2,1,k)      0             Betrag(3,2,k);
         Betrag(3,1,k)    Betrag(3,2,k)         0      ];
   Guwn = diag([Betraguw1n(1,1,k), Betraguw2n(1,1,k), Betraguw3n(1,1,k)]);
   Gyfn = diag([Betragyf1n(1,1,k), Betragyf2n(1,1,k), Betragyf3n(1,1,k)]);
   Gufn = Guwn;
   Pmatrixn = Gufn*V2;
   lPn(k,1) = lambdaP(Pmatrixn);
   Vn(k,1) = max([Pmatrixn(1,2)+Pmatrixn(1,3), Pmatrixn(2,1)+Pmatrixn(3,1)]);
   Vn(k,2) = max([Pmatrixn(2,1)+Pmatrixn(2,3), Pmatrixn(1,2)+Pmatrixn(3,2)]);
   Vn(k,3) = max([Pmatrixn(3,1)+Pmatrixn(3,2), Pmatrixn(1,3)+Pmatrixn(2,3)]);
   Vyw11n(k,1) = [1 0 0]*Gyfn*V2*inv(eye(3,3) - Gufn*V2)*Guwn*[1; 0; 0];
   Vyw12n(k,1) = [1 0 0]*Gyfn*V2*inv(eye(3,3) - Gufn*V2)*Guwn*[0; 1; 0];
end;
for k=1:length(BetragGw1n)
    Gdrucko(k)=BetragGw1n(k)+Vyw11n(k);
    Gdrucku(k)=BetragGw1n(k)-Vyw11n(k);
end
semilogx(W, 20*log10(Gdrucko(:)), 'b', W, 20*log10(Gdrucku(:)), 'b');
latexxlabel('$$\omega$$ in $$\frac{\rm rad}{\rm s}$$');
yticks([-60 0 ]);
latexylabel('$$|G_{\rm w11}|$$ in dB');
rotateY;
latextitle('Erster Regelkreis nach Verst\"arkungserh\"ohung');
grid('on')
axis([0.001 100 -60 20]);
xticks([0.0001 0.001 0.01 0.1 1 10 100]);
hold off;


%%
%
%     Perronwurzel nach Verstärkungserhöhung
%
%
figure(15);
subplot(4,1,[1,2]);
semilogx(W, lPn, 'b');
latexxlabel('$$\omega$$ in $$\frac{\rm rad}{\rm s}$$');
latexylabel('$$\lambda_{\rm P}$$');
rotateY;
latextitle('Perronwurzel nach Verst\"arkungserh\"ohung');
axis([0.001 100 0 1]);
xticks([0.0001 0.001 0.01 0.1 1 10 100]);
yticks([0 1]);
grid('on');

%%
%
%     neue Fehlerschranken V11 und V12
%
%
figure(16);
subplot(4,1,[1,2]);
semilogx(W, 20*log10(Vyw11n), 'b--', W, 20*log10(Vyw12n), 'b');
latexxlabel('$$\omega$$ in $$\frac{\rm rad}{\rm s}$$');
latexylabel('$$|V|$$ in dB');
rotateY;
latextitle('Amplitudengang der neuen Fehlerschranken');
latextext(20, -80, '$${V}_{12}$$');
latextext(30, -20, '$${V}_{11}$$');
axis([0.001 100 -100 20]);
xticks([0.0001 0.001 0.01 0.1 1 10 100]);
yticks([-100 -50 0]);
grid('on');

%%
%
%    F\"uhrungs\"ubergangsfunktion des ersten Regelkreises
%    nach Verstärkungserhöhung
%
%
figure(17);
T=0:0.03:4;
subplot(4,1,[1,2]);
Bg1n = -Bf1n;
Cg1n = Cg1;
Dg1n = Dg1;
RegelkreisNetz1=ss(Ag1, Bg1, Cg1, Dg1);
Y10=step(RegelkreisNetz1, T);
%     neuer dezentraler Regler
Ardezn=[Arn zeros(nri, 2*nri)
       zeros(nri, nri) Arn zeros(nri, nri)
       zeros(nri, 2*nri) Arn];
Brdezn=[Brn zeros(nri, 2)
       zeros(nri, 1) Brn zeros(nri, 1)
       zeros(nri, 2) Brn];
Crdezn=[Crn zeros(1, 2*nri)
       zeros(1, nri) Crn zeros(1, nri)
       zeros(1, 2*nri) Crn];
Drdezn=diag([Drn Drn Drn]);
%    geschlossener Regelkreis
Agn=[A-B*Drdezn*C  B*Crdezn
    -Brdezn*C     Ardezn ];
Bgn=[B*Drdezn
    Brdezn];
Cgn=Cg;
Dgn=zeros(3,3);
ErsterRegelkreis=ss(Agn, Bgn, Cgn, Dgn);
Y1=step(ErsterRegelkreis,  T);
plot(T, Y1(:,1,1), 'b', T, Y10(:,1,1), 'b--');
latexxlabel('$$t$$ in s');
yticks([0 1]);
latexylabel('$$h_{\rm w11}$$');
rotateY;
latextitle('F\"uhrungs\"ubergangsfunktion bei TR 1 - - und dezentralem Regler --');
axis([0 4 -0.2 1.4]);
xticks([0 1 2 3 4]);


%%
%
%   Vergleich der Abschätzungen
%
%
figure(18)
subplot(4,1,[1,2]);
semilogx(W, Vn, 'b--', W, lPn, 'b');
axis([0.001, 100, 0, 1]);
latexxlabel('$$\omega$$ in $$\frac{\rm rad}{\rm s}$$');
yticks([0 1]);
xticks([0.0001 0.001 0.01 0.1 1 10 100]);
latexylabel('$$\lambda_{\rm P}$$');
rotateY;
latextitle('Relative Radien der Gershgorinkreise');
%text(0.005, 0.4, 'TR_2');
%text(0.005, 0.63, 'TR_1');
%text(0.005, 0.75, 'TR_3');
%text(0.005, 0.45, '\lambda_P');
grid('on');



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