%   Script FKLMOTR.M zur Herstellung von Bildern
%
%   Reglerentwurf im Bodediagramm
%   Entwurf einer Geschwindigkeitsregelung
%   für einen Gleichstrommotor
%
% J. Lunze
% 7.12.1995
% für 2. Auflage: 1.11.1998
% für 3. Auflage: 4.4.2001
% für 4. Auflage: 21.9.2003
% für 5. Auflage: 18.4.2005
% für 6. Auflage: 2.12.2006 (veränderte Motordaten)
% für 7. Auflage: 28.4.2008
% für 8. Auflage: 2.2.2010
% für 9. Auflage: 23.4.2012
% für 10. Auflage: 25.2.2014
% für 11. Auflage: 3.5.2016
% für 12. Auflage: 05.08.2019
%
%   
clear
close all
%
%%    Gleichstrommotor als Regelstrecke
%
figure(1);
Bspmotor;   % liefert das Modell "Motor"
phimin=-300;      % Grenze für Darstellung des Phasenwinkels in allen Diagrammen
%
%  Drehzahlgeber
Td = 0.008;       % s
Drehzahlgeber=tf([1], [Td 1]);
MotorMitDrehzahl=series(Motor, Drehzahlgeber);
O=logspace(-1, 3, 300);
[BetragMotor, PhaseMotor]=bode(MotorMitDrehzahl, O);
%
subplot(2, 1, 1);
semilogx(O, 20*log10(BetragMotor(:)), 'b');
hold on
latextitle('Gleichstrommotor');
xticks([0.1 1 10 100 1000]);
xleer;
latexylabel('$$|G|$$ in dB');
rotateY;
axis([0.1 1000 -100 40]);
yticks([-100 -50 0]);
grid('on');
%
subplot(2, 1, 2);
semilogx(O, PhaseMotor(:),'b');
latexxlabel('$$\omega$$ in rad/s');
grid('on');
latexylabel('$$\phi$$ in Grad');
rotateY;
axis([0.1 1000 phimin 0]);
xticks([0.1 1 10 100 1000]);
yticks([-360 -270 -180 -90 0]);
grid('on');

%%   Verwendung eines PI-Reglers
%    mit kp=4, TN=0.0777
%
figure(2);
%  Reglerentwurf für verzögerte Störung
%
%  Modellerweiterung:
%   [zd, nd] = Übertragungsfunktion Drehmoment -> Drehzahl
%   [zdappr, ndappr] = PT1-Approximation von [zd, nd]
%   
zd=-1/2/pi*[LA RA/kT];
nd=ns;
zdappr=[-0.145];
ndappr=[0.09 1];
MotorStoerUFappr=tf(zdappr, ndappr);
TSigma=0.09;
% Güteforderung: schneller Störabbau mit kleinem em und delta e
%   Summenzeitkonstante der Störung: 0.09
%
% Entwurf für kleinstmögliches delta e, also d=0.7
%  ->  Knickpunktabstand > 2
%
%    Verwendung eines PI-Reglers
%    mit kp=6, TN=0.0777
kP=6;
TI=0.0777;
zr=kP*[TI 1];
nr=[TI 0];
PIRegler=tf(zr, nr);
offeneKette=series(MotorMitDrehzahl, PIRegler);
[BetragGo, PhaseGo]=bode(offeneKette, O);
[BetragG, PhaseG]=bode(MotorMitDrehzahl, O);
subplot(2, 1, 1);
semilogx(O, 20*log10(BetragGo(:)), 'b', O, 20*log10(BetragG(:)), 'b--');
axis([0.1 1000 -100 40]);
yticks([-100 -50 0]);
latextitle('Regelstrecke mit PI-Regler');
latexylabel('$$|G_0|$$ in dB');
rotateY;
xleer;
grid('on');
%
subplot(2, 1, 2);
semilogx(O, PhaseGo(:), 'b', O, PhaseG(:), 'b--');
grid('on');
latexxlabel('$$\omega$$ in rad/s');
latexylabel('$$\phi_0$$ in Grad');
rotateY;
axis([0.1 1000 phimin 0]);
xticks([0.1 1 10 100 1000]);
yticks([-360 -270 -180 -90 0]);

