%   Script OFENSTB.M
%
%   Steuerbarkeit eines Ofens
%
% J. Lunze
% 24.6.2002
% für 2. Auflage: 29.12.2006
% für 3. Auflage: 1.12.2011
% für 5. Auflage: 11.01.2020
%
%
clear
close all
bspofen;    %  Definition von Ofen (A, B, C, D);
%%    Umsteuerung des Ofens
%
figure(1)
%        Berechnung von W_S(Te)
%
Tes=2.5;                  % Endzeit für Zustandsüberführung
dt1=0.0001;             % für Berechnung von W ist sehr kleine Schritteweite notw.
dt2=0.005;              % Schrittweite für Simulation
Th=Tes*1.2;             % Simulationshorizont
T=[0:dt2:Th];
%
BB=B*B';
expA=expm(A*Tes);
W=BB+expA*BB*expA';
W=W/2;
for k=1:Tes/dt1-1
  t=k*dt1;
  expAt=expm(A*t);
  W=W + expAt*BB*expAt';
end;
W=W*dt1;      % Steuerbarkeits-Gramsche
%
Umin=800;
%  statischer Zustand bei u=Umin (konstant)
x0=-inv(A)*B*Umin;
%  Endzustand
xe1=[8; 5];
%  Steuerung, die Endzustand xe(2,1) aufrechterhält
ksx=-inv(A)*B;
Ukonst=xe1(2,1)/ksx(2,1);
xdiff=expA*x0-xe1;
Wx=inv(W)*xdiff;
for k=1:Tes/dt2+1
   t=(k-1)*dt2;
   U(k,:)=-B'*expm(A*(Tes-t))'*Wx;
end;
for k=1:ceil((Th-Tes)/dt2)
   U(Tes/dt2+k+1,:)=Ukonst;
end
Ulow=min(1.1*min(U), 0);  %  Untere Grenze des Darstellungsbereichs
Uhigh=1.1*max(U);
[Y, T, X]=lsim(Ofen, U, T, x0);
%
subplot(4,1,[1,2]);
plot(T, U(:,1), 'b', [-1 0 0],[Umin Umin U(1,1)], 'b',...
    [0 0],[Ulow, Uhigh],'k:', [Tes, Tes],[Ulow, Uhigh],'b--');
latextitle(['Ofen  $$x_0=$$',num2str(x0'),'$$\rightarrow x_e=$$',num2str(xe1')]);
xleer;
xticksempty([-1 0 1 2 3]);
latexylabel('$$u$$ in W');
rotateY;
axis([-1 Th Ulow, Uhigh]);
yticks([-5000 0 5000]);
%
subplot(4,1,[3,4]);
plot(T, X, 'b');
hold on;
plot([-1 0 0],[x0(1,1), x0(1,1), X(1,1)], 'b',... 
     [-1 0 0],[x0(2,1), x0(2,1), X(1,2)], 'b');
plot([0 0], [-1 20], 'k:', [Tes Tes], [-1 20], 'b--');
plot(Tes, xe1(1,1), 'bo');
plot(Tes, xe1(2,1), 'bo');
latexylabel('$$T_1, T_2$$ in K');
rotateY;
axis([-1 Th 0 1.1*max(max(X))]);
yticks([0 5 10]);
xticks([-1 0 1 2 3]);
latexxlabel('$$t$$ in s');
hold off;


%%    Umsteuerung des Ofens in kürzerer Zeit
%
figure(2)
%        Berechnung von W_S(Te)
%
Tes=5;                  % Endzeit für Zustandsüberführung
dt1=0.0001;             % für Berechnung von W ist sehr kleine Schritteweite notw.
dt2=0.005;              % Schrittweite für Simulation
Th=Tes*1.2;             % Simulationshorizont
T=[0:dt2:Th];
U=[];
%
BB=B*B';
expA=expm(A*Tes);
W=BB+expA*BB*expA';
W=W/2;
for k=1:Tes/dt1-1
  t=k*dt1;
  expAt=expm(A*t);
  W=W + expAt*BB*expAt';
