%   Script NichtminReaktor.M
%
%   Reaktor mit minimalphasigem und nichtminimalphasigem Verhalten
%
% J. Lunze
% 26.1.2010
%
% für 8. Auflage: 30.1.2010
% für 9. Auflage: 13.2.2012
% für 10. Auflage: 5.4.2014
% für 12. Auflage: 23.03.2020
%
%   
clear
close all

%   Reaktormodelle
F=0.1;      % m^3/min
V1=1;       % m^3
V2=0.7;     % m^3
A=[-F/V1  0; F/V2  -2*F/V2];
b1=[F/V1; F/V2];
b2=[F/V1; -F/V2];
c=[0  1];
d=0;
Reaktor1=ss(A, b1, c, d);
Reaktor2=ss(A, b2, c, d);

%  Transformation in E/A-Normalform
%
%    Reaktor 1
%
r=1;
n=2;
SS1=ctrb(A, b1);
sRT1=[0 1]*inv(SS1);
TE1=inv([c; sRT1]);
AE1=inv(TE1)*A*TE1;
bE1=inv(TE1)*b1;
cE1=c*TE1;
Reaktor1EA=ss(AE1, bE1, cE1, d);
%   Nulldynamik
ANN1=AE1(r+1:n, r+1:n);
bq1=c*b1;
bN1=AE1(r+1:n,1);            % Eingang ist y(t)=w(t)
eNT1=AE1(r,r+1:n);
cN1=-1/bq1*eNT1;              % Ausgabe ist uVN= -1/bq eNT*xN = uVI
dN1=0;
Nulldynamik1=ss(ANN1, bN1, cN1, dN1);
eIT1=AE1(r, 1:r);

%
%
%  Sollwert
we=4;
Tbeg=0;
Tend=20;
T=Tbeg:0.01:Tend;
We=ones(length(T),1)*we;
x0=[2; we];


%%     Reaktor 1
%
%  Berechnung der Stellgröße für Anfangszustand x0
%  uV1(t) = -1/bq1 eNT1*xN - 1/bq1 eIT1 we
%       = uVN1 + uVI1
figure(1)
uVI1=-1/bq1*eIT1*We;
xN01=[0 1]*inv(TE1)*x0;
[uVN1, T, XN]=lsim(Nulldynamik1, We, T, xN01);  % Bestimmung von uVN mit XN0=0
uV1=uVN1+uVI1;
subplot(6,1,[1,2])
plot(T, uV1, 'b');
hold on
plot([Tbeg 0 0], [0 0 uV1(1,1)], 'k');
axis([Tbeg, Tend, -0.5, 7]);
yticks([0 2 4 6 8]);
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm m^3}\;$$');
rotateY;
hold off
xleer;
%
subplot(6,1,[3,4])
[Y1, T, X1]=lsim(Reaktor1,uV1,T,x0);
plot(T, Y1, 'b', [-5 0 0],[0 0 Y1(1,1)], 'k');
hold on
plot(0, we, 'bo');
latexylabel('$$y$$ in $$\frac{\rm mol}{\rm m^3}\;$$');
rotateY;
axis([Tbeg, Tend, -0.5 5]);
yticks([0 2 4 6 8]);
xticks([0 5 10 15 20]);
latexxlabel('$$t$$ in min');

%%  Darstellung im Zustandsraum
%
figure(2)
subplot(3,4,[1 2 5 6]);
plot(X1(:, 1), X1(:, 2), 'b');
hold on;
plot([-5 6],[0 0],'k:', [0 0],[-5 5],'k:');
latexxlabel('$$x_1$$ in $$\frac{\rm mol}{\rm m^3}$$');
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm m^3}$$');
latextext(2.5, 3, '$$\rightarrow$$');
rotateY;
plot(x0(1,1), x0(2,1), 'b*');
us1=we/dcgain(Reaktor1);
xe=-inv(A)*b1*us1;
plot(xe(1,1), xe(2,1), 'bo');
axis([-1 6 -0.5 5]);
xticks([0 2 4 6]);
yticks([0 2 4 6]);
hold off;
 
