%   Script OPTRDE.M zur Herstellung von Bildern
%
%   Optimalreglerentwurf für den Dampferzeuger
%
% verwendet BSPDE.M, EPSBILD.M
%
% J. Lunze
% 29.9.1996            Daten wie 16.3.1996
% für 2. Auflage: 16.12.2001
% für 3. Auflage: 19.4.2004
% für 5. Auflage: 18.10.2007
% für 6. Auflage: 5.1.2010
% für 7. Auflage 25.5.2012
% für 10. Auflage: 9.9.2019
%
%
clear
close all
Bspde;
figure(1);
Q=C'*C;
R=eye(m,m);
x0=ones(n, 1);
T=[0:1:150];
for q=[0.1, 1, 10, 30]
  K=lqr(A, B, q*Q, R);
  subplot(2, 2, 1);
  DampferzeugerAusgang=ss(A-B*K, B, C, D);
  [Y, T]=initial(DampferzeugerAusgang, x0, T);
  plot(T, Y(:,1,1), 'b', [0 150], [0 0], 'k:', [-10, 0], [0, 0], 'b', [0 0],[-0.5 2.5],'k:');
  hold on
  subplot(2, 2, 3);
  plot(T, Y(:,2,1), 'b', [0 150], [0 0], 'k:', [-10, 0], [0, 0], 'b', [0 0],[-0.5 0.5],'k:');
  hold on
  DampferzeugerStellgr=ss(A-B*K, B, K, D);
  [U, T]=initial(DampferzeugerStellgr, x0, T);
  subplot(2, 2, 2);
  plot(T, U(:,1,1), 'b', [0 150], [0 0], 'k:', [-10, 0], [0, 0], 'b', [0 0],[-0.5 2.5],'k:');
  hold on
  subplot(2, 2, 4);
  plot(T, U(:,2,1), 'b', [0 150], [0 0], 'k:', [-10, 0], [0, 0], 'b', [0 0],[-0.5 0.5],'k:');
  hold on
end
subplot(2, 2, 1);
axis([0 60 -0.005 0.05]);
xleer;
latexylabel('$$y_1$$');
rotateY;
yticks([0 0.02 0.04]);
latextitle('Dampferzeugers mit Optimalregler');
hold off;
subplot(2, 2, 3);
axis([0 60 -0.005 0.002]);
latexxlabel('$$t$$ in min');
latexylabel('$$y_2$$');
rotateY;
yticks([-0.004 -0.002 0 0.002]);
xticks([0 20 40 60]);
hold off;
xticks([0 50]);
%
subplot(2, 2, 2);
axis([0 60 -0.03 0.2]);
xleer;
yticks([0 0.1 0.2]);
latexylabel('$$u_1$$');
rotateY;
latextitle('$$Q=qI, q=0.1, 1, 10, 30$$');
hold off;
%
subplot(2, 2, 4);
axis([0 60 -0.0002 0.004]);
latexxlabel('$$t$$ in min');
latexylabel('$$u_2$$');
rotateY;
yticks([0  0.004]);
xticks([0 20 40 60]);
hold off;
xticks([0 50]);


%%   Führungsverhalten mit Vorfilter
%
%
figure(2)
q=10;
K=lqr(A, B, q*Q, R);
Ag=A-B*K;
Vorfilter=ss(Ag, B, C, D);
Kyw=dcgain(Vorfilter);
V=inv(Kyw);
DampferzeugerFuehrg=ss(Ag, B*V, C, D);
Y=step(DampferzeugerFuehrg, T);
subplot(2, 2, 1);
plot(T, Y(:,1,1), 'b', [0 150],[0 0],'k:');
latextitle('$$Q=10 I$$');
latexylabel('$$h_{\rm w11}$$');
rotateY;
yticks([0 1]);
xleer;
axis([0 150 -0.2 1.2]);
subplot(2, 2, 2);
plot(T, Y(:,1,2), 'b', [0 150],[0 0],'k:');
latextitle('');
latexylabel('$$h_{\rm w12}$$');
rotateY;
yticks([0 1]);
xleer;
axis([0 150 -0.2 1.2]);
subplot(2, 2, 3);
plot(T, Y(:,2,1), 'b', [0 150],[0 0],'k:');
latextitle('');
latexylabel('$$h_{\rm w21}$$');
rotateY;
yticks([0 1]);
latexxlabel('$$t$$ in min');
axis([0 150 -0.2 1.2]);
xticks([0 100]);
subplot(2, 2, 4);
plot(T, Y(:,2,2), 'b', [0 150],[0 0],'k:');
latextitle('');
latexylabel('$$h_{\rm w22}$$');
rotateY;
yticks([0 1]);
latexxlabel('$$t$$ in min');
axis([0 150 -0.2 1.2]);
xticks([0 100]);


