%   Script STBREAK.M zur Herstellung von Bildern
%
%   Steuerbarkeit zweier R\"uhrkesselreaktoren
%
% verwendet BSPREA3.M, EPSBILD.M
%
% J. Lunze
% 30.9.1996
% für 3. Auflage: 15.3.2004
% für 5. Auflage: 2.1.2008
% für 7. Auflage 25.5.2012
% für 10. Auflage: 10.9.2019
%
%
clear
close all
Bsprea3;


%%        Berechnung von W(Te)
%
figure(1);
Te=2;                  % Endzeit für Zustands\"uberf\"uhrung
dt=Te/1000;
BB=B*B';
expA=expm(A*Te);
W=BB+expA*BB*expA';
W=W/2;
for k=1:Te/dt-1
  t=k*dt;
  expAt=expm(A*t);
  W=W + expAt*BB*expAt';
end;
W=W*dt;
%
Th=3;                   % Simulationshorizont
T=[0:dt:Th];
xe=[1 5]';
x0=[0 0]';
xdiff=expA*x0-xe;
Wx=inv(W)*xdiff;
for k=1:Th/dt+1
   t=(k-1)*dt;
   U(k,:)=-B'*expm(A*(Te-t))'*Wx;
end;
ReaktorX=ss(A, B, eye(2,2), zeros(2,1));   % y=x
[Y, T, X]=lsim(ReaktorX, U, T, x0);
%
subplot(6,1,[1,2]);
plot(T, Y(:,1), 'b',[-0.3 0], [0 0],'b');
hold on;
plot([0 0], [-10 12], 'k:', [2 2], [-10 12], 'k:', [-0.3 3 ], [0 0], 'k:');
plot(Te, 1, 'bo');
latextitle('Umsteuerung des Reaktors $$x_0=(0 0) \rightarrow x_e=(1 5)$$');
xleer;
latexylabel('$$x_1$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.3 2.5 -5 12]);
yticks([0 5 10]);
hold off;
%
subplot(6,1,[3,4]);
plot(T, Y(:,2), 'b', [-0.3 0], [0 0], 'b');
hold on;
plot(Te, 5, 'bo');
plot([0 0], [-10 12], 'k:', [2 2], [-10 12], 'k:', [-0.3 3 ], [0 0], 'k:');
%text(Te, xe(2,1), 'b*');
xleer;
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.3 2.5 -5 12]);
yticks([0 5 10]);
hold off;
subplot(6,1,[5,6]);
plot(T, U(:,1,1), 'b', [-0.3 0 0], [0 0 U(1,1)], 'b');
hold on;
plot([-0.3 3], [0 0], 'k:', [2 2], [-150 100], 'k:', [0 0], [-150 100], 'k:');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([-0.3 2.5 -150 50]);
xticks([0 1 2 3]);
yticks([-100 -50 0 50]);
hold off;


%%
%
%  Dieselbe Zustands\"uberf\"uhrung mit k\"urzerer Endzeit
%
%
figure(2);
%
%        Berechnung von W(Te)
%
Te=1.6;                  % Endzeit für Zustands\"uberf\"uhrung
                         % vgl. dieselben Betrachtungen bei Abtastsystemen
dt=Te/1000;
BB=B*B';
expA=expm(A*Te);
W=BB+expA*BB*expA';
W=W/2;
for k=1:Te/dt-1
  t=k*dt;
  expAt=expm(A*t);
  W=W + expAt*BB*expAt';
end;
W=W*dt;
%
Th=3;                   % Simulationshorizont
T=[0:dt:Th];
xe=[1 5]';
x0=[0 0]';
xdiff=expA*x0-xe;
Wx=inv(W)*xdiff;
U=[];
for k=1:Th/dt+1
   t=(k-1)*dt;
   U(k,:)=-B'*expm(A*(Te-t))'*Wx;
end;
[Y, T, X]=lsim(ReaktorX, U, T, x0);
%
subplot(6,1,[1,2]);
plot(T, Y(:,1), 'b',[-0.3 0], [0 0],'b');
hold on;
plot([0 0], [-10 12], 'k:', [Te Te], [-10 12], 'k:', [-0.3 3 ], [0 0], 'k:');
plot(Te, 1, 'bo');
latextitle('Umsteuerung des Reaktors $$x_0=(0 0) \rightarrow x_e=(1 5)$$');
xleer;
latexylabel('$$x_1$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.3 2.5 -5 12]);
yticks([0 5 10]);
hold off;
subplot(6,1,[3,4]);
plot(T, Y(:,2), 'b', [-0.3 0], [0 0], 'b');
hold on;
plot(Te, 5, 'bo');
plot([0 0], [-10 12], 'k:', [Te Te], [-10 12], 'k:', [-0.3 3 ], [0 0], 'k:');
%text(Te, xe(2,1), 'b*');
xleer;
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.3 2.5 -5 12]);
yticks([0 5 10]);
hold off;
%
subplot(6,1,[5,6]);
plot(T, U(:,1,1), 'b', [-0.3 0 0], [0 0 U(1,1)], 'b');
hold on;
plot([-0.3 3], [0 0], 'k:', [Te Te], [-150 100], 'k:', [0 0], [-150 100], 'k:');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([-0.3 2.5 -150 50]);
xticks([0 1 2 3]);
yticks([-100 -50 0 50]);
hold off;