%%    Reaktor 2
%
SS2=ctrb(A, b2);
sRT2=[0 1]*inv(SS2);
TE2=inv([c; sRT2]);
AE2=inv(TE2)*A*TE2;
bE2=inv(TE2)*b2;
cE2=c*TE2;
Reaktor2EA=ss(AE2, bE2, cE2, d);
%   Nulldynamik
ANN2=AE2(r+1:n, r+1:n);
bq2=c*b2;
bN2=AE2(r+1:n,1);            % Eingang ist y(t)=w(t)
eNT2=AE2(r,r+1:n);
cN2=-1/bq2*eNT2;              % Ausgabe ist uVN= -1/bq eNT*xN = uVI
dN2=0;
Nulldynamik2=ss(ANN2, bN2, cN2, dN2);
eIT2=AE2(r, 1:r);
vT2=[eIT2 -1];
%  Berechnung der Stellgröße für Anfangszustand x0
%  u(t) = -1/bq eNT*xN - 1/bq eIT1 we
%       = uVN + uVI
figure(3)
uVI2=-1/bq2*eIT2*We;
xN02=[0 1]*inv(TE2)*x0;
uVN2=lsim(Nulldynamik2, We, T, xN02);  % Bestimmung von uVN mit XN0=0
uV2=uVN2+uVI2;
subplot(6,1,[1,2]);
plot(T, uV2, 'b');
hold on
plot([Tbeg 0 0], [0 0 uV2(1,1)], 'b');
axis([Tbeg, Tend, -20, 1]);
yticks([-20 -10 0]);
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm m^3}\;$$');
rotateY;
hold off
xleer;
%
subplot(6,1,[3,4])
[Y2, T, X2]=lsim(Reaktor2,uV2,T,x0);
plot(T, Y2, 'b', [-5 0 0],[0 0 Y2(1,1)],'b');
hold on
plot(0, we, 'bo');
latexylabel('$$y$$ in $$\frac{\rm mol}{\rm m^3}\;$$');
rotateY;
axis([Tbeg, Tend, -0.5 5])
xticks([0 5 10 15 20]);
yticks([0 2 4]);
latexxlabel('$$t$$ in min');


%%  Darstellung im Zustandsraum
%
figure(2)
subplot(3,4,[3 4 7 8]);
plot(X2(:, 1), X2(:, 2));
hold on;
plot([-5 6],[0 0],'k:', [0 0],[-5 5],'k:');
latexxlabel('$$x_1$$ in $$\frac{\rm mol}{\rm m^3}$$');
latextext(0.5, 3, '$$\leftarrow$$');
yleer;
plot(xe(1,1), xe(2,1), 'bo');
plot(x0(1,1), x0(2,1), 'b*');
axis([-1 6 -0.5 5]);
xticks([0 2 4 6 8]);
hold off;
 

%%    Zusammengefasste Darstellung
%
figure(4)
subplot(6,1,[1,2])
plot(T, uV2,'b--');
hold on
plot([Tbeg 0 0], [0 0 uV2(1,1)],'b--');
plot(T, uV1, 'b');
plot([Tbeg 0 0], [0 0 uV1(1,1)]), 'b';
axis([Tbeg, Tend, -20, 7]);
yticks([-20 -10 0]);
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm m^3}\;$$');
rotateY;
xleer;
hold off;
%
subplot(6,1,[3,4])
[Y2, T, X2]=lsim(Reaktor2,uV2,T,x0);
plot(T, Y2, 'b',  [-5 0],[Y2(1,1) Y2(1,1)], 'b');
hold on
plot(0, we, 'bo');
latexylabel('$$y$$ in $$\frac{\rm mol}{\rm m^3}\;$$');
rotateY;
axis([Tbeg, Tend, -0.5 5])
xticks([0 5 10 15 20]);
yticks([0 2 4]);
latexxlabel('$$t$$ in min');
 

%%  Dasselbe für we=0
% 
%  Sollwert
we=0;
We=ones(length(T),1)*we;
x0=[-3; we];

%
%     Reaktor 1
%
%  Berechnung der Stellgröße für Anfangszustand x0
%  uV1(t) = -1/bq1 eNT1*xN - 1/bq1 eIT1 we
%       = uVN1 + uVI1
figure(5)
uVI1=-1/bq1*eIT1*We;
xN01=[0 1]*inv(TE1)*x0;
[uVN1, T, XN]=lsim(Nulldynamik1, We, T, xN01);  % Bestimmung von uVN mit XN0=0
uV1=uVN1+uVI1;
subplot(6,1,[1,2]);
plot(T, uV1, 'b');
hold on
plot([Tbeg 0 0], [0 0 uV1(1,1)], 'b');
axis([Tbeg, Tend, -0.5, 3.5]);
yticks([0 1 2 3]);
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm m^3}$$');
rotateY;
hold off
xleer;
%
subplot(6,1,[3,4])
[Y1, T, X1]=lsim(Reaktor1,uV1,T,x0);
plot(T, Y1, 'b',[-5 0 0],[0 0 Y1(1,1)], 'b');
hold on
plot(0, we, 'bo');
latexylabel('$$y$$ in $$\frac{\rm mol}{\rm m^3}$$');
rotateY;
axis([Tbeg, Tend, -1 1])
xticks([0 5 10 15 20]);
yticks([-1 0 1]);
latexxlabel('$$t$$ in min');