%%                  
%
%   Übergangsfunktion des PI-geregelten Motors
%
%   
figure(3);
subplot(3,1,[1,2]);
Regelkreis=feedback(offeneKette, 1);
[Y, T]=step(Regelkreis);
plot(T, Y, 'b');
latextitle('F\"uhrungs\"ubergangsfunktion des Motors');
latexxlabel('$$t$$ in s');
latexylabel('$$h_{\rm w}$$');
rotateY;
grid('on');
axis([0 1.5 0 1.2]);
xticks([0 0.5 1 1.5]);
yticks([0 0.5 1]);


%%                  
%
%    Schar der Übergangsfunktionen des geregelten Motors
%    bei verändertem Trägheitsmoment  J = 0.1 ... 0.5
%
%   
figure(4);
subplot(2,1,1)
T=0:0.02:6;
for J=0.1:0.1:0.5
  zsneu = [kT/2/pi];
  nsneu = [LA*J (LA*kL + RA*J) (RA*kL + km*kT)];
  Motorneu=tf(zsneu, nsneu);
  Motorneu1=series(Motorneu, Drehzahlgeber);
  offeneKetteneu=series(Motorneu1, PIRegler);
  Regelkreisneu=feedback(offeneKetteneu, 1);
  Y=step(Regelkreisneu, T);
  plot(T, Y, 'b');
  hold on;
end
latextitle('PI-geregelter Motor f\"ur $$J=0,1 ... 0,5$$');
axis([0 3.5 0 1.3]);
xticks([0 1 2 3 4 5]);
yticks([0 0.5 1]);
latexxlabel('$$t$$ in s');
latexylabel('$$h_{\rm w}$$');
rotateY;
grid('on');

%%                  
%
%  Ortskurve der offenen Kette
%  bei unterschiedlichem Trägheitsmoment
%
%   
figure(5);
for J=0.1:0.05:0.5
  zsneu = [kT/2/pi];
  nsneu = [LA*J (LA*kL + RA*J) (RA*kL + km*kT)];
  Motorneu=tf(zsneu, nsneu);
  Motorneu1=series(Motorneu, Drehzahlgeber);
  offeneKetteneu=series(Motorneu1, PIRegler);
  [Real, Imag]=nyquist(offeneKetteneu, O);
  plot(Real(:), Imag(:), 'b');
  hold on;
end
plot([-3 0.5],[0 0], 'k:');
plot([-3 0.5],[0 0],'k:',[0 0],[-3 0.5],'k:');
plot([-1 -1], [-0.1 0.1], 'k');
%grid('on');
latextitle('Offene Kette bei $$J=0,1 ... 0,5$$');
latexylabel('Im');
rotateY;
latexxlabel('Re');
axis([-3 0.5 -3 0.5]);
axis('square');
xticks([-3 -2 -1 0]);
yticks([-3 -2 -1 0]);
hold off;

%%  Robustheitsabschätzung mit oberer Fehlerschranke
%
figure(6); 
% Berechnung der oberen Fehlerschranke
Wr=logspace(-1, 3, 40);
[GReal, GImag]=nyquist(MotorMitDrehzahl, Wr);
k=0;
for J=0.12:0.03:0.5
  k=k+1;
  zsneu = [kT/2/pi];
  nsneu = [LA*J (LA*kL + RA*J) (RA*kL + km*kT)];
  Motorneu=tf(zsneu, nsneu);
  dG=bode(Motorneu*Drehzahlgeber-MotorMitDrehzahl, Wr);
  BetragFM(:,k)=(abs(dG(:)));
