%   Script PERFREGTEMP.M
%
%   Bestimmung einer Solltrajektorie für einen Rührkessel
%     Eingang: Heizleistung
%     Ausgang: Temperatur der Flüssigkeit
%     (wie in Beispiel 12.1)
%
% J. Lunze
% für 5. Auflage: 6.6.2005
% für 6. Auflage: 4.12.2006
% für 7. Auflage: 28.4.2008
% für 8. Auflage: 24.1.2010
% für 9. Auflage: 13.2.2012
% für 10. Auflage: 27.1.2014
% für 12. Auflage: 23.03.2020
%
%   
clear
close all
Bspreak;    %  liefert "Reaktor" = Temperaturverhalten des Reaktors


%%  Umsteuerung Temperatur von T=0 nach T_e=30 Grad
%
figure(1)
Te=8;       % Endzeit in Minuten
we=30;      %  Endwert der Temperatur
%    Beschreibung der Trajektorie durch ein Polynom 3. Ordnung mit aw0=0
%    Berechnung der Polynomkoeffizienten aw=[aw3; aw2; aw1] als Lösung der
%    Gl. M aw= N
%    Forderung: w(0)=0; wPunkt(0)=0, w(Te)=30, WPunkt(Te)=0
aw0=0;
aw1=0;
M=[Te^3 Te^2; 3*Te^2 2*Te];
N=[30; 0];
aVektor=inv(M)*N;
aw3=aVektor(1,1);
aw2=aVektor(2,1);
WPoly=[aw3, aw2, aw1, aw0];   % Vektor der Polynomkoeffizienten von w(t)
T=[0:0.01:Te];
W=polyval(WPoly, T);
WPunktPoly=[3*aw3, 2*aw2, aw1];
WPunkt=polyval(WPunktPoly, T);
WZweipunktPoly=[6*aw3, 2*aw2];
WZweipunkt=polyval(WZweipunktPoly, T);
Tbeg=-1;
%
subplot(8,1,[1,2]);
plot(T, W, 'b')
hold on
plot(0, 0, 'bo');
plot(Te, we, 'bo');
plot([Tbeg 0 0],[0 0 W(1,1)], 'b', [Te, Te+1], [we, we], 'b');
axis([Tbeg, Te+1 -2 33]);
xleer;
yticks([0, 30]);
latexylabel('$$w$$ in $$^o$$C');
rotateY;
hold off
%
subplot(8,1,[3,4]);
plot(T, WPunkt, 'b')
hold on
plot([Tbeg 0 0],[0 0 WPunkt(1,1)], 'b', [Te, Te+1], [0 0], 'b');
latexylabel('$$\frac{\mathrm{d}w}{\mathrm{d}t}$$ in $$\frac{\rm K}{\rm s}$$');
yticks([0, 15]);
axis([Tbeg, Te+1 -1 18]);
xleer;
rotateY;
hold off
%
subplot(8,1,[5,6]);
plot(T, WZweipunkt, 'b');
hold on
LengthWZweipunkt=length(WZweipunkt);
plot([Tbeg 0 0],[0 0 WZweipunkt(1,1)], 'b',...
     [Te Te Te+1], [WZweipunkt(LengthWZweipunkt),0, 0], 'b');
latexylabel('$$\frac{\mathrm{d}^2w}{\mathrm{d}t^2}$$ in $$\frac{\rm K}{\rm s^2}\;\;$$');
yticks([-20, 0, 20]);
axis([Tbeg, Te+1 -23 23]);
rotateY;
xleer;
hold off
%
%  Berechnung der Stellgröße für diese Solltrajektorie
%  Koeffizienten des charakteristischen Polynoms des Reaktors
[z, n]=tfdata(Reaktor, 'v');
U=n(1)*WZweipunkt+n(2)*WPunkt+n(3)*W;
U=U/z(3);
ks=dcgain(Reaktor);
%
subplot(8,1,[7,8]);
plot(T, U, 'b', [0 Te+1],[25 25],'k:');
hold on
LengthU=length(U);
plot([Tbeg 0 0], [0 0 U(1,1)], 'b', [Te Te Te+1], [U(LengthU), we/ks, we/ks], 'b');
axis([Tbeg, Te+1, -3 30]);
yticks([0 30]);
xticks([0 2 4 6 8 10]);
latexylabel('$$u$$ in kW');
rotateY;
hold off
latexxlabel('$$t$$ in min');
 