%%  Darstellung im Zustandsraum
%
figure(6)
subplot(2,2,1);
plot(X1(:, 1), X1(:, 2));
hold on;
plot([-5 5],[0 0],'k:', [0 0],[-5 5],'k:');
latexxlabel('$$x_1$$ in $$\frac{\rm mol}{\rm m^3}$$');
latexylabel('$$x_2$$ in $$\frac{\rm mol}{\rm m^3}$$');
rotateY;
plot(x0(1,1), x0(2,1), 'b*');
us1=we/dcgain(Reaktor1);
xe=-inv(A)*b1*us1;
plot(xe(1,1), xe(2,1), 'bo');
axis([-5 1 -3.5 3.5]);
xticks([-4 -2 0]);
yticks([-2 0 2]);
axis('square');
hold off;


%%   Reaktor 2
%
figure(7)
uVI2=-1/bq2*eIT2*We;
xN02=[0 1]*inv(TE2)*x0;
uVN2=lsim(Nulldynamik2, We, T, xN02);  % Bestimmung von uVN mit XN0=0
uV2=uVN2+uVI2;
subplot(6,1,[1,2]);
plot(T, uV2, 'b');
hold on
plot([Tbeg 0 0], [0 0 uV2(1,1)], 'b');
axis([Tbeg, Tend, -3.5, 0.5]);
yticks([-3 -2 -1 0]);
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm m^3}\;$$');
rotateY;
hold off
xleer;
%
subplot(6,1,[3,4])
[Y2, T, X2]=lsim(Reaktor2,uV2,T,x0);
plot(T, Y2, 'b', [-5 0 0],[0 0 Y2(1,1)], 'b');
hold on
plot(0, we, 'bo');
latexylabel('$$y$$ in $$\frac{\rm mol}{\rm m^3}\;$$');
rotateY;
axis([Tbeg, Tend, -1 1])
xticks([0 5 10 15 20]);
yticks([-1 0 1]);
latexxlabel('$$t$$ in min');
hold off


%%  Darstellung im Zustandsraum
%
figure(6)
subplot(2,2,2);
plot(X2(:, 1), X2(:, 2), 'b');
hold on;
plot([-5 5],[0 0],'k:', [0 0],[-5 5],'k:');
latexxlabel('$$x_1$$ in $$\frac{\rm mol}{\rm m^3}$$');
yleer;
plot(xe(1,1), xe(2,1), 'bo');
plot(x0(1,1), x0(2,1), 'b*');
axis([-5 1 -3.5 3.5]);
xticks([-4 -2 0]);
axis('square');
hold off;
 

%%                  
%
%   Zusammengefasste Darstellung
%
%   
figure(8)
subplot(6,1,[1,2]);
plot(T, uV2,'b--');
hold on
plot([Tbeg 0 0], [0 0 uV2(1,1)],'b--');
plot(T, uV1, 'b');
plot([Tbeg 0 0], [0 0 uV1(1,1)], 'b');
axis([Tbeg, Tend, -4, 4]);
yticks([-4 0 4]);
latexylabel('$$u$$ in $$\frac{\rm mol}{\rm m^3}\;$$');
rotateY;
xleer;
hold off;
%
subplot(6,1,[3,4])
[Y2, T, X2]=lsim(Reaktor2,uV2,T,x0);
plot(T, Y2, 'b', [-5 0],[Y2(1,1) Y2(1,1)], 'b');
hold on
plot(0, we, 'bo');
latexylabel('$$y$$ in $$\frac{\rm mol}{\rm m^3}\;$$');
rotateY;
axis([Tbeg, Tend, -1 1]);
xticks([0 5 10 15 20]);
yticks([-1 0 1]);
latexxlabel('$$t$$ in min');



%% Herstellung der Bilddateien *.eps
%
Epsbildc15('nichtminReak');