end
GFM=max(BetragFM')';
%
[Real, Imag]=nyquist(offeneKette);
plot(Real(:), Imag(:), 'b--');
hold on
latextitle('Toleranzband f\"ur die Ortskurve der offenen Kette');
latexxlabel('Re');
latexylabel('Im');
rotateY;
hold on;
plot([-3 3],[0 0],'k:', [0 0],[-3 0.5], 'k:');
%                      Robustheitsabschätzung
kW = size(Wr,2);
[BetragR, PhaseR]=bode(PIRegler, Wr);
Betrag=GFM(:).*BetragR(:);
[BetragG, PhaseG, O] = bode(offeneKette, Wr);
for k=1:kW
    k1 = 0;
    RealG = BetragG(k)*cos(PhaseG(k)*pi/180);
    ImagG = BetragG(k)*sin(PhaseG(k)*pi/180);
    for x1=-Betrag(k):Betrag(k)/100:Betrag(k)
        k1=k1+1;
        K1(k1, 1)= x1+RealG;
        K1(k1, 2)= ImagG + sqrt(Betrag(k)^2 - x1^2);
        K1(402-k1, 1) = x1+RealG;
        K1(402-k1, 2) = ImagG - sqrt(Betrag(k)^2 - x1^2);
    end
    plot(K1(:, 1), K1(:, 2), 'b');
end
plot([-1 -1],[-0.05 0.05],'k', [-1.05 -0.95],[0 0],'k');
axis([-1.2 0.7 -1.6 0.3]);
xticks([-1 -0.5 0 0.5 1]);
yticks([-1.5 -1 -0.5 0]);
grid('on');
axis('square');
hold off;

%%  Darstellung des relativen Fehlers des Regelstreckenmodells
%
%   
figure(7);
subplot(2,1,1)
[Betrag, Phase]=bode(MotorMitDrehzahl, O);
relF=GFM(:)./Betrag(:);
semilogx(O, relF(:), 'b');
latextitle('Relativer Fehler des N\"aherungsmodells');
latexxlabel('$$\omega$$ in rad/s');
axis([0.1 1000 0 1]);
xticks([0.1 1 10 100 1000]);
yticks([0 0.5 1]);
latexylabel('$$\bar{G}_{\rm M}$$');
rotateY;
grid('on');

%%    Wiederholung des Bodediagramms der offenen Kette
%
figure(8);
[GoBetrag, GoPhase]=bode(offeneKette, O);
subplot(2, 1, 1);
semilogx(O, 20*log10(GoBetrag(:)), 'b');
latextitle('Bodediagramm der offenen Kette mit PI-Regler $$(k_P=4, T_I=0.0777)$$');
latexylabel('$$|G_0|$$ in dB');
rotateY;
xleer;
grid('on');
axis([0.1 1000 -100 50]);
yticks([-100 -50 0 50]);
%
subplot(2, 1, 2);
semilogx(O, GoPhase(:), 'b');
axis([0.1 1000 phimin 0]);
grid('on');
latexxlabel('$$\omega$$ in rad/s');
xticks([0.1 1 10 100 1000]);
latexylabel('$${\rm arg}(G_0)$$ in Grad');
rotateY;
xticks([0.1 1 10 100 1000]);
yticks([-360 -270 -180 -90 0]);

%%   Übergangsfunktionen des PI-geregelten Motors
%
figure(9);
subplot(2,1,1)
[Y, T]=step(Regelkreis, 1.5);
plot(T, Y, 'b');
hold on;
Stoerverhalten=feedback(1, offeneKette);
Stoerverhalten2=series(Stoerverhalten, MotorStoerUF);
[Y, T]=step(Stoerverhalten2);
plot(T, Y, 'b');
latextitle('\"Ubergangsfunktionen des PI-geregelten Motors');
latexxlabel('$$t$$ in s');
latexylabel('$$h_{\rm w}, h_{\rm d}$$');
rotateY;
grid('on');
axis([0 1.5 -0.2 1.2]);
xticks([0 0.5 1 1.5 2]);
yticks([0 0.5 1]);

%%   Schar der Übertragungsfunktionen des Motors
%    bei verändertem Trägheitsmoment  J = 0.1 ... 0.5
%
%   
figure(10);
% Näherungsmodell
[BetragNM, PhaseNM]=bode(MotorMitDrehzahl, Wr);
subplot(2,1,1)
semilogx(Wr, 20*log10(BetragNM(:)), 'b--');
hold on
%
figure(11)
[RealNM, ImagNM]=nyquist(MotorMitDrehzahl, Wr);
plot(RealNM(:), ImagNM(:), 'b');
hold on
k=0;                    % Zähler für das berechnete Modell
for J=0.15:0.05:0.5
  k=k+1;
  zsneu = [kT/2/pi];
  nsneu = [LA*J (LA*kL + RA*J) (RA*kL + km*kT)];
  Motorneu=tf(zsneu, nsneu);
  figure(10)
  [Betrag, Phase]=bode(Motorneu*Drehzahlgeber, Wr);
  subplot(2,1,1)
  semilogx(Wr, 20*log10(Betrag(:)), 'b');
  subplot(2,1,2);
  [BetragFM, PhaseFM]=bode(MotorMitDrehzahl-Motorneu*Drehzahlgeber, Wr);
  BetragFehler(k,1,:)=BetragFM(1,1,:);       %  Sammeln aller Fehlermodelle
  semilogx(Wr, 20*log10(BetragFM(:)), 'b');
  hold on;
  figure(11)
  [Real, Imag]=nyquist(Motorneu*Drehzahlgeber);
  plot(Real(:), Imag(:), 'b');
end
figure(10)
subplot(2,1,1)
latextitle('Motor f\"ur $$J=0.1 ... 0.5$$');
axis([0.1 1000 -120 -20]);
yticks([-100 -60 -20]);
xleer;
latexylabel('$$|G|$$ in dB');
rotateY;
grid('on');
hold off
%
subplot(2,1,2)
semilogx(Wr, 20*log10(GFM(:)), 'b--');
latexxlabel('$$\omega$$ in rad/s');
xticks([0.1 1 10 100 1000]);
latexylabel('$$|\delta G|$$ in dB');
rotateY;
grid('on');
axis([0.1 1000 -120 -20]);
xticks([0.1 1 10 100 1000]);
yticks([-100 -60 -20]);
hold off
%  obere Fehlerschranke G_A
GA(1,1,:)=max(BetragFehler);

figure(11)
latextitle('Motor f\"ur $$J=0.1 ... 0.5$$');
axis([-0.005 0.03 -0.03 0.005]);
axis('square');
xticks([-0.005 0 0.005 0.01 0.015 0.02 0.025]);
yticks([-0.03 -0.02 -0.01 0 0.01]);
latexxlabel('Re');
latexylabel('Im');
rotateY;
grid('on');
hold off

%%  Schar der Ortskurven des Motors und Toleranzband
%
%   
figure(12);
% Näherungsmodell
plot(RealNM(:), ImagNM(:),'b--');
hold on
plot([-3 0.5],[0 0], 'k:', [0 0],[-3 0.5],'k:');
%                      Robustheitsabschätzung
kW = size(Wr,2);
for k=1:kW
    k1 = 0;
    for x1=-GFM(k):GFM(k)/100:GFM(k)
        k1=k1+1;
        K1(k1, 1)= x1+RealNM(k);
        K1(k1, 2)= ImagNM(k) + sqrt(GFM(k)^2 - x1^2);
        K1(402-k1, 1) = x1+RealNM(k);
        K1(402-k1, 2) = ImagNM(k) - sqrt(GFM(k)^2 - x1^2);
    end
    plot(K1(:, 1), K1(:, 2), 'b');
end
plot([-1 -1], [-5 5], 'b');
plot([-1 -1], [-5 5], 'k:');
latextext(-1, 0, '|');
axis([-0.02 0.08 -0.07 0.03]);
xticks([-0.02 0 0.02 0.04 0.06 0.08]);
yticks([-0.06 -0.04 -0.02 0 0.02 0.04]);
axis('square');
grid('on');
latextitle('Toleranzband f\"ur die Ortskurve des Motors');
latexxlabel('Re');
latexylabel('Im');
rotateY;
hold off;

%%   Darstellung der Bedingung für robuste Stabilität
%
%   
figure(13)
%   Prüfung der Robustheitsbedingung  GA < |\frac{1+Ghat*K}{K}|
subplot(2,1,1)
semilogx(Wr, 20*log10(GA(:)),'b');
hold on
Stellverhalten=feedback(PIRegler,MotorMitDrehzahl);
GKGBetrag=bode(Stellverhalten, Wr);
k=0;
for l=1:length(GKGBetrag)
    k=k+1;
    GKGBetraginv(k)=1/GKGBetrag(k);
end
GKGmax=20*log10(GKGBetraginv(end));
semilogx([Wr 10000], [20*log10(GKGBetraginv(:))' GKGmax],'b--'); %Verlängerung der Kurve
grid('on');
latexxlabel('$$\omega$$ in rad/s');
latextext(1, -50, '$$\left|\bar{G}_{\rm A}({\rm j} \omega)\right|$$')
latextext(150, -40, '$$\left|\frac{1+\hat{G}({\rm j}\omega)K({\rm j}\omega)}{K({\rm j}\omega)}\right|$$')
rotateY;
latextitle('Robustheitsanalyse des Drehzahlregelkreises');
axis([0.1 10000 -70 0]);
xticks([0.1 1 10 100 1000]);
yticks([-50 -20 0]);

%%   Wurzelortskurve des P-geregelten Motors
%
%   
figure(14)
subplot(3,3,[1, 2, 4, 5]);
rlocprintblue(MotorUF);
latextitle('P-geregelter Gleichstrommotor');
axis([-80 10 -45 45]);
xticks([-80 -40 0]);
yticks([-40 0 40]);
rotateY;


%%    Bodediagramm des Motors ohne Drehzahlgeber
%
%   
figure(15);
O=logspace(-1, 3, 300);
[BetragMotor, PhaseMotor]=bode(Motor, O);
subplot(2, 1, 1);
semilogx(O, 20*log10(BetragMotor(:)), 'b');
hold on
latextitle('Gleichstrommotor');
xticks([0.1 1 10 100 1000]);
xleer;
latexylabel('$$|G|$$ in dB');
rotateY;
yticks([-80 -40 0]);
axis([0.1 1000 -100 0]);
grid('on');
%
subplot(2, 1, 2);
semilogx(O, PhaseMotor(:), 'b');
latexxlabel('$$\omega$$ in rad/s');
latexylabel('$$\phi$$ in Grad');
rotateY;
axis([0.1 1000 phimin 0]);
xticks([0.1 1 10 100 1000]);
yticks([-360 -270 -180 -90 0]);
grid('on');


%%  Wiederholung der Robustheitsabschätzung aus Figure(6)
%   mit größerer oberer Fehlerschranke
%
figure(16)
% Berechnung der oberen Fehlerschranke
Wr=logspace(-1, 3, 40);
[GReal, GImag]=nyquist(MotorMitDrehzahl, Wr);
k1=0;
BetragFM=[];
for J=0.12:0.03:1   % statt 0.5
  k1=k1+1;
  zsneu = [kT/2/pi];
  nsneu = [LA*J (LA*kL + RA*J) (RA*kL + km*kT)];
  Motorneu=tf(zsneu, nsneu);
  dG=bode(Motorneu*Drehzahlgeber-MotorMitDrehzahl, Wr);
  BetragFM(:,k1)=(abs(dG(:)));
end
GFM=max(BetragFM')';
%
[Real, Imag]=nyquist(offeneKette);
plot(Real(:), Imag(:), 'b--');
hold on
axis([-1.2 0.7 -1.6 0.3]);
grid('on');
latextitle('Toleranzband f\"ur die Ortskurve der offenen Kette');
latexxlabel('Re');
latexylabel('Im');
rotateY;
plot([-3 3],[0 0], 'k:', [0 0],[-3 0.5], 'k:');
%                      Robustheitsabschätzung
kW = size(Wr,2);
[BetragR, PhaseR]=bode(PIRegler, Wr);
Betrag=GFM(:).*BetragR(:);
[BetragG, PhaseG, O] = bode(offeneKette, Wr);
for k=1:kW
    k1 = 0;
    RealG = BetragG(k)*cos(PhaseG(k)*pi/180);
    ImagG = BetragG(k)*sin(PhaseG(k)*pi/180);
    for x1=-Betrag(k):Betrag(k)/100:Betrag(k)
        k1=k1+1;
        K1(k1, 1)= x1+RealG;
        K1(k1, 2)= ImagG + sqrt(Betrag(k)^2 - x1^2);
        K1(402-k1, 1) = x1+RealG;
        K1(402-k1, 2) = ImagG - sqrt(Betrag(k)^2 - x1^2);
    end
    plot(K1(:, 1), K1(:, 2), 'b');
end
plot([-1 -1],[-0.05 0.05],'k', [-1.05 -0.95],[0 0],'k');
axis([-1.2 0.7 -1.6 0.3]);
xticks([-1 -0.5 0 0.5 1]);
yticks([-1.5 -1 -0.5 0]);
grid('on');
axis('square');
hold off;

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