%%
%
%   Verbesserte Wichtung der Eingänge
%   Q=C'*diag(1, 20)*C
%
figure(3);
Q=C'*diag([1, 30])*C;
K=lqr(A, B, Q, R);
Ag=A-B*K;
Vorfilter=ss(Ag, B, C, D);
Kyw=dcgain(Vorfilter);
V=inv(Kyw);
DampferzeugerFuehrg2=ss(Ag, B*V, C, D);
Y=step(DampferzeugerFuehrg2, T);
subplot(2, 2, 1);
plot(T, Y(:,1,1), 'b', [0 150],[0 0],'k:');
latextitle('$$Q=C*diag(1, 30)*C$$');
latexylabel('$$h_{\rm w11}$$');
rotateY;
yticks([0 1]);
xleer;
axis([0 150 -0.2 1.2]);
subplot(2, 2, 2);
plot(T, Y(:,1,2), 'b', [0 150],[0 0],'k:');
latextitle('');
latexylabel('$$h_{\rm w12}$$');
rotateY;
yticks([-1 0]);
xleer;
axis([0 150 -1.3 0.1]);
subplot(2, 2, 3);
plot(T, Y(:,2,1), 'b', [0 150],[0 0],'k:');
latextitle('');
latexylabel('$$h_{\rm w21}$$');
rotateY;
yticks([0 1]);
latexxlabel('$$t$$ in min');
axis([0 150 -0.2 1.2]);
xticks([0 100]);
subplot(2, 2, 4);
plot(T, Y(:,2,2), 'b', [0 150],[0 0],'k:');
latextitle('');
latexylabel('$$h_{\rm w22}$$');
rotateY;
yticks([0 1]);
latexxlabel('$$t$$ in min');
axis([0 150 -0.2 1.2]);
xticks([0 100]);


%                         Stellgrößen
figure(4);
DampferzeugerStellgr2=ss(Ag, B*V, K, zeros(m,m));
U=step(DampferzeugerStellgr2,T);
latextitle('Stellgr\"o\ssen');
subplot(2, 2, 1);
plot(T, U(:,1,1), 'b', [0 150],[0 0],'k:');
latextitle('');
latexylabel('$$u_1$$');
rotateY;
xleer;
axis([0 150 -3 1.7]);
yticks([-2 0 2]);
%
subplot(2, 2, 2);
plot(T, U(:,1,2), 'b', [0 150],[0 0],'k:');
latextitle('');
latexylabel('$$u_1$$');
rotateY;
xleer;
axis([0 150 -3 1.7]);
yticks([-2 0 2]);
%
subplot(2, 2, 3);
plot(T, U(:,2,1), 'b', [0 150],[0 0],'k:');
latextitle('');
latexylabel('$$u_2$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([0 150 -0.1 1.1]);
xticks([0 100]);
yticks([0 1]);
%
subplot(2, 2, 4);
plot(T, U(:,2,2), 'b', [0 150],[0 0],'k:');
latextitle('');
latexylabel('$$u_2$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([0 150 -0.1 1.1]);
xticks([0 100]);
yticks([0 1]);
 


%%   Ortskurve der Determinante der R\"uckf\"uhrdifferenzmatrix
%   für Q=q*C'*C mit q=1,  10, 30
%
figure(5);
W=logspace(-4, 0, 200);
for q=[1 10 30]
    K=lqr(A, B, q*C'*C, R);
    Rueckfuehrdifferenzmatrix=ss(A, B, K, eye(m,m));
    [Real, Imag]=nyquist(Rueckfuehrdifferenzmatrix, W);
    for k=1:length(Real)
       detF(k)=det(Real(:,:,k)+i*Imag(:,:,k));
    end
    plot(real(detF), imag(detF), 'b', real(detF), -imag(detF),'b--');
    hold on;
    plot(detF(1:1), 0, 'bo');
end
k1=0;
for x1=-1:1/80:1
     k1=k1+2;
     K1(k1, 1)= x1;
     K1(k1, 2)= sqrt(1 - x1^2);
     K1(k1+1, 1) = x1;
     K1(k1+1, 2) = -sqrt(1 - x1^2);
end;
plot(K1(:,1), K1(:,2), 'b', K1(:,2), K1(:,1), 'b', [-2 14], [0 0], 'k:', [0 0], [8 -8], 'k:');
axis([-2 14 -8 8]);
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Ortskurve von $$\det F$$ der Dampferzeugerregelung');
text(4, -4, '\leftarrow');
text(4,  4, '\rightarrow');
latextext(1.5, -1.5, '$$q=1$$');
latextext(4.5, -1.5, '$$q=10$$');
latextext(11, -1.5, '$$q=30$$');
hold off;
xticks([0 5 10 15]);
yticks([-5 0 5]);



%%   Ortskurve der Determinante der R\"uckf\"uhrdifferenzmatrix
%   für Q=10*C'*C
%
figure(6);
W=logspace(-5, 0, 400);
K=lqr(A, B, 200*C'*C, R);
Rueckfuehrdifferenzmatrix=ss(A, B, K, eye(m,m));
[Real, Imag]=nyquist(Rueckfuehrdifferenzmatrix, W);
for k=1:length(Real)
    detF(k)=det(Real(:,:,k)+i*Imag(:,:,k));
end
plot(real(detF), imag(detF), 'b', real(detF), -imag(detF),'b--');
hold on;
plot(detF(1,1), 0, 'bo');
plot([0 0],[-30 60],'k:', [-30 60],[0 0],'k:');
axis([-5 55 -30 30]);
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
latextitle('Ortskurve von det F der Dampferzeugerregelung');
text(25, -25, '\leftarrow');
text(25,  25, '\rightarrow');
plot(1, 0, 'b+');
latextext(2, 1, '1');
hold off;
xticks([0 20 40 60]);
yticks([-20 0 20]);


%%     Ausgabe der Bilder als .eps-Dateien
%
Epsbildc15('optrde');