%%   Solltrajektorie noch einmal größer darstellen
%
figure(2)
subplot(2,1,1)
plot(T, W, 'b');
hold on
plot([Tbeg 0 0],[0 0 W(1,1)], 'b',[Te, Te+1], [we, we], 'b');
plot(0, 0, 'bo');
plot(Te, we, 'bo');
latexylabel('$$w$$ in $$^o$$C');
rotateY;
latexxlabel('$$t$$ in min');
axis([Tbeg, Te+1 -2 33]);
xticks([0 2 4 6 8 10]);
yticks([0 15 30]);
hold off
 
%%
%
figure(3)
%  Umsteuerung Temperatur von T=0 nach T_e=30 Grad
Te=3;       % Endzeit in Minuten
%    Beschreibung der Trajektorie durch ein Polynom 3. Ordnung mit aw0=0
%    Berechnung der Polynomkoeffizienten aw=[aw3; aw2; aw1] als Lösung der
%    Gl. M aw= N
%    Forderung: w(0)=0; wPunkt(0)=0, w(Te)=30, WPunkt(Te)=0
aw0=0;
aw1=0;
M=[Te^3 Te^2; 3*Te^2 2*Te];
N=[we; 0];
aVektor=inv(M)*N;
aw3=aVektor(1,1);
aw2=aVektor(2,1);
WPoly=[aw3, aw2, aw1, aw0];   % Vektor der Polynomkoeffizienten von w(t)
T=[0:0.01:Te];
W=polyval(WPoly, T);
WPunktPoly=[3*aw3, 2*aw2, aw1];
WPunkt=polyval(WPunktPoly, T);
WZweipunktPoly=[6*aw3, 2*aw2];
WZweipunkt=polyval(WZweipunktPoly, T);
WDreipunktPoly=[6*aw3];
WDreipunkt=polyval(WDreipunktPoly, T);
WVierpunkt=zeros(size(WDreipunkt));
Tend=9;   % Achsenskalierung wie bisher
%
subplot(8,1,[1,2]);
plot(T, W, 'b')
hold on
plot(0, 0, 'bo');
plot(Te, we, 'bo');
plot([Tbeg 0 0],[0 0 W(1,1)], 'b',[Te, Tend], [we, we], 'b');
axis([Tbeg, Tend -2 33]);
xleer;
yticks([0 30]);
latexylabel('$$w$$ in $$^o$$C');
rotateY;
hold off
%
subplot(8,1,[3,4]);
plot(T, WPunkt, 'b')
hold on
plot([Tbeg 0 0],[0 0 WPunkt(1,1)], 'b', [Te, Tend], [0 0], 'b');
latexylabel('$$\frac{\mathrm{d}w}{\mathrm{d}t}$$ in $$\frac{\rm K}{\rm s}\;\;$$');
yticks([0 15]);
axis([Tbeg, Tend -1 18]);
xleer;
rotateY;
hold off
%
subplot(8,1,[5,6]);
plot(T, WZweipunkt, 'b');
hold on
LengthWZweipunkt=length(WZweipunkt);
plot([Tbeg 0 0],[0 0 WZweipunkt(1,1)], 'b',...
    [Te Te Tend], [WZweipunkt(LengthWZweipunkt),0, 0], 'b');
latexylabel('$$\frac{\mathrm{d}^2w}{\mathrm{d}t^2}$$ in $$\frac{\rm K}{\rm s^2}\;\;$$');
yticks([-20 0 20]);
axis([Tbeg, Tend -23 23]);
rotateY;
xleer;
hold off
%
%  Berechnung der Stellgröße für diese Solltrajektorie
%  Koeffizienten des charakteristischen Polynoms des Reaktors
U=n(1)*WZweipunkt+n(2)*WPunkt+n(3)*W;
U=U/z(3);
ks=dcgain(Reaktor);
subplot(8,1,[7,8]);
plot(T, U, 'b', [0 Tend],[25 25],'k:');
hold on
LengthU=length(U);
plot([Tbeg 0 0], [0 0 U(1,1)], 'b', [Te Te Tend], [U(LengthU), we/ks, we/ks], 'b');
axis([Tbeg, Tend, -3 30]);
latexylabel('$$u$$ in kW');
rotateY;
yticks([0 30]);
xticks([0 2 4 6 8 10]);
latexxlabel('$$t$$ in min');
rotateY;
hold off

%%
%   Berechnung der E/A-Normalform
%
[A, b, c, d]=ssdata(Reaktor);
TE=inv([c; c*A]);
AE=inv(TE)*A*TE
bE=inv(TE)*b
cE=c*TE
d=0;




%%
%
Epsbildc15('perfregtemp');
