%   Script ZRFMAGN.M zur Herstellung von Bildern
%
%   Entwurf einer Zustandsr\"uckf\"uhrung für die Magnetschwebebahn
%
% verwendet BSPMAGN.M, EPSBILD.M
%
% J. Lunze
% 15.8.1996
% für 3. Auflage: 25.3.2004
% für 5. Auflage: 18.10.2007
% für 6. Auflage: 3.2.2010
% für 7. Auflage: 2Te.2012
% für 8. Auflage: 3.4.2014
% für 9. Auflage: 5.12.2017
% für 10. Auflage: 12.9.2019
%
%
clear
close all
Bspmagn;    % definiert Transrapid sowie A, B, C, D

%%   Analyse der Eigenbewegung des Regelkreises
%
figure(1);
sigmazrf=eig(A);
sigmazrf(3)=-5;           %=-sigmazrf(3)/3;  % instabiler Eigenwert
K=place(A, B, sigmazrf);
x0=[5; 0; 0; 0; 0];       % Anfangsauslenkung: 5 mm
subplot(6,1,[1,2]);
Te=0.8;
T=[0:0.01:Te];
Regelkreis1=ss(A-B*K, B, [1 0 0 0 0], 0);
[X1, T1, X]=initial(Regelkreis1, x0, T);
plot(T, X1, 'b', [0 Te], [0 0], 'k:');
xleer;
latexylabel('$$x_1$$ in mm');
rotateY;
yticks([0 4]);
axis([0 Te -1 5.2]);                    % Massstab wie in BOBMAGN.M
latextitle('Eigenbewegungen der geregelten Magnetschwebebahn');
subplot(6,1,[3,4]);
plot(T, X(:, 2), 'b', [0 Te], [0 0], 'k:');
xleer;
latexylabel('$$x_2$$ in mm');
rotateY;
axis([0 Te -0.2 1.7]);
yticks([0 1.5]);
subplot(6,1,[5,6]);
Regelkreis3=ss(A-B*K, B, K, 0);
U=initial(Regelkreis3, x0, T);
plot(T, U, 'b', [0 Te], [0 0], 'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$u$$ in V');
rotateY;
axis([0  Te -10000 30000]);
yticks([0 20000]);
xticks([0 0.2 0.4 0.6 0.8]);


Regelkreis4=ss(A-B*K, B, [0 1 0 0 0], 0);
[X2, T, X]=initial(Regelkreis4, x0, T);    % x_2 zum späteren Vergleich

%%  Reglerentwurf mit veränderten Polen
%
sigmazrf(4)=-3+4.844i;
sigmazrf(5)=-3-4.844i;
sigmazrf
K=place(A, B, sigmazrf);
figure(2);
% 
subplot(6,1,[1,2]);
Regelkreis1=ss(A-B*K, B, [1 0 0 0 0], 0);
[X1n, T1, X]=initial(Regelkreis1, x0, T);
plot(T, X1n, 'b', T, X1, 'b--', [0 Te], [0 0], 'k:');
xleer;
latexylabel('$$x_1$$ in mm');
rotateY;
yticks([0 4]);
axis([0 Te -1 5.2]);                    % Massstab wie in BOBMAGN.M
latextitle('Eigenbewegungen bei der zweiten Regelung');
%
subplot(6,1,[3,4]);
Regelkreis4=ss(A-B*K, B, [0 1 0 0 0], 0);
[X2n, T, X]=initial(Regelkreis4, x0, T);
plot(T, X2n, 'b', T, X2, 'b--', [0 Te], [0 0], 'k:');
xleer;
latexylabel('$$x_2$$ in mm');
rotateY;
axis([0 Te -0.2 1.7]);
yticks([0 1.5]);
%
subplot(6,1,[5,6]);
Regelkreis3=ss(A-B*K, B, K, 0);
[Un, T, X]=initial(Regelkreis3, x0, T);
plot(T, Un, 'b', T, U, 'b--', [0 Te], [0 0], 'k:');
latexxlabel('$$t$$ in s');
latexylabel('$$u$$ in V');
rotateY;
axis([0  Te -10000 30000]);
yticks([0 20000]);
xticks([0 0.2 0.4 0.6 0.8]);


%%  Reglerentwurf mit veränderten Polen
%
sigmazrf=[-20+20i -20-20i, -5 -1.253+4.8442i -1.253-4.8442i]
K=place(A, B, sigmazrf);
figure(3);

subplot(6,1,[1,2]);
Regelkreis1=ss(A-B*K, B, [1 0 0 0 0], 0);
[X1n, T1, X]=initial(Regelkreis1, x0, T);
plot(T, X1n, 'b', T, X1, 'b--', [0 Te], [0 0], 'k:');
xleer;
latexylabel('$$x_1$$ in mm');
rotateY;
yticks([0 4]);
axis([0 Te -1 5.2]);                    % Massstab wie in BOBMAGN.M
latextitle('Eigenbewegungen bei der dritten Regelung');
%
subplot(6,1,[3,4]);
plot(T, X2n, 'b', T, X2, 'b--', [0 Te], [0 0], 'k:');
xleer;
latexylabel('$$x_2$$ in mm');
rotateY;
axis([0 Te -0.2 1.7]);
yticks([0 1.5]);
%
subplot(6,1,[5,6]);
Regelkreis3=ss(A-B*K, B, K, 0);
[Un, T]=initial(Regelkreis3, x0, T);
plot(T, Un, 'b', T, U, 'b--', [0 Te], [0 0], 'k:');
axis([0  Te -10000 30000]);
yticks([0 20000]);
xticks([0 0.2 0.4 0.6 0.8]);
latexxlabel('$$t$$ in s');
latexylabel('$$u$$ in V');
rotateY;


%%   Stabilitätsanalyse mit dem Nyquistkriterium
%
figure(4)
%  Regelkreis besteht aus Regelstrecke (A, B, I, 0) (wobei y=x ist)
%  Zustandsr\"uckf\"uhrung u = - Kx
%  ergibt Darstellung der R\"uckf\"uhrdifferenzmatrix durch
%           xdot = A x + BKy
%              y = I x + I y
Rueckfuehrdifferenzmatrix=ss(A, B*K, eye(n, n), eye(n, n));
[Real, Imag]=nyquist(Rueckfuehrdifferenzmatrix);
for k=1:length(Real)
    detF(k)=det(Real(:,:,k)+i*Imag(:,:,k));
end
plot(real(detF), imag(detF), 'b', real(detF), -imag(detF),'b--');
hold on;
plot([0 0],[-30 60],'k:', [-30 60],[0 0],'k:');
axis([-1.8 1.8 -1.8 1.8]);
yticks([-1 0 1]);
xticks([-1 0 1]);
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Ortskurve von det F der Magnetschwebebahn mit ZRF');
text(0.6, 1.3, '\leftarrow');
text(0.6,  -1.3, '\rightarrow');
plot(1, 0, 'b*');
plot(0, 0, 'b+');
text(1.1, 0.1, '1');
plot(detF(1), 0, 'bo');
hold off;


%%
%   .EPS- und .PS- Files
%
Epsbildc15('zrfmagn');