end;
W=W*dt1;      % Steuerbarkeits-Gramsche
%
Wx=inv(W)*xdiff;
for k=1:Tes/dt2+1
   t=(k-1)*dt2;
   U(k,:)=-B'*expm(A*(Tes-t))'*Wx;
end;
for k=1:ceil((Th-Tes)/dt2)
   U(Tes/dt2+k+1,:)=Ukonst;
end
Ulow=min(1.1*min(U), 0);  %  Untere Grenze des Darstellungsbereichs
Uhigh=1.1*max(U);
[Y, T, X]=lsim(Ofen, U, T, x0);
%
subplot(4,1,[1,2]);
plot(T, U(:,1), 'b', [-1 0 0],[Umin Umin U(1,1)], 'b', [0 0],[Ulow, Uhigh],'k:', [Tes, Tes],[Ulow, Uhigh],'b--');
latextitle(['Ofen  $$x_0=$$',num2str(x0'),'$$\rightarrow x_e=$$',num2str(xe1')]);
xleer;
xticksempty([-1 0  2 4 6 8]);
latexylabel('$$u \;\;$$');  % in W
rotateY;
axis([-1 Th Ulow, Uhigh]);
yticks([-5000 0 5000]);
%
subplot(4,1,[3,4]);
plot(T, X, 'b');
hold on;
plot([-1 0 0],[x0(1,1), x0(1,1), X(1,1)], 'b', [-1 0 0],[x0(2,1), x0(2,1), X(1,2)], 'b');
plot([0 0], [-1 20], 'k:', [Tes Tes], [-1 20], 'b--');
plot(Tes, xe1(1,1), 'bo');
plot(Tes, xe1(2,1), 'bo');
latexylabel('$$T_1, T_2$$');   % in K
rotateY;
axis([-1 Th 0 15]);
yticks([0 5 10 15]);
xticks([-1 0  2 4 6 8]);
latexxlabel('$$t$$ in s');
hold off;



%%   Umsteuerung des Ofens bei verändertem Zielzustand
%
figure(3)
%
Umin=800;
%  statischer Zustand bei u=Umin (konstant)
x0=-inv(A)*B*Umin;
xe2=2*x0;
xdiff=expA*x0-xe2;
Wx=inv(W)*xdiff;
for k=1:Tes/dt2+1
   t=(k-1)*dt2;
   U(k,:)=-B'*expm(A*(Tes-t))'*Wx;
end;
for k=1:ceil((Th-Tes)/dt2)
   U(Tes/dt2+k+1,:)=2*Umin;
end
Ulow=min(1.1*min(U), 0);  %  Untere Grenze des Darstellungsbereichs
Uhigh=1.1*max(U);
[Y, T, X]=lsim(Ofen, U, T, x0);
%
subplot(4,1,[1,2]);
plot(T, U(:,1), 'b', [-1 0 0],[Umin Umin U(1,1)], 'b', [0 0],[Ulow, Uhigh],'k:', [Tes, Tes],[Ulow, Uhigh],'b--');
latextitle(['Ofen  $$x_0=$$',num2str(x0'),'$$\rightarrow x_e=$$',num2str(xe2')]);
xleer;
xticksempty([ 0  2 4 6 8]);
latexylabel('$$u$$ in W');
rotateY;
axis([-1 Th Ulow, Uhigh]);
yticks([0 1000 2000]);
%
subplot(4,1,[3,4]);
plot(T, X, 'b');
hold on;
plot([-1 0 0],[x0(1,1), x0(1,1), X(1,1)], 'b', [-1 0 0],[x0(2,1), x0(2,1), X(1,2)], 'b');
plot([0 0], [-1 12], 'k:', [Tes Tes], [-1 12], 'b--');
plot(Tes, xe2(1,1), 'bo');
plot(Tes, xe2(2,1), 'bo');
latexylabel('$$T_1, T_2$$ in K');
rotateY;
axis([-1 Th 0 1.1*max(max(X))]);
yticks([0 2 4 6]);
xticks([ 0  2 4 6 8]);
latexxlabel('$$t$$ in s');
hold off;


%%  Vergleich mit der Steuerung u=2*Umin  ->  xe=2*x0
%
figure(4)
U2=2*Umin*ones(size(T));
[Y2, T, X2]=lsim(Ofen, U2, T, x0);
%
subplot(4,1,[1,2]);
plot(T, U2(:,1), 'b', [-1 0 0],[Umin Umin U2(1,1)], 'b', [0 0],[Ulow, Uhigh],'k:', [Tes, Tes],[Ulow, Uhigh],'b--');
latextitle('Ofen mit verdoppelter Stellgr\"osse');
xleer;
xticksempty([ 0  2 4 6 8]);
latexylabel('$$u$$ in W');
rotateY;
axis([-1 Th Ulow, Uhigh]);
yticks([0 1000 2000]);
%
subplot(4,1,[3,4]);
plot(T, X2, 'b');
hold on;
plot([-1 0 0],[x0(1,1), x0(1,1), X2(1,1)], 'b', [-1 0 0],[x0(2,1), x0(2,1), X2(1,2)], 'b');
plot([0 0], [-1 12], 'k:', [Tes Tes], [-1 12], 'b--');
plot([3, Th], [xe2(1,1), xe2(1,1)], 'k:');
plot([3, Th], [xe2(2,1),xe2(2,1)], 'k:');
latexylabel('$$T_1, T_2$$ in K');
rotateY;
axis([-1 Th 0 1.1*max(max(X))]);
yticks([0 2 4 6]);
xticks([ 0  2 4 6 8]);
latexxlabel('$$t$$ in s');
hold off;



%%  Umsteuerung mit einem Dead-beat Regler
%
figure(5)
%x0=[0;0];
%Umin=0;
Ta=Tes/2;   % Abtastzeit
[Ad, Bd]=c2d(A, B, Ta);
%    Berechnung des Dead-beat-Reglers
M=[          Bd(1,1)                        Bd(2,1);
    Ad(1,2)*Bd(2,1)-Ad(2,2)*Bd(1,1)   Ad(2,1)*Bd(1,1)-Ad(1,1)*Bd(2,1)];
r=[Ad(1,1)+Ad(2,2);
   Ad(1,2)*Ad(2,1)-Ad(1,1)*Ad(2,2)];
K=(inv(M)*r)';
Aquer=Ad-Bd*K;
V=inv(C*inv(eye(2,2)-Aquer)*Bd);
%   Sollwert
w=5.6;
%  Ausgabegleichung für u
Cd=-K;
Dd=V;
%   Dead-beat-Stellgröße
[Ud, X]=dstep(Aquer, Bd*V, Cd, Dd, 1, 3);
%   Endzustand
%   dieser Endzustand kann durch u=Ud(2,1) erhalten werden
xe3=inv(eye(2,2)-Aquer)*Bd*V*w;
%   Kontinuierliches Stellsignal
Udb(1:Ta/dt2)=Ud(1,1);
Udb(Ta/dt2+1:2*Ta/dt2)=Ud(2,1);
Udb(2*Ta/dt2+1:Th/dt2+1)=Ud(3,1);
Udb=w*Udb;
[Ydb, T, Xdb]=lsim(Ofen, Udb, T, x0);
%   Grafische Darstellung später
%
%    Stellgröße mit der Gramschen Matrix berechnet
xdiff=expA*x0-xe3;
Wx=inv(W)*xdiff;
for k=1:Tes/dt2+1
   t=(k-1)*dt2;
   U(k,:)=-B'*expm(A*(Tes-t))'*Wx;
