%   Script ZRFREAK.M zur Herstellung von Bildern
%
%   Zustandsr\"uckf\"uhrung für Reaktoren
%
% verwendet BSPREA3.M, EPSBILD.M
% aus RT1: RLOCPRINT
%
% J. Lunze
% 12.5.1999              Daten wie 13.9.1996
% für 2. Auflage: 16.12.2001
% für 3. Auflage: 19.4.2004
% für 5. Auflage: 18.10.2007
% für 6. Auflage: 5.1.2010
% für 7. Auflage: 31.3.2012
% für 8. Auflage: 3.4.3014
% für 10. Auflage: 9.9.2019
%
%
clear
close all
%
%    Zustandsrückführung für zwei Reaktoren
%
Bsprea3;

%%  Messgrößen
figure(1);
sigma = [-1 -2];
K = place(A, B, sigma);
Ag=A-B*K;
Bg=B;
Cg=C;
Dg=D;
x0=[1 0];
Reaktor=ss(Ag, Bg, Cg, Dg);
[Y1, T]=initial(Reaktor, x0);
Reaktor2=ss(Ag, Bg, -K, zeros(m,m));
U1=initial(Reaktor2, x0, T);
sigma = [-1.5 -2.5];
K = place(A, B, sigma);
Ag=A-B*K;
Bg=B;
Cg=C;
Dg=D;
x0=[1 0];
Reaktor=ss(Ag, Bg, Cg, Dg);
Y2=initial(Reaktor, x0, T);
Reaktor2=ss(Ag, Bg, -K, zeros(m,m));
[U2, T, X]=initial(Reaktor2, x0, T);
subplot(4,1,[1,2]);
plot(T, Y1,'b', [-1 0 0],[0 0 Y1(1,1)]);
hold on
plot(T, Y2, 'b--', [-1 0 0],[0 0 Y2(1,1)],'b--', [-1 5], [0 0],'k:');
xleer;
latexylabel('$$y$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
yticks([0 0.5]);
axis([-0.5 5 -0.1 0.7]);
yticks([0 0.5]);
latextitle('Zwei Reaktoren mit Zustandsr\"uckf\"uhrung');
subplot(4,1,[3,4]);
plot(T, U1,'b', [-1 0 0],[0 0 U1(1,1)],'b', T, U2,...
            'b--', [-1 0 0],[0 0 U2(1,1)], 'b--', [0, 5], [0, 0], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
yticks([-5 0]);
axis([-0.5 5 -5.5 0.5]);
xticks([0 1 2 3 4]);
yticks([-5 0 5]);


%%
%
%   Wurzelortskurve für Ausgangsrückführung
%
%
figure(2)
Reaktor=ss(A, B, C, D);
rlocprintblue(Reaktor);
latextitle('Wurzelortskurve der Reaktorregelung');
hold on;
Reaktorregelung=ss(A, -B, C, -D);
rlocprintblue(Reaktorregelung);
axis([-4 4 -4 4]);
axis('square');
xticks([-4 -2 0 2 4]);
yticks([-4 -2 0 2 4]);
hold off;
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextext(2, 0.5, '$$k<0 \; \rightarrow$$');
latextext(-3.5, 0.5, '$$\leftarrow \; k<0$$');
latextext(-3, 3, '$$k>0 \; \uparrow$$');
latextext(-3, -3, '$$k>0 \; \downarrow$$');


%%
%
%    Näherung der ZRF durch ARF
%
echo  off
figure(3);
Regelstrecke=ss(A, B, C, D);
P=pzmap(Regelstrecke);
for k=1:length(P)
    plot(real(P(k)), imag(P(k)), 'bx');
    hold on;
end
sigma=[-1 -2];
K=place(A, B, sigma);
[V, Lambda]=eig(A-B*K);
for w=[0.5, 1, 2, 3, 5, 10]
   W=diag([1 w]);
   CVW=C*V*W;
   CVWplus=CVW'*inv(CVW*CVW');
   Ky=K*V*W*CVWplus;
   Naeherung=ss(A-B*Ky*C, B, C, D);
   P=pzmap(Naeherung);
   for k=1:length(P)
       plot(real(P(k)), imag(P(k)), 'b*');
   end
end;
w=100;
W=diag([1 w]);
CVW=C*V*W;
CVWplus=CVW'*inv(CVW*CVW');
Ky=K*V*W*CVWplus;
Naeherung=ss(A-B*Ky*C, B, C, D);
P=pzmap(Naeherung);
for k=1:length(P)
    plot(real(P(k)), imag(P(k)), 'bd');
