%   Script DSTBREA.M zur Herstellung von Bildern
%
%   Steuerbarkeit zweier R\"uhrkesselreaktoren
%
% verwendet BSPREA3.M, EPSBILD.M
%
% J. Lunze
% 6.1.1997               Daten v. 30.9.1996  aus STBREAK.M
% für 3. Auflage: 15.3.2004
% für 5. Auflage: 18.10.2007
% für 7. Auflage 31.3.2012
% für 10. Auflage: 9.9.2019
%
%
clear
close all

Bsprea3;        % -> (A, B, C, D)
Ta=0.3;
[Ad, Bd]=c2d(A, B, Ta);



%%    Umsteuerung des Reaktors mit Abtastzeit Ta=0,3
%
figure(1);
Te=n;
BB=Bd*Bd';
W=zeros(n,n);
Adj=eye(n,n);
for k=0:Te-1
  W=W + Adj*BB*Adj';
  Adj=Adj*Ad;
end;
%
Th=Te+2;
T=[0:Th];
xe=[1 5]';
x0=[0 0]';
xdiff=Ad^Te*x0-xe;
Wx=inv(W)*xdiff;
for k=1:Th
   U(k,:)=-Bd'*expm(A*(Te-k)*Ta)'*Wx;
end;
U(k+1,:)=U(k,:);
%U(k+2,:)=U(k,:);
%U(k+3,:)=U(k,:);
[Y, X]=dlsim(Ad, Bd, eye(n,n), zeros(n,m), U);
subplot(6,1,[1,2]);
stairs([-1 T], [0 Y(:,1)'], 'Color','b','LineWidth',1.5);
hold on;
plot(T, Y(:,1), 'bo', [0 0], [-100 100], 'k:', [2 2], [-100 100], 'k:', [-1 4 ], [0 0], 'k:');
%text(Te, xe(1,1), 'b*');
latextitle('Umsteuerung des Reaktors $$x_0=(0 0) \rightarrow x_e=(1 5)$$ bei $$T_a=0.3$$ min');
set(gca, 'XTick', [-1:3]);
xleer;
latexylabel('$$x_1$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
hold on
plot(2, 1, 'b*');
axis([-0.5 3.2 -20 20]);
yticks([-20 0 20]);
hold off;
subplot(6,1,[3,4]);
stairs([-1 T], [0 Y(:,2)'], 'Color','b','LineWidth',1.5);
hold on;
plot(T, Y(:,2), 'bo', [0 0], [-100 100], 'k:', [2 2], [-100 100], 'k:', [-1 4 ], [0 0], 'k:');
%text(Te, xe(2,1), 'b*');
set(gca, 'XTick', [-1:3]);
xleer;
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
hold on
plot(2, 5, 'b*');
axis([-0.5 3.2 -10 10]);
yticks([-10 0 10]);
hold off;
subplot(6,1,[5,6]);
stairs([-1 T], [0 U'], 'Color','b','LineWidth',1.5);
hold on;
plot([-1 4], [0 0], 'k:', [2 2], [-1000 1000], 'k:', [0 0], [-1000 1000], 'k:');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
set(gca, 'XTick', [-1:3]);
latexxlabel('$$k$$');
axis([-0.5 3.2 -160 160]);
yticks([-150 0 150]);
xticks([0 1 2 3]);
hold off;


%%   Wiederholung bei Abtastzeit Ta=0,8
%
%
figure(2);
Ta=0.8;
[Ad, Bd]=c2d(A, B, Ta);
BB=Bd*Bd';
W=zeros(n,n);
Adj=eye(n,n);
for k=0:Te-1
  W=W + Adj*BB*Adj';
  Adj=Adj*Ad;
end;
%
xdiff=Ad^Te*x0-xe;
Wx=inv(W)*xdiff;
for k=1:Th/Ta-1
   U(k,:)=-Bd'*expm(A*(Te-k)*Ta)'*Wx;
end;
U(k+1,:)=U(k,:);
%U(k+2,:)=U(k,:);
%U(k+3,:)=U(k,:);
[Y, X]=dlsim(Ad, Bd, eye(n,n), zeros(n,m), U);
subplot(6,1,[1,2]);
stairs([-1 T], [0 Y(:,1)'], 'Color','b','LineWidth',1.5);
hold on;
plot(T, Y(:, 1),'bo', [0 0], [-100 100], 'k:', [2 2], [-100 100], 'k:', [-1 4 ], [0 0], 'k:');
%text(Te, xe(1,1), 'b*');
latextitle('Umsteuerung des Reaktors $$x_0=(0 0) \rightarrow x_e=(1 5)$$ mit $$T_a$$=0,8 min');
set(gca, 'XTick', [-1:3]);
xleer;
latexylabel('$$x_1$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
hold on
plot(2, 1, 'b*');
axis([-0.5 3.2 -20 20]);
yticks([-20 0 20]);
hold off;
subplot(6,1,[3,4]);
stairs([-1 T], [0 Y(:,2)'], 'Color','b','LineWidth',1.5);
hold on;
plot(T, Y(:, 2), 'bo', [0 0], [-100 100], 'k:', [2 2], [-100 100], 'k:', [-1 4 ], [0 0], 'k:');
%text(Te, xe(2,1), 'b*');
set(gca, 'XTick', [-1:3]);
xleer;
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
hold on
plot(2, 5, 'b*');
axis([-0.5 3.2 -10 10]);
yticks([-10 0 10]);
hold off;
subplot(6,1,[5,6]);
stairs([-1 T], [0 U'], 'Color','b','LineWidth',1.5);
hold on;
plot([-1 4], [0 0], 'k:', [2 2], [-1000 1000], 'k:', [0 0], [-1000 1000], 'k:');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
set(gca, 'XTick', [-1:3]);
latexxlabel('$$k$$');
axis([-0.5 3.2 -55 55]);
yticks([-50 0 50]);
xticks([0 1 2 3]);
hold off;


%%
%
%    Verhalten des kontinuierlichen Systems bei Ta=0,8
%
%
figure(3);
dt=0.02;
Tk=Ta*[0:dt:Th+1-dt];
for k=1:Th+1
for k1=1:1/dt
    Uk((k-1)/dt+k1,:)=U(k,:);
end;
end;
kontSystem=ss(A, B, eye(n,n), zeros(n,m));
[Y, Tk, X]=lsim(kontSystem, Uk, Tk);
subplot(6,1,[1,2]);
plot([-Ta Tk'], [0 Y(:,1,1)'], 'b');
hold on;
plot([0 0], [-100 100], 'k:', [2*Ta 2*Ta], [-100 100], 'k:', [-1 3 ], [0 0], 'k:');
plot([Ta Ta], [-100 100], 'k:', [2*Ta 2*Ta], [-100 100], 'k:', [3*Ta 3*Ta], [-100 100], 'k:');
%text(Te, xe(1,1), 'b*');
latextitle('Verhalten des kontinuierlichen Reaktors bei $$T_a$$=0,8 min');
xleer;
latexylabel('$$x_1$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.5 3.2 -20 20]);
yticks([-20 0 20]);
plot(2*Ta, 1, 'b*');
hold off;
subplot(6,1,[3,4]);
plot([-Ta Tk'], [0 Y(:,2)'], 'b');
hold on;
plot([0 0], [-100 100], 'k:', [2*Ta 2*Ta], [-100 100], 'k:', [-1 3.5 ], [0 0], 'k:');
plot([Ta Ta], [-100 100], 'k:', [2*Ta 2*Ta], [-170 170], 'k:', [3*Ta 3*Ta], [-170 170], 'k:');
%text(Te, xe(2,1), 'b*');
xleer;
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.5 3.2 -10 10]);
yticks([-10 0 10]);
plot(2*Ta, 5, 'b*');
hold off;
subplot(6,1,[5,6]);
stairs([-Ta Tk'], [0 Uk'], 'Color','b','LineWidth',1.5);
hold on;
plot([-1 3.5], [0 0], 'k:', Ta*[2 2], [-1000 1000], 'k:', [0 0], [-1000 1000], 'k:');
plot([Ta Ta], [-170 170], 'k:', [2*Ta 2*Ta], [-170 170], 'k:', [3*Ta 3*Ta], [-170 170], 'k:');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([-0.5 3.2 -55 55]);
yticks([-50 0 50]);
xticks([0 1 2 3]);
hold off;


%%   Wiederholung bei Abtastzeit Ta=1.1 min
%
%
figure(4);
Ta=1.1;
[Ad, Bd]=c2d(A, B, Ta);
BB=Bd*Bd';
W=zeros(n,n);
Adj=eye(n,n);
for k=0:Te-1
  W=W + Adj*BB*Adj';
  Adj=Adj*Ad;
end;
%
xdiff=Ad^Te*x0-xe;
Wx=inv(W)*xdiff;
for k=1:Th/Ta-1
   U(k,:)=-Bd'*expm(A*(Te-k)*Ta)'*Wx;
end;
U(k+1,:)=U(k,:);
%U(k+2,:)=U(k,:);
%U(k+3,:)=U(k,:);
Y=dlsim(Ad, Bd, eye(n,n), zeros(n,m), U);
%
subplot(6,1,[1,2]);
stairs([-1 T], [0 Y(:,1)'], 'Color','b','LineWidth',1.5);
hold on;
plot(T, Y(:,1),'bo', [0 0], [-100 100], 'k:', [2 2], [-100 100], 'k:', [-1 4], [0 0], 'k:');
plot(2, 1, 'b*');
latextitle('Umsteuerung des Reaktors $$x_0=(0 0) \rightarrow x_e=(1 5)$$ mit $$T_a$$=1,1 min');
set(gca, 'XTick', [-1:3]);
xleer;
latexylabel('$$x_1$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.5 3.2 -20 20]);
yticks([-20 0 20]);
hold off;
%
subplot(6,1,[3,4]);
stairs([-1 T], [0 Y(:,2)'], 'Color','b','LineWidth',1.5);
hold on;
plot(T, Y(:,2), 'bo', [0 0], [-100 100], 'k:', [2 2], [-100 100], 'k:', [-1 4], [0 0], 'k:');
plot(2, 5, 'b*');
set(gca, 'XTick', [-1:3]);
xleer;
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.5 3.2 -10 10]);
yticks([-10 0 10]);
hold off;
subplot(6,1,[5,6]);
stairs([-1 T], [0 U'], 'Color','b','LineWidth',1.5);
hold on;
plot([-1 4], [0 0], 'k:', [2 2], [-1000 1000], 'k:', [0 0], [-1000 1000], 'k:');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
set(gca, 'XTick', [-1:3]);
latexxlabel('$$k$$');
axis([-0.5 3.2 -55 55]);
yticks([-50 0 50]);
xticks([0 1 2 3]);
hold off;


%%
%
%    Verhalten des kontinuierlichen Systems bei Ta=1.1
%
%
figure(5);
dt=0.02;
Tk=Ta*[0:dt:Th+1-dt];
for k=1:Th+1
for k1=1:1/dt
   Uk((k-1)/dt+k1,:)=U(k,:);
end;
end;
kontSystem=ss(A, B, eye(n,n), zeros(n,m));
[Y, Tk, X]=lsim(kontSystem, Uk, Tk);
subplot(6,1,[1,2]);
plot([-Ta Tk'], [0 Y(:,1)'], 'b');
hold on;
plot([0 0], [-100 100], 'k:', [2*Ta 2*Ta], [-100 100], 'k:', [-1 4], [0 0], 'k:');
plot([Ta Ta], [-100 100], 'k:', [2*Ta 2*Ta], [-100 100], 'k:', [3*Ta 3*Ta], [-100 100], 'k:');
plot(2*Ta, 1, 'b*');
latextitle('Verhalten des kontinuierlichen Reaktors bei $$T_a$$=1,1 min');
xleer;
latexylabel('$$x_1$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.5 3.2 -20 20]);
yticks([-20 0 20]);
hold off;
subplot(6,1,[3,4]);
plot([-Ta Tk'], [0 Y(:,2)'], 'b');
hold on;
plot([0 0], [-100 100], 'k:', [2*Ta 2*Ta], [-100 100], 'k:', [-1 4], [0 0], 'k:');
plot([Ta Ta], [-100 100], 'k:', [2*Ta 2*Ta], [-100 100], 'k:', [3*Ta 3*Ta], [-100 100], 'k:');
plot(2*Ta, 5, 'b*');
xleer;
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
axis([-0.5 3.2 -10 10]);
yticks([-10 0 10]);
hold off;
subplot(6,1,[5,6]);
stairs([-Ta Tk'], [0 Uk'], 'Color','b','LineWidth',1.5);
hold on;
plot([-1 4], [0 0], 'k:', Ta*[2 2], [-1000 1000], 'k:', [0 0], [-1000 1000], 'k:');
plot([Ta Ta], [-150 150], 'k:', [2*Ta 2*Ta], [-150 150], 'k:', [3*Ta 3*Ta], [-150 150], 'k:');
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm l}$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([-0.5 3.2 -55 55]);
yticks([-50 0 50]);
xticks([0 1 2 3 4]);
hold off;
%

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