end;
%  Berechnung der konstanten Stellgröße, mit der der Endzustand xe
%  erhalten bleibt
xec=-inv(A)*B;
Ukonst=xe3(1,1)/xec(1,1);
for k=1:ceil((Th-Tes)/dt2)
   U(Tes/dt2+k+1,:)=Ukonst;
end
Ulow=min(1.1*min(U), 0);  %  Untere Grenze des Darstellungsbereichs
Uhigh=1.1*max(U);
[Y, T, X]=lsim(Ofen, U, T, x0);
%
subplot(4,1,[1,2]);
plot(T, U(:,1), 'b', [-1 0 0],[Umin Umin U(1,1)], 'b', [0 0],[Ulow, Uhigh],'k:', [Tes, Tes],[Ulow, Uhigh],'b--');
latextitle(['Ofen $$x_0=$$',num2str(x0'),'$$\rightarrow x_e=$$',num2str(xe3')]);
xleer;
xticksempty([0 2 4 6 8]);
latexylabel('$$u \;\;$$');    % in W
rotateY;
axis([-1 Th Ulow, Uhigh]);
yticks([0 5000]);
%
subplot(4,1,[3,4]);
plot(T, X, 'b');
hold on;
plot([-1 0 0],[x0(1,1), x0(1,1), X(1,1)], 'b', [-1 0 0],[x0(2,1), x0(2,1), X(1,2)], 'b');
plot([0 0], [-20 20], 'k:', [Tes Tes], [-20 20], 'b--');
plot(Tes, xe3(1,1), 'bo');
plot(Tes, xe3(2,1), 'bo');
latexylabel('$$T_1, T_2$$');    % in K
rotateY;
axis([-1 Th 0 1.1*max(max(X))]);
yticks([0 5 10 15]);
xticks([0 2 4 6 8]);
latexxlabel('$$t$$ in s');
hold off;


%%  Darstellung des Ergebnisses beim Dead-beat-Regler in derselben Weise
%
figure(6)
%
subplot(4,1,[1,2]);
plot(T, Udb, 'b', [-1 0 0],[Umin Umin Udb(1,1)], 'b', [0 0],[Ulow, Uhigh],'k:', [Tes, Tes],[Ulow, Uhigh],'b--');
latextitle(['Ofen $$x_0=$$',num2str(x0'),'$$\rightarrow x_e=$$',num2str(xe3')]);
xleer;
xticksempty([0 2 4 6 8]);
latexylabel('$$u$$ in W');
rotateY;
axis([-1 Th Ulow, Uhigh]);
yticks([0 5000]);
%
subplot(4,1,[3,4]);
plot(T, Xdb, 'b');
hold on;
plot([-1 0 0],[x0(1,1), x0(1,1), Xdb(1,1)], 'b', [-1 0 0],[x0(2,1), x0(2,1), Xdb(1,2)], 'b');
plot([0 0], [-20 20], 'k:', [Tes Tes], [-20 20], 'b--');
plot(Tes, xe3(1,1), 'bo');
plot(Tes, xe3(2,1), 'bo');
latexylabel('$$T_1, T_2$$ in K');
rotateY;
axis([-1 Th 0 15]);
yticks([0 5 10 15]);
xticks([0 2 4 6 8]);
latexxlabel('$$t$$ in s');
hold off;


%% Zielzustände, in denen das System verharren kann
%
figure(7)
%
subplot(4,1,[1,2]);
plot([0 15], [0 15*ksx(2,1)/ksx(1,1)],'b--');
hold on
plot(xe1(1,1), xe1(2,1),'bx');
plot(xe2(1,1), xe2(2,1),'bx');
plot(xe3(1,1), xe3(2,1),'bx');
latextitle('Zielzust\"ande, in denen der Ofen verharren kann');
latexxlabel('$$T_1$$ in K');
latexylabel('$$T_2$$ in K');
rotateY;
hold off
axis([0 15 0 7.5]);
yticks([0 2 4 6 8]);
xticks([0 5 10 15]);

%%     Bilder in .EPS-Dateien
%
Epsbildc15('ofenst');