end
plot([-3 1], [0 0], 'k:', [0 0], [-5 5], 'k:');
latexylabel('Im');
rotateY;
latexxlabel('Re');
latextitle('Eigenwerte der Reaktoren mit Ausgangsr\"uckf\"uhrung');
axis([-2.1 0.1 -1.05 1.05]);
axis('square')
xticks([-2 -1 0]);
yticks([-1 0 1]);
hold off;
latextext(-1.8, -0.2, '$$w_{22} \; \rightarrow$$');
latextext(-0.9, -0.2, '$$\leftarrow \; w_{22}$$');


%%
%
%               Eigenbewegung mit Ausgangsr\"uckf\"uhrung
%
%
figure(4);
Ag=A-B*K;
Reaktor=ss(Ag, Bg, Cg, Dg);
[Y1, T, X]=initial(Reaktor, x0);
Reaktor2=ss(Ag, Bg, -K, zeros(m,m));
[U1, T, X]=initial(Reaktor2, x0, T);
Ag=A-B*Ky*C;
Reaktor=ss(Ag, Bg, Cg, Dg);
[Y2, T]=initial(Reaktor, x0, T);
Reaktor3=ss(Ag, Bg, -Ky*C, zeros(m,m));
U2=initial(Reaktor3, x0, T);
subplot(4,1,[1,2]);
plot(T, Y2, 'b', [-1 0 0],[0 0 Y2(1,1)]);
hold on
plot(T, Y1, 'b--', [-1 0 0],[0 0 Y1(1,1)],'b--', [-1 5], [0 0],'k:');
xleer;
latexylabel('$$y$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
yticks([0 0.5]);
axis([-0.5 5 -0.1 0.7]);
latextitle('Zwei Reaktoren mit Ausgangsr\"uckf\"uhrung');
subplot(4,1,[3,4]);
plot(T, U2,'b', [-1 0 0],[0 0 U2(1,1)], T, U1, 'b--', [-1 0 0],[0 0 U1(1,1)], 'b--', [0, 5], [0, 0], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
yticks([-2 0]);
axis([-0.5 5 -2.5 0.5]);


%%
%
%    Näherung der ZRF durch ARF
%
echo  off
figure(5);
sigma=[-1.5 -2.5];
K=place(A, B, sigma);
[V, Lambda]=eig(A-B*K);
for w=logspace(-1, 2, 20)
   W=diag([w 1]);
   CVW=C*V*W;
   CVWplus=CVW'*inv(CVW*CVW');
   Ky=K*V*W*CVWplus;
   Naeherung2=ss(A-B*Ky*C, B, C, D);
   P=pzmap(Naeherung2);
   for k=1:length(P)
       plot(real(P(k)), imag(P(k)), 'bx');
       hold on;
   end
end;
plot([0 0],[-5 5],'k:', [-5 5],[0 0],'k:');
latexylabel('Im');
rotateY;
latexxlabel('Re');
latextitle('Eigenwerte der Reaktoren mit Ausgangsr\"uckf\"uhrung');
hold off;
latextext(-2.3, -0.2, '$$w_{22}\; \rightarrow$$');
latextext(-0.7, -0.2, '$$\leftarrow \; w_{22}$$');
axis([-2.7 0.3 -1.5 1.5]);
axis('square');
xticks([-2 -1 0]);
yticks([-1 0 1]);


%%
%
%  Stabilitätspr\"ufung mit dem Nyquistkriterium
%
%
figure(6)
%  Reaktor (A, B, I, 0) mit Zustandsr\"uckf\"uhrung K
%  R\"uckf\"uhrdifferenzmatrix
A0 = A;
B0 = B*K;
C0 = eye(n,n);
D0 = eye(n,n);
W=logspace(-3, 2, 150);
Wsize=max(size(W));
Rueckfuehrdifferenzmatrix=ss(A0, B0, C0, D0);
[Real, Imag]=nyquist(Rueckfuehrdifferenzmatrix, W);
for k=1:Wsize
   detF(k)=det(Real(:,:,k)+i*Imag(:,:,k));
end;
plot(detF, 'b');
hold on;
plot(detF','b--');
plot([0 0],[-10 10],'k:',[-10 10],[0 0],'k:');
plot(det(dcgain(Rueckfuehrdifferenzmatrix)), 0, 'o');
plot(1, 0, 'b*');
axis([-1 6 -3.5 3.5]);
axis('square');
xticks([0 2 4 6]);
yticks([-2 0 2]);
text(3, -2.5, '\leftarrow');
text(3, 2.5, '\rightarrow');
latextitle('$$\det F(s)$$ f\"ur die Reaktoren mit ZRF');
latexylabel('Im');
rotateY;
latexxlabel('Re');




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