%   Script OFENBEOB.M
%
%   Steuerbarkeit und Beobachtbarkeit einer Beschichtungsanlage
%
% J. Lunze
% 15.4.2003
% für 2. Auflage: 29.12.2006
% für 3. Auflage: 4.12.2011
% für 4. Auflage: 10.11.2015
% für 5. Auflage: 10.01.2020
%
%
clear
close all
bspofen;    %  Definition des Ofens (A, B, C, D);

%%    Berechnung des Anfangszustands mit zwei Methoden
%
figure(1);
Te=3;       % Beobachtungshorizont
dt1=0.0001;   % Schrittweite zur Berechnung von W_B
dt2=0.005;    % Schrittweite für Simulation
T=[0:dt2:Te];
kend=length(T);
Umin=800;
U=Umin+1600*[0:dt2:Te/2];
Umax=max(U);
U=[U Umax*ones(1,Te/2/dt2)];
%  statischer Zustand bei u=Umin (konstant)
x0=-inv(A)*B*Umin
[Y, T, X]=lsim(Ofen, U, T, x0);
%
subplot(4,1,[1,2]);
plot(T, U, 'b', [-10 0],[Umin Umin], 'b', [0 0],[-Umax 1.1*Umax],'k:', [-Te/10 Te],[0 0],'k:');
latexylabel('$$u$$ in W');
rotateY;
latextitle('Temperatur von Werkst\"uck und Halter');
xleer;
xticksempty([0 1 2 3]);
axis([-Te/10 Te -Umax/10 1.1*Umax]);
yticks([0 1000 2000 3000]);
%
subplot(4,1,[3,4]);
plot(T, X, 'b', [-10 0],[X(1,1) X(1,1)],'b', [-10 0],[X(1,2) X(1,2)],'b', [0 0],[-5 20],'k:');
latexylabel('$$T_1, T_2$$ in K');
rotateY;
latexxlabel('$$t$$ in s');
axis([-Te/10 Te 0 1.1*max(max(X))]);
yticks([0 2 4 6 8]);
xticks([0 1 2 3]);
latextext(1.5, 7, '$$T_1$$');
latextext(1.5, 3, '$$T_2$$');



%%  Berechnung des Anfangszustandes aus zwei Messpunkten
%
figure(2);
subplot(4,1,[1,2]);
[Yerzw, T, Xerzw]=lsim(Ofen, U, T, [0;0]);
Yfrei=Y-Yerzw;
plot(T, Yfrei, 'b');
hold on;
latexylabel('$$y_{\rm frei}$$ in K');
rotateY;
latexxlabel('$$t$$ in s');
latextitle('Eigenbewegung des Ofens');
axis([0 Te 0.9*min(Yfrei) 1.1*max(Yfrei)]);
yticks([0 0.5 1]);
xticks([0 1 2 3]);
%
%     zwei Meßpunkte herausgreifen bei t1=0,1, t2=1.0
t1=0.1;
t2=1;
k1=ceil(t1/dt2+1);
[t1, Yfrei(k1, 1)]
plot(t1, Yfrei(k1,1), 'bo');
k2=ceil(t2/dt2+1);
[t2, Yfrei(k2,1)]
plot(t2, Yfrei(k2,1), 'bo');
plot([t1 t1], [0 Yfrei(k1,1)], 'k:', [0 t1], [Yfrei(k1,1) Yfrei(k1,1)], 'k:');
plot([t2 t2], [0 Yfrei(k2,1)], 'k:', [0 t2], [Yfrei(k2,1) Yfrei(k2,1)], 'k:');
hold off;
M=[C*expm(A*k1*dt2); C*expm(A*k2*dt2)];
x0hat1=inv(M)*[Yfrei(k1,1); Yfrei(k2,1)]


%%   Berechnung des Anfangszustandes
%        Berechnung von W_B(Te)
%
CC=C'*C;
expA=expm(A*Te);
WB=CC+expA'*CC*expA;
WB=WB/2;
for k=1:Te/dt1-1
  t=k*dt1;
  expAt=expm(A*t);
  WB=WB + expAt'*CC*expAt;
end;
WB=WB*dt1;
%   Berechnung von int_0^Te exp(A't) C' y_frei(tau) dtau
kmax=Te/dt2+1;
Int = C'*Yfrei(1,1)/2 + expA'*C'*Yfrei(kmax,1)/2;
for k=1:kmax
  t=k*dt2;
  expAt=expm(A*t);
  Int=Int + expAt'*C'*Yfrei(k,1);
end;
Int=Int*dt2;
x0hat2=inv(WB)*Int

%%  Beobachterentwurf
%
figure(3)
LambdaB=5*[-1, -1.5];
LT=acker(A', C', LambdaB);
L=LT';
Ab=[A zeros(n,n); L*C A-L*C];
Bb=[B;B];
Cb=[0 1  0 -1];       % Differenz y-y^
Db=[0];
%  Zustände: x1, x2  Temperaturen im Ofen
%            x3, x4  beobachtete Temperaturen
SystemBeobachter=ss(Ab, Bb, Cb, Db);
%  Ausgangsgröße bei x0:
y0=C*x0;
%  Beobachterzustand aus Annahme, dass alle Zustände auf diesem Wert liegen
x0b1=[x0; y0; y0];
%  Konstante Eingangsgröße
Uc=Umin*ones(size(T));
[Yd1, T, X1]=lsim(SystemBeobachter, Uc, T, x0b1);
%
subplot(6,1,[1,2]);
%
plot(T, Uc, 'b',[-Te/10 Te],[0 0],'k:');
latexylabel('$$u \;\;$$');  % in W
rotateY;
latextitle('Beobachtung der Werkst\"ucktemperatur');
xleer;
axis([0 1 -Umax/10 1.1*Umax]);
yticks([0 1000 2000 3000]);
%
subplot(6,1,[3,4]);
plot(T, X1(:,1),'b--', T, X1(:,3), 'b', T, X1(:,2),'b--', T, X1(:,4), 'b');
hold on
plot([0 Te], [0 0], 'k:');
xleer
latexylabel('$$T_1, T_2$$');   % in K
rotateY;
axis([0 1 0 2.5]); 
yticks([0 1 2]);
latextext(0.8, 1.5, '$$T_1$$');
latextext(0.8, 0.5, '$$T_2$$');
hold off
%
subplot(6,1,[5,6]);
plot(T, Yd1, 'b');
hold on
plot([-1 3],[0 0],'k:');
latexylabel('$$y-\hat{y}$$');  % in K
rotateY;
latexxlabel('$$t$$ in s');
axis([0 1 -0.5 1])
yticks([-0.5 0 0.5 1]);
xticks([0 0.2 0.4 0.6 0.8 1]);

%
%  für zweiten Anfangszustand des Beobachters
%  Beobachterzustand bei null
figure(4)
x0b2=[x0; 0; 0];
[Yd2, T, X2]=lsim(SystemBeobachter, Uc, T, x0b2);
%
subplot(6,1,[1,2]);
plot(T, Uc, 'b', [-Te/10 Te],[0 0],'k:');
latexylabel('$$u \;\;$$');
rotateY;
latextitle('Beobachtung der Werkst\"ucktemperatur');
xleer;
axis([0 1 -Umax/10 1.1*Umax]);
yticks([0 1000 2000 3000]);
%
subplot(6,1,[3,4]);
plot(T, X2(:,1),'b--', T, X2(:,3), 'b',T, X2(:,2),'b--', T, X2(:,4), 'b');
hold on
plot([0 Te], [0 0], 'k:');
xleer
latexylabel('$$T_1, T_2$$');
rotateY;
axis([0 1 0 2.5]);
yticks([0 1 2]);
latextext(0.8, 1.5, '$$T_1$$');
latextext(0.8, 0.5, '$$T_2$$');
hold off
%
subplot(6,1,[5,6]);
plot(T, Yd2, 'b');
hold on
plot([-1 3],[0 0],'k:');
latexylabel('$$y-\hat{y}$$');
rotateY;
latexxlabel('$$t$$ in s');
axis([0 1 -0.5 1])
yticks([-0.5 0 0.5 1]);
xticks([0 0.2 0.4 0.6 0.8 1]);

%%
figure(5)
%  Eingangsgröße wie im ersten Bild
[Yd3, T, X3]=lsim(SystemBeobachter, U, T, x0b2);
%
subplot(6,1,[1,2]);
plot(T, U, 'b', [-10 0],[Umin Umin],'b', [0 0],[-Umax 1.1*Umax],'k:', [-Te/10 Te],[0 0],'k:');
latexylabel('$$u$$ in W');
rotateY;
latextitle('Beobachtung der Werkst\"ucktemperatur');
xleer;
xticksempty([0 1 2 3]);
axis([-Te/10 Te -Umax/10 1.1*Umax]);
yticks([0 1000 2000 3000]);
%
subplot(6,1,[3,4]);
plot(T, X3(:,3),  'b', T, X3(:,4), 'b', [0 Te], [0 0], 'k:');
hold on
plot(T, X, 'b--',[-10 0],[X(1,1) X(1,1)],'b--',[-10 0],[X(1,2) X(1,2)], 'b--',[0 0],[-5 20],'k:');
latexylabel('$$T_1, T_2$$ in K');
rotateY;
axis([-Te/10 Te 0 1.1*max(max(X))]);
yticks([0 2 4 6 8]);
xleer
xticksempty([0 1 2 3]);
latextext(1.2, 7, '$$T_1$$');
latextext(1.5, 4, '$$T_2$$');
%
%
subplot(6,1,[5,6]);
plot(T, Yd3, 'b', [-10 0 0], [0 0 Yd3(1,1)], 'b');
hold on
plot([-1 Te],[0 0],'k:');
latexylabel('$$y-\hat{y}$$ in K');
rotateY;
latexxlabel('$$t$$ in s');
axis([-Te/10 Te -0.5 1])
yticks([-0.5 0 0.5 1]);
xticks([0 1 2 3]);
hold off

%% Beobachter mit Impulsstörung
%
figure(6)
%  Beobachter
AB=A-L*C;
BU=B;   %  Eingang: Wärmestrom u
BY=L;   %  Eingang: gemessene Temperatur des Werkstückhalters
CB=-C;
DB=[0 1];  % Ausgang: y - y^
Beobachter=ss(AB, [BU, BY], CB, DB);
%  Eingangsgröße fällt kurzzeitig zwischen Zeitpunkten T10 und T20 aus
T10=1.5;
T20=1.65;
k10=T10/dt2;
k20=T20/dt2;
Uimpuls=U;
for k=k10:k20
    Uimpuls(1,k)=0;
end
%  Ofen bei veränderter Eingangsgröße
[Yimpuls, T, Ximpuls]=lsim(Ofen, Uimpuls, T, x0);
[YBimpuls, T, XBimpuls]=lsim(Beobachter, [U', Yimpuls], T);
%
subplot(6,1,[1,2]);
plot(T, Uimpuls,'b--');
hold on
plot(T, U, 'b', [-10 0], [U(1,1), U(1,1)], 'b');
plot([0 0], [-1000 5000],'k:');
xleer
xticksempty([0 1 2 3]);
latexylabel('$$u \;\;$$');
rotateY;
axis([-Te/10 Te -Umax/10 1.1*Umax]);
yticks([0 1000 2000 3000]);
latextitle('Beobachtung der Werkst\"ucktemperatur mit impulsf\"ormiger St\"orung');
%
subplot(6,1,[3,4]);
plot(T, Ximpuls, 'b--',[-10 0],[Ximpuls(1,1) Ximpuls(1,1)],'b--',[-10 0],[Ximpuls(1,2) Ximpuls(1,2)], 'b--',[0 0],[-5 20],'k:');
hold on
plot(T, XBimpuls(:,1), 'b', T, XBimpuls(:,2), 'b', [0 Te], [0 0], 'k:');
xleer
xticksempty([0 1 2 3]);
latexylabel('$$T_1, T_2$$');
rotateY;
axis([-Te/10 Te 0 1.1*max(max(X))]);
yticks([0 2 4 6 8]);
latextext(1.2, 7, '$$T_1$$');
latextext(1.2, 3, '$$T_2$$');
hold off
%
subplot(6,1,[5,6]);
plot(T, YBimpuls, 'b', [-10 0 0], [0 0 YBimpuls(1,1)], 'b', [-Te/10 Te], [0 0],'k:')
latexxlabel('$$t$$ in s');
latexylabel('$$y-\hat{y}$$');
rotateY;
axis([-Te/10 Te -0.5 1]);
yticks([-0.5 0 0.5 1]);
xticks([0 1 2 3]);



%%  Beobachter mit sprungförmiger Störung
%
figure(7)
%  Eingangsgröße hat systematischen Fehler Udelta nach T30=1
T30=1;
k30=T30/dt2;
Udelta=0.25*max(U);
Usprung=U;
for k=k30:kend
    Usprung(1,k)=Usprung(1,k)-Udelta;
end
%  Ofen bei veränderter Eingangsgröße
[Ysprung, T, Xsprung]=lsim(Ofen, Usprung, T, x0);
[YBsprung, T, XBsprung]=lsim(Beobachter, [U', Ysprung], T);
%
subplot(6,1,[1,2]);
plot(T,Usprung, 'b--', [-10 0],[Usprung(1,1), Usprung(1,1)], 'b--', T, U, 'b')
hold on
plot([0 0], [-1000 5000],'k:');
xleer
xticksempty([0 1 2 3]);
latexylabel('$$u \;\;$$');
rotateY;
axis([-Te/10 Te -Umax/10 1.1*Umax]);
yticks([0 1000 2000 3000]);
latextitle('Beobachtung der Werkst\"ucktemperatur mit sprungf\"ormiger St\"orung');
%
subplot(6,1,[3,4]);
plot(T, Xsprung, 'b--',[-10 0],[Xsprung(1,1) Xsprung(1,1)],'b--',...
    [-10 0],[Xsprung(1,2) Xsprung(1,2)], 'b--',[0 0],[-5 20],'k:');
hold on
plot(T, XBsprung(:,1), 'b', T, XBsprung(:,2), 'b', [0 Te], [0 0], 'k:');
xleer
xticksempty([0 1 2 3]);
latexylabel('$$T_1, T_2$$');
rotateY;
axis([-Te/10 Te 0 1.1*max(max(X))]);
yticks([0 2 4 6 8]);
hold off
%
subplot(6,1,[5,6]);
plot(T, YBsprung, 'b', [-10 0 0], [0 0 YBsprung(1,1)], 'b', [-Te/10 Te], [0 0],'k:')
latexxlabel('$$t$$ in s');
latexylabel('$$y-\hat{y}$$');
rotateY;
axis([-Te/10 Te -0.5 1]);
yticks([-0.5 0 0.5 1]);
xticks([0 1 2 3]);



%%  Beobachter mit stochastischer Störung
%
figure(8)
Ys=Y+randn(size(Y));
[YB, TB, XB]=lsim(Beobachter, [U', Ys], T);
%
subplot(4,1,[1,2]);
plot(T, X, 'b--',[-10 0],[X(1,1) X(1,1)],'b--',[-10 0],[X(1,2) X(1,2)], 'b--',[0 0],[-5 20],'k:');
hold on
plot(T, XB(:,1), 'b', T, XB(:,2), 'b', [0 Te], [0 0], 'k:');
xleer
xticksempty([0 1 2 3]);
latexylabel('$$T_1, T_2$$');
rotateY;
axis([-Te/10 Te 0 1.1*max(max(X))]);
yticks([0 2 4 6 8]);
latextitle('Beobachtung der Werkst\"ucktemperatur mit stochastischer St\"orung');
latextext(1.2, 7, '$$T_1$$');
latextext(1.2, 3, '$$T_2$$');
hold off
%
subplot(4,1,[3,4]);
plot(T, Ys, 'b')
latexxlabel('$$t$$ in s');
latexylabel('$$y \;\;$$');
rotateY;
axis([-Te/10 Te 0 1.1*max(max(X))]);
yticks([0 2 4 6 8]);
xticks([0 1 2 3]);



%%  Schnellerer Beobachter
%
figure(9)
LambdaB2=10*[-1, -1.5];
LT2=acker(A', C', LambdaB2);
L2=LT2';
Ab2=[A zeros(n,n); L2*C A-L2*C];
Bb=[B;B];
Cb=[0 1  0 -1];       % Differenz y-y^
Db=[0];
SystemBeobachter2=ss(Ab2, Bb, Cb, Db);
[Yd4, T, X4]=lsim(SystemBeobachter2, Uc, T, x0b1);
%
% subplot(4,1,[3,4]);
% plot(T, X1(:,3), 'b-.', T, X4(:,1),'b--', T, X4(:,3),  [0 Te], [0 0], 'k:');
% latexxlabel('$$t$$ in s');
% latexylabel('$$T_1$$ in K');
% rotateY;
% axis([0 1 0 2.5]);
%yticks([]);
%xticks([]);
%  Beobachter 2
AB2=A-L2*C;
BU=B;    %  Eingang: Wärmestrom u
BY=L2;   %  Eingang: gemessene Temperatur des Werkstückhalters
CB=-C;
DB=[0 1];  % Ausgang: y - y^
Beobachter2=ss(AB2, [BU, BY], CB, DB);
[YB2, T, XB2]=lsim(Beobachter2, [U', Ys], T);
%
subplot(4,1,[1,2]);
plot(T, X, 'b--',[-10 0],[X(1,1) X(1,1)],'b--',[-10 0],[X(1,2) X(1,2)], 'b--',[0 0],[-5 20],'k:');
hold on
plot(T, XB2(:,1), 'b', T, XB2(:,2), 'b', [0 Te], [0 0], 'k:');
latextitle('Schnellerer Beobachter der Werkst\"ucktemperatur');
latexylabel('$$T_1, T_2$$');
xleer
xticksempty([0 1 2 3]);
rotateY;
axis([-Te/10 Te 0 1.1*max(max(X))]);
yticks([0 2 4 6 8]);

%
subplot(4,1,[3,4]);
plot(T, Ys, 'b')
latexxlabel('$$t$$ in s');
latexylabel('$$y \;\;$$');
rotateY;
axis([-Te/10 Te 0 1.1*max(max(X))]);
yticks([0 2 4 6 8]);
xticks([0 1 2 3]);



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