%   Script WOKENTW.M
%
%   Wurzelortskurve und Entwurf mit WOK
%
% J. Lunze
% 24.1.1996
% für 3. Auflage: 10.2.2001
% für 4. Auflage: 18.9.2003
% für 5. Auflage: 18.4.2005
% für 7. Auflage: 28.4.2008
% für 8. Auflage: 1.1.2010
% für 9. Auflage: 16.2.2012
% für 10. Auflage: 27.1.2014
% für 12. Auflage: 5.08.2019
%
%   
clear
close all

%% Beispiel: Invertiertes Pendel
%
figure(1)
%   
g = 9.81;   % m/s^2
m = 0.3;    % kg
l = 0.5;    % m
            % Zeit in s
A1 = [
           0        1
     3*g/2/l  0];
b1 = [0
     3/m/l/l];
c1=[1 0];             % Winkel phi
d1=[0];
%
subplot(3,3,[1, 2, 4, 5]);
Pendel=ss(A1, b1, c1, d1);
rlocprintblue(Pendel);
latextitle('Wurzelortskurve des P-geregelten Pendels');
axis([-9 9 -9 9]);
xticks([-5 0 5]);
yticks([-5 0 5]);
latexxlabel('Re');
latexylabel('Im');
latextext(-6, 1, '$$s_1$$');
latextext(5, 1, '$$s_2$$');
rotateY;



%% PD-Regler k_P \frac{T_D s + 1}{T_1 s + 1},
%  der den stabilen Pol vollständig kompensiert
%  d.h. Nullstelle des PD-Reglers kompensiert den Pol der Strecke
%
figure(2);
subplot(3,3,[1, 2, 4, 5]);
T1=0.125;               % gewählter neuer Pol der offenen Kette
lambda=eig(A1);
TD=1/abs(lambda(1));
zPD = [TD, 1];
nPD = [T1 1];
PDRegler=tf(zPD, nPD);
OffeneKette = series(PDRegler, Pendel);
OffeneKettemin=minreal(OffeneKette);
rlocprintblue(OffeneKettemin);
hold on;
latextitle('Wurzelortskurve des PD-geregelten Systems');
kP=0.736;                     % kritische Verstärkung
                              % hier erhält man k_P, nicht k, weil
                              % die charakteristische Gleichung nicht
                              % in Standardform überführt wurde
[r1, kP]=rlocus(OffeneKettemin, kP);
plot(real(r1(1)), imag(r1(1)), 'bs');
plot(real(r1(2)), imag(r1(2)), 'bs');
axis([-9 9 -9 9]);
xticks([-5 0 5]);
yticks([-5 0 5]);
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
hold off;



%% PD-Regler, der den stabilen Pol vollständig kompensiert
%  und Forderungen an das Übergangsverhalten erfüllt
%
figure(3)
subplot(3,3,[1, 2, 4, 5]);
T1=1/9.5;               % gewählter neuer Pol
nPD = [T1 1];
PDRegler=tf(zPD, nPD);
OffeneKette = series(PDRegler, Pendel);
OffeneKettemin=minreal(OffeneKette);
rlocprintblue(OffeneKettemin);
hold on;
latextitle('Wurzelortskurve des PD-geregelten Systems');
phiimag=[15 0 -15];
phireal=[-15 0 -15];
plot(phireal, phiimag, 'k:');
kP=0.8535;               %Verstärkung zur Erfüllung der Güteforderungen
[r1, kP]=rlocus(OffeneKettemin, kP);
plot(real(r1(1)), imag(r1(1)), 'bs');
plot(real(r1(2)), imag(r1(2)), 'bs');
axis([-10 6 -8 8]);
xticks([-10 -5 0 5]);
yticks([-5 0 5]);
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
hold off;



%% PD-Regler, der den stabilen Pol nicht kompensiert
%  d.h.  20%''Fehler'' von TD
%
figure(4);
subplot(3,3,[1, 2, 4, 5]);
TD1=0.8*TD;
zPD1 = [TD1, 1];
PDRegler1=tf(zPD1, nPD);
OffeneKette1= series(PDRegler1, Pendel);
rlocprintblue(OffeneKette1);
hold on;
latextitle('Wurzelortskurve des PD-geregelten Systems');
phiimag=[15 0 -15];
phireal=[-15 0 -15];
plot(phireal, phiimag, 'k:');
axis([-10 6 -8 8]);
xticks([-10 -5 0 5]);
yticks([-5 0 5]);
axis('square');
[r, kP]=rlocus(OffeneKette1, kP);
plot(real(r(1)), imag(r(1)), 'bs');
plot(real(r(2)), imag(r(2)), 'bs');
plot(real(r(3)), imag(r(3)), 'bs');
latexxlabel('Re');
latexylabel('Im');
rotateY;
hold off;



%%  Übergangsverhalten des geregelten Systems
%   
figure(5);
%                   Zustandsraummodell des Reglers
Ar=-1/T1;
br=1/T1;
cr=-kP*(1-TD/T1);
dr=-kP*TD/T1;
%                   Zustandsraummodell des geregelten Pendels
Ag=[A1+b1*dr*c1 b1*cr;
    br*c1 Ar];
bg=[b1; 0];
cgphi=[1 0 0];
cgu=-[kP*TD/T1 0 kP*(1-TD/T1)];
dg=[0];
T=[0:0.05:4];
Kreis=ss(Ag, bg, cgphi, dg);
subplot(6,1,[1,2]);
Y=initial(Kreis, [0.1; 0; 0], T);
plot(T, Y, 'b', [0 10], [0 0], 'k:');
hold on;
%
subplot(6,1,[3,4]);
StellgrKreis=ss(Ag, bg, cgu, dg);
U=initial(StellgrKreis, [0.1; 0; 0], T);
plot(T, U, 'b', [0 10], [0 0], 'k:');
hold on;
%
% geregeltes System mit TD1 anstelle von TD
cr=-kP*(1-TD1/T1);
dr=-kP*TD1/T1;
Ag=[A1+b1*dr*c1 b1*cr;
    br*c1 Ar];
bg=[b1; 0];
cgphi=[1 0 0];
cgu=-[kP*TD1/T1 0 kP*(1-TD1/T1)];
dg=[0];
Kreis=ss(Ag, bg, cgphi, dg);
%
subplot(6,1,[1,2]);
Y=initial(Kreis, [0.1; 0; 0], T);
plot(T, Y, 'b--');
latextitle('Verhalten des geregelten Pendels bei $$\phi(0)=0.1$$ rad ');
latexylabel('$$\phi$$ in rad');
xleer;
rotateY;
axis([0 3 -0.08 0.12]);
yticks([0 0.1]);
hold off;
%
subplot(6,1,[3,4]);
StellgrKreis=ss(Ag, bg, cgu, dg);
U=initial(StellgrKreis, [0.1; 0; 0], T);
plot(T, U, 'b--');
latexylabel('$$u$$ in N');
rotateY;
latexxlabel('$$t$$ in s');
axis([0 3 -0.2 0.08]);
xticks([0 1 2 3]);
yticks([-0.2 0]);
hold off;

%%     Bilder auf .EPS- und .PS-Dateien
%
Epsbildc15('wokentw');
