%
% ACTDEMO Digital Hydraulic Actuator Design Demonstration
%

% R. Y. Chiang & M. G. Safonov 3/88
% Copyright (c) 1988 by the MathWorks, Inc.
% All Rights Reserved.
% ---------------------------------------------------------------
clear
clc
disp('    ')
disp(' ')
disp('          << Demo #1: SISO Digital Hydraulic Actuator Design Example >>')
disp('  ')
disp(' ')
disp('  ')
disp('                         -------             -----        -----')
disp('        +        /      | H-Inf |     /     |     |      | Hyd.|') 
disp('      --->(X)----  -----| Cont- | ----  --->| ZOH |----->|     |--------->')
disp('           ^ - Ts:0.01  | roller|   Ts:0.01 |     |      | ACT.|    |')
disp('           |      sec    -------       sec   -----        -----     |')
disp('           |                                                        |')
disp('           |________________________________________________________|')
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
format short e
num = 9000;
den = [1 30 700 1000];
[a,b,c,d] = tf2ss(num,den);
clc
disp('  ')
disp('  ')
disp('   Hydraulic Actuator Open Loop:')
disp(' ')
disp('   ')
disp('                                 9000')
disp('             G(s) = -----------------------------')
disp('                     s^3 + 30 s^2 + 700 s + 1000')
disp('  ')
disp(' ')
disp('      Poles of the open loop plant:')
poleg = roots(den)
disp('  ')
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
clc
disp(' ')
disp(' ')
disp('                              << W-Plane Design >>')
disp('  ')
disp('  ')
disp('           Let us attach the plant with a Z.O.H and convert it into ')
disp(' ')
disp('           W-plane (sampling period: 0.01 sec) ....')
disp(' ')
disp('  ')
disp('         ------------------------------------------------------------')
disp('           [az,bz] = c2d(a,b,0.01);   % Attach the plant with Z.O.H')
disp('         ------------------------------------------------------------')
[az,bz] = c2d(a,b,0.01);
disp('           [ag,bg,cg,dg] = bilin(az,bz,c,d,-1,`Tustin`,0.01);')
disp('                           % Convert the plant with Z.O.H to W-plane')
disp('         ------------------------------------------------------------')
[ag,bg,cg,dg] = bilin(az,bz,c,d,-1,'Tustin',0.01);
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
clc
disp(' ')
disp(' ')
disp('     - - - Computing Bode plot of the open loop plant (in s & w) - - -')
w = logspace(-3,5,50);
svg = bode(a,b,c,d,1,w); svg = 20*log10(svg);
svw = bode(ag,bg,cg,dg,1,w); svw = 20*log10(svw);
disp(' ')
disp(' ')
disp(' ')
disp(' ')
disp('                             (strike a key to see the plot ...)')
pause
clg
plot(111)
semilogx(w,svg,w,svw)
title('SISO Hydraulic Actuator Open Loop (in s & w domain)')
xlabel('Frequency - Rad/Sec')
ylabel('SV - db')
text(150,-30,'Nyquist Freq.: 100pi')
grid
pause
clc
disp('  ')
disp('                       << Design Specifications >>  ')
disp(' ')
disp('      1). Robustness Spec. : closed loop bandwidth -- 30 r/s')
disp('          Associated Weighting:')
disp('     ')
disp('                    -1     3.16(1 + s/300)')
disp('                  W3(s) = -----------------')
disp('                             (1 + s/10)')
disp('   ')
disp(' ')
disp('      2). Performance Spec.: sensitivity reduction of at least 100:1')
disp('                             up to approx. 1 r/s')
disp('          Associated Weighting:')
disp(' ')          
disp('                   -1         -1   0.01(1 + s)^2')
disp('                  W1(s) = Gam  * ---------------')
disp('                                   (1 + s/30)^2')
disp('   ')
disp('          where "Gam" in our design goes from 1 --> 1.5')
nuw3i = 3.16*[1/300 1]; dnw3i = [1/10 1];
svw3i = bode(nuw3i,dnw3i,w); svw3i = 20*log10(svw3i);
nuw1i = 0.01*[1 2 1]; dnw1i = conv([1/30 1],[1/30 1]);
svw1i = bode(nuw1i,dnw1i,w); svw1i = 20*log10(svw1i);
disp('  ')
disp('  ')
disp('                   (strike a key to see the plot of the weightings ...)')
pause
clg
plot(111)
semilogx(w,svw1i,w,svw3i)
grid
title('MIMO LSS Design Example -- Design Specifications')
xlabel('Frequency - Rad/Sec')
ylabel('1/W1 & 1/W3 - db')
text(0.005,-20,'Sensitivity Spec.-- 1/W1(s)')
text(100,0,'Robustness Spec.')
text(1000,-10,'1/W3(s)')
pause
clc
disp('                      << Problem Formulation >>')
disp(' ')
disp('      Form an augmented plant P(s) with these two weighting functions:')
disp(' ')
disp('                 1). W1 penalizing error signal "e"')
disp('  ')
disp('                 2). W3 penalizing plant output "y"')
disp(' ')
disp('      and find a stabilizing controller F(s) such that the Hinf-norm')
disp('      of TF Ty1u1 is minimized and less than one, i.e.')
disp('  ')
disp('                         min |Ty1u1|   < 1,')
disp('                         F(s)       inf')
disp(' ')
disp('      where ')
disp('                       |               -1|')                         
disp('               Ty1u1 = | Gam*W1*(I + GF) |  = | Gam * W1 * S  |')
disp('                       |               -1|    |  W3 * (I - S) |')
disp('                       |  W3*GF*(I + GF) |')
disp('  ')
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
clc
disp('  ')
disp(' ')
disp('                              << DESIGN PROCEDURE >>')
disp('  ')
disp('            * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')
disp('            *    [Step 1]. Do plant augmentation (run AUGTF.M or      *')
disp('            *              AUGSS.M)                                   *')
disp('            *                                                         *')
disp('            *    [Step 2]. Do H-Inf synthesis (run HINF.M)            *')
disp('            *                                                         *')
disp('            *    [Step 3]. Redo the plant augmentation by setting     *')
disp('            *              new "Gam" --> 1.5 and rerun HINF.M         *')
disp('            * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
clc
disp('  ')
disp('           Assign the cost coefficient "Gam" --> 1 ')
disp('       ')
disp('           this will serve as the baseline design ....')
disp('  ')
disp('    ---------------------------------------------------------------')
disp('     % Plant augmentation of the actuator:')
disp('       w1 = [Gam*dnw1i;nuw1i];')
disp('       w3 = [dnw3i;nuw3i];')
disp('       w2 = [];')
disp('       ss_g  = mksys(ag,bg,cg,dg);')
disp('       [TSS_]=augtf(ss_g,w1,w2,w3);')
disp('    ---------------------------------------------------------------')
Gam = input('           Input cost coefficient "Gam" = ');
w1 = [Gam*dnw1i;nuw1i];
w3 = [dnw3i;nuw3i];
w2 = [];
ss_g = mksys(ag,bg,cg,dg);
[TSS_] = augtf(ss_g,w1,w2,w3);
disp('  ')
disp('     - - - State-Space TSS_:=(A,B1,..,D21,D22) is ready for')
disp('           the Small-Gain problem - - -')
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
clc
disp(' ')
disp(' ')
disp('   --------------------------------------------------------------- ')
disp('    [ss_cp,ss_cl,hinfo] = hinf(TSS_);')
disp('    [acp,bcp,ccp,dcp] = branch(ss_cp);')
disp('    [acl,bcl,ccl,dcl] = branch(ss_cl);')
disp('                    % Running file HINF.M for H-Inf optimization')
disp('   ---------------------------------------------------------------')
[ss_cp,ss_cl,hinfo] = hinf(TSS_);
[acp,bcp,ccp,dcp] = branch(ss_cp);
[acl,bcl,ccl,dcl] = branch(ss_cl);
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
pltopt           % Preparing singular values for plotting
svw1i1 = svw1i; hsvs1 = svs; hsvt1 = svt; hsvtt1 = svtt;
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
clc      
disp('  ')
disp('  ')
disp('  ')
disp('     After a few iterations, we found a new Gam of 1.5 can push the')
disp('  ')
disp('     H-Inf cost function close to its limit. ')
disp('  ')
disp('  ')
disp('            Input "Gam" --> 1.5, and try HINF again .....')
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
Gam = input('           Input cost coefficient "Gam" = ');
w1 = [Gam*dnw1i;nuw1i];
[TSS_] = augtf(ss_g,w1,w2,w3);
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
[ss_cp,ss_cl,hinfo] = hinf(TSS_);
[acp,bcp,ccp,dcp] = branch(ss_cp);
[acl,bcl,ccl,dcl] = branch(ss_cl);
disp('  ')
disp('  ')
disp('                                 (strike a key to continue ...)')
pause
pltopt
svw1i2 = svw1i; hsvs2 = svs; hsvt2 = svt; hsvtt2 = svtt;
disp('  ')
disp('  ')
disp('                  (strike a key to see the plots of the comparison ...)')
pause
clg
plot(111)
semilogx(w,svw1i1,w,hsvs1,w,svw1i2,w,hsvs2)
title('H-Inf W-Plane Actuator Design -- 1/W1 & Sensitivity Func.')
xlabel('Frequency - Rad/Sec')
ylabel('SV - db')
grid
text(0.002,10,'H-Inf (Gam = 1) ---> H-Inf (Gam = 1.5)')
pause
semilogx(w,svw3i,w,hsvt1,w,hsvt2)
title('H-Inf W-Plane Actuator Design -- 1/W3 & Comp. Sens. Func.')
xlabel('Frequency - Rad/Sec')
ylabel('SV - db')
grid
text(0.002,-30,'H-Inf (Gam = 1) ---> H-Inf (Gam = 1.5)')
pause
semilogx(w,hsvtt1,w,hsvtt2)
title('H-Inf W-Plane Actuator Design -- Cost function Ty1u1')
xlabel('Frequency - Rad/Sec')
ylabel('SV - db')
grid
text(0.002,-10,'H-Inf (Gam = 1) ---> H-Inf (Gam = 1.5)')
pause
clc
disp('  ')
disp('  ')
disp('               << H-Inf Controller (Gam = 1.5) >>')
disp('    ')
disp('    Poles of Controller :')
polecp = eig(acp)
disp(' ')
disp('                                 (strike a key to continue ...)')
pause
disp(' ')
disp('    Transmission Zeros of the Controller :')
tzcp = tzero(acp,bcp,ccp,dcp)
disp(' ')
disp('                                 (strike a key to continue ...)')
pause
svcp = sigma(acp,bcp,ccp,dcp,1,w); svcp = 20*log10(svcp);
semilogx(w,svcp)
title('Final 8-State H-Inf Controller')
xlabel('Rad/Sec')
ylabel('SV (db)')
pause
clc
disp('    Poles of the Cost Function :')
poletyu = eig(acl)
disp('                                 (strike a key to continue ...)')
pause
%
% ----- End of ACTDEMO.M --- RYC/MGS %
