%  Example 10.5.  Mu-synthesis controller design.
%  Note the file dk_defin.m must be inserted into the subs folder in the 
%  mutools toolbox.  This file provides data for the DK iteration algorithm.
clear
%  Define the plant model.
A=[0  1
   0.5 0];
Bu=[0 0.1]'; Bw=[0 0.1]';
Cm=[1 0]; Cy=[1 0];
%  Define the bounds.
Ww=0.5;
Wv=0.05;
%  Generate the normalized plant model in standard form.
B1=[Ww*sqrt(2)*Bw zeros(2,1) Bu]; B2=Bu;
C1=[5*Cy; zeros(1,2); zeros(1,2)]; C2=Cm;
D11=zeros(3);
D12=[0;0.01;0.5];D21=[0 Wv*sqrt(2) 0];D22=0;
P=pck(A,[B1 B2],[C1;C2],[D11 D12;D21 D22]);
%  Set variables for DK-iteration algorithm.
  %  Number of measurements.
nmeas=1;
  %  Number of controls.
ncont=1;
  %  Block structure for the perturbation.
blk=[2 2;1 1];
  %  Frequency vector for evaluating the SSV.
omega=logspace(-1,3,30);
%  Perform DK-iteration.
DK_DEF_NAME='my_defin';
dkit
%  The DK-iteration algorithm yields a controllers for each iteration.
%    Request the desired controller.
iter=input('Enter the iteration number of the desired controller ', 's');
K=eval(['k_dk',iter,'suffix']);
%  Generate the state model of the mu-synthesis controller.
[Ac,Bc,Cc,Dc]=unpck(K);
%  Evaluate the robustness of the mu-synthesis controller.
%  Generate the state model of the loop transfer function.
[n,m]=size(Bc);
Al=[  A   zeros(2,n)
    Bc*Cm   Ac    ];
Bl=[    Bu
    zeros(n,1)];
Cl=[zeros(1,2) -Cc];
%  Generate Nyquist plot.
set(0,'DefaultAxesFontName','times')
set(0,'DefaultAxesFontSize',16)
set(0,'DefaultTextFontName','times')
figure(1)
clf
axes('Position',[.15 .2 .75 .75]);
w=logspace(-2,3,100);
[re,im]=nyquist(Al,Bl,Cl,0,1,w);
plot(re,im,'r-',re,-im,'r-',-1,0,'r+',re(50),im(50),...
   'r<',re(50),-im(50),'r>')
xlabel('Real axis')
ylabel('Imaginary axis')
axis([-2 2 -2 2])
axis('square')
grid
%  Generate the closed loop model for simulation.
Acl=[A B2*Cc;Bc*C2 Ac];
Bcl=[B1;Bc*D21];
Ccl=[C1 D12*Cc];
Dcl=zeros(3);
%  Define the time vector for simulation.
t=0:.02:5;
%  Generate the step responses for the closed loop system.
[y1,x]=step(Acl,Bcl,Ccl,Dcl,1,t);
[y2,x]=step(Acl,Bcl,Ccl,Dcl,2,t);
[y2,x]=step(Acl,Bcl,Ccl,Dcl,2,t);
%  Plot the step responses.
figure(2)
clf
subplot(211)
plot(t,y1(:,1)/5,'r-',t,y2(:,1)/5,'r--')
xlabel('Time (sec)')
ylabel('Angle (deg)')
grid
legend('{\ity}({\itt}), {\itw}({\itt}) = 1({\itt})','{\ity}({\itt}), {\itv}({\itt}) = 1({\itt})')
subplot(212),
plot(t,y1(:,2)*100,'r-',t,y2(:,2)*100,'r--')
xlabel('Time (sec)')
ylabel('Control torque (N-m)')
grid
legend('{\itu}({\itt}), {\itw}({\itt}) = 1({\itt})','{\itu}({\itt}), {\itv}({\itt}) = 1({\itt})')
