%   Script PerfregSchwingkreis.M
%
%   Bestimmung einer Solltrajektorie für einen Reihenschwingkreis
%     Eingangsspannung u_1
%     Ausgangsspannung u_2
%
% J. Lunze
% 25.1.2010
%
% für 8. Auflage: 26.1.2010
% für 9. Auflage: 10.2.2012
% für 10. Auflage: 27.1.2014
% für 12. Auflage: 24.07.2019
%
%   
clear
close all;
BspSchwingkreis;           % liefert A, b, c, d
                           % RLC=ss(A, b, c, d);
                           % x0=[1; 1];

%  Transformation in E/A-Normalform
r=1;
n=2;
SS=ctrb(A, b);
sRT=[0 1]*inv(SS);
TE=inv([c; sRT]);
AE=inv(TE)*A*TE;
bE=inv(TE)*b;
cE=c*TE;
RLCinEAForm=ss(AE, bE, cE, d);

%   Nulldynamik
ANN=AE(r+1:n, r+1:n);
bq=c*b;
bN=AE(r+1:n,1);            % Eingang ist y(t)=w(t)
eNT=AE(r,r+1:n);
cN=-1/bq*eNT;              % Ausgabe ist uVN= -1/bq eNT*xN = uVI
dN=0;
Nulldynamik=ss(ANN, bN, cN, dN);
eIT=AE(r, 1:r);
vT=[eIT -1];


%%   Umsteuerung durch konstante Eingangsgröße
%
figure(1);
ks=dcgain(RLC);
we=1;    % Sollposition nach Arbeitspunktwechsel
ue=we/ks;
dt=0.0002;
Tend=0.04;
T=0:dt:Tend;
Tbeg=-0.005;
subplot(6,1,[1,2]);
plot([Tbeg 0 0 Tend], [0 0 ue ue], 'b');
xleer;
yticks([0 1]);
axis([Tbeg, Tend, -0.1, 1.1]);
latextitle('Verhalten des Reihenschwingkreises bei konstanter Eingangsgr\"osse');
latexylabel('$$u$$ in V');
rotateY;
%
subplot(6,1,[3,4]);
Y1=step(RLC,T)*ue;
plot(T, Y1, 'b');
hold on
plot([Tbeg 0 0],[0 0 Y1(1,1)], 'b');
plot([Tbeg 0 0 Tend],[0 0 1 1],'b--'); % Sollwert
latexylabel('$$y$$ in V');
rotateY;
axis([Tbeg Tend -0.3 1.7]);
xticks([0 0.01 0.02 0.03 0.04]);
yticks([0 1]);
latexxlabel('$$t$$ in s');
hold off;
 

%%   Berechnung einer Solltrajektorie
%
figure(2)
%  Umsteuerung der Ausgangsspannung von y0=0 in ye=1
Te=0.01;       % Endzeit in Sekunden
we=1;          %  Endwert der Spannung
%    Beschreibung der Trajektorie durch ein Polynom 3. Ordnung mit aw0=0
%    Berechnung der Polynomkoeffizienten
%    aw=[aw3; aw2; aw1; aw0] als Lösung der
%    Gl. m M aw= N
%    Forderung: w(0)=0; wPunkt(0)=0,
%               w(Te)=we, 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:dt:Te];
W=polyval(WPoly, T);                   % w(t)
WPunktPoly=[3*aw3, 2*aw2, aw1];
WPunkt=polyval(WPunktPoly, T);         % wdot(t)
%
subplot(6,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, Tend,  -0.2 1.2]);
xleer;
yticks([0, 1]);
latexylabel('$$w$$ in V');
rotateY;
hold off
%
subplot(6,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 V}{\rm s}$$');
axis([Tbeg, Tend -10 200]);
latexxlabel('$$t$$ in s');
yticks([0, 200]);
xticks([0 0.01 0.02 0.03 0.04]);
rotateY;
hold off


%%  Berechnung der Stellgröße für diese Solltrajektorie
%   u(t) = -1/bq eNT*xN - 1/bq vT x_ISoll
%        = uVN + uVI
%
figure(3)
%  Zeitintervall  0....Te
uVI=-1/bq*(vT(1,1)*W' + vT(1,2)*WPunkt'); % uVI als Spaltenvektor
[uVN, T, XN]=lsim(Nulldynamik, W, T);  % Bestimmung von uVN mit XN0=0
uV=uVN+uVI;
xNe=XN(length(XN),:);    % xN(te) -  Zustand der Nulldynamik am
                         % Ende der Umsteuerung
subplot(6,1,[1,2]);
plot(T, uV, 'b');
hold on
%plot(T, uVI, 'b--');      % Anteil uVI durch Integratordynamik
plot([Tbeg 0 0], [0 0 uV(1,1)], 'b');
% Zeitintervall  Te....Tend
% hier ist y(t) = we
T2=Te:dt:Tend;
Ye=ones(size(T2))*we;                    % y ist konstant
uVI2=-1/bq*eIT(1,1)*we*ones(size(T2))';  % konstanter Anteil an u(t)
u2=lsim(Nulldynamik,Ye,T2,xNe);
u2=uVI2+u2;
plot(T2, u2, 'b');
plot([Te, Te],[uV(length(uV),1), u2(1,1)], 'b');
plot([Te Te Tend],[uV(length(uV),1) we/ks, we/ks],'b--');      %  u=konst. im Arbeitspunkt
axis([Tbeg, Tend, -0.1, 1.2]);
xleer;
yticks([0 1]);
latexylabel('$$u$$ in V');
rotateY;
hold off
%
subplot(6,1,[3,4]);
%  Zeitintervall 0...Te
[Y, T, X]=lsim(RLC,uV,T);
xe=X(length(X),:);    % x(te) -  Zustand der Integratordynamik am
                         % Ende der Umsteuerung
plot(T, Y, 'b', [-5 0 0],[0 0 Y(1,1)], 'b');
hold on
plot(0, 0, 'bo');
plot(Te, we, 'bo');
% Zeitintervall Te....Tend
Y2=lsim(RLC,u2,T2,xe);
plot(T2, Y2, 'b');
latexylabel('$$y$$ in V');
rotateY;
%  im Vergleich dazu das Verhalten bei kosntanter Eingangsröße
u2konst=we/ks*ones(1,length(T2));
Y2k=lsim(RLC,u2konst,T2,xe);
plot(T2,Y2k,'b--');
plot(T2, u2konst,'b--');   % Verhalten bei konst. Eingang
axis([Tbeg Tend -0.3 1.7]);
yticks([0, 1]);
xticks([0 0.01 0.02 0.03 0.04]);
latexxlabel('$$t$$ in s');


%%
%
Epsbildc15('PerfregRSK');