%%
%
%  Dieselbe Zustands\"uberf\"uhrung mit st\"uckweise konstantem Eingang
%
%
figure(3);
Te=1.6;     % Endzeit
T1=1;       % u(t) = u0 für 0 < t < T1
            % u(t) = u1 für T1 < t < Te
p0=-inv(A)*expm(A*(Te-T1))*B +inv(A)*expm(A*Te)*B;
p1=-inv(A)*B +inv(A)*expm(A*(Te-T1))*B;
P=[p0 p1];
xe=[1 5]';   % wie bisher
x0=[0 0]';
U=inv(P)*(xe-(expm(-A*Te)*x0));
u0=U(1,1);
u1=U(2,1);

Th=Te;                   % Simulationshorizont
T=[0:dt:Th];
k=0;
U=u0*ones(1,length(T));
k1=T1/dt;
for k=k1:length(T)
    U(1,k)=u1;
end
[Y, T, X]=lsim(ReaktorX, U, T, x0);
subplot(6,1,[1,2]);
plot(T, Y(:,1), 'b', [-0.3 0], [0 0], 'b');
hold on;
plot([0 0], [-10 12], 'k:', [Te Te], [-10 12], 'k:', [-0.3 3 ], [0 0], 'k:');
plot(Te, 1, 'bo');
latextitle('Umsteuerung des Reaktors $$x_0=(0 0) \rightarrow x_e=(1 5)$$');
xleer;
latexylabel('$$x_1$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.3 2.5 -5 12]);
yticks([0 5 10]);
hold off;
subplot(6,1,[3,4]);
plot(T, Y(:,2), 'b',[-0.3 0], [0 0], 'b');
hold on;
plot(Te, 5, 'bo');
plot([0 0], [-10 12], 'k:', [Te Te], [-10 12], 'k:', [-0.3 3 ], [0 0], 'k:');
%text(Te, xe(2,1), 'b*');
xleer;
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.3 2.5 -5 12]);
yticks([0 5 10]);
hold off;
subplot(6,1,[5,6]);
plot(T, U, 'b', [-0.3 0 0], [0 0 U(1,1)], 'b');
hold on;
plot([-0.3 3], [0 0], 'k:', [Te Te], [-150 100], 'k:', [0 0], [-150 100], 'k:');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([-0.3 2.5 -150 50]);
xticks([0 1 2 3]);
yticks([-100 -50 0 50]);
hold off;

%%   Bewegung im Zustandsraum
%
figure(4)             %  Bewegung im Zustandsraum
plot(X(:,1), X(:,2), 'b')
hold on
latexxlabel('$$x_1$$ in $${\rm mol/l}$$');   % ... sonst wird Maßeinheit im eps-Bild abgeschnitten
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-1 11 -1 7]);
latextext(x0(1,1)-0.5, x0(2,1)-0.5, '$$\mathbf{x}_0=\mathbf{0}$$');
latextext(xe(1,1)-0.5, xe(2,1)-0.5, '$$\mathbf{x}_{\rm e}=(1\;\;5)^{\rm T}$$');
plot(x0(1,1), x0(2,1), 'bo');
plot(xe(1,1), xe(2,1), 'bo');
plot(X(k1,1), X(k1,2), 'b*');
annotation('textarrow',[0.3839 0.8289],[0.3284 0.8155]);
annotation('textarrow',[0.8214 0.439],[0.8343 0.744]);
hold off
xticks([0 5 10]);
yticks([0 5 10]);


%%  Richtungsvektoren bei unterschiedlicher Zeit t_1
%
figure(5)
for T1=[0.01, 0.05, 0.1, 0.2, 0.5, 1]
    p0=-inv(A)*expm(A*(Te-T1))*B +inv(A)*expm(A*Te)*B;
    p1=-inv(A)*B +inv(A)*expm(A*(Te-T1))*B;
    P=[p0 p1];  
    x0=[0 0]';
    U=inv(P)*(xe-(expm(-A*Te)*x0));
    u0=U(1,1);
    u1=U(2,1);
    x1=expm(A*T1)*x0 - inv(A)*B*u0+inv(A)*expm(A*T1)*B*u0;
    plot([x0(1,1), x1(1,1), xe(1,1)], [x0(2,1), x1(2,1), xe(2,1)], 'b');
    hold on
    plot(x1(1,1), x1(2,1), 'b*');
    %text('Interpreter','latex', 'String','$$t_1=$$ num2str(T1)', 'Position', [x1(1,1)+0.5, x1(2,1)]);
    latextext( x1(1,1)+0.5, x1(2,1), num2str(T1));
end
latexxlabel('$$x_1$$ in $${\rm mol/l}$$');
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
plot(x0(1,1), x0(2,1), 'bo');
plot(xe(1,1), xe(2,1), 'bo');
latextext(x0(1,1)-0.5, x0(2,1)-0.5, '$$\mathbf{x}_0=\mathbf{0}$$');
latextext(xe(1,1)-0.5, xe(2,1)+1.5, '$$\mathbf{x}_{\rm e}=(1\;\;5)^{\rm T}$$');
axis([-1 23 -1 7]);
xticks([0 10  20]);
yticks([0 5 10 20]);
hold off

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