home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / ACCDEMO.M next >
Encoding:
Text File  |  1993-03-11  |  10.2 KB  |  277 lines
  1. %
  2. % ---------------------------------------------------------------
  3. %  ACCDEMO.M is a script file that demonstrates the designs of 
  4. %     1990 ACC benchmark problem.
  5. % ---------------------------------------------------------------
  6. %
  7.  
  8. % R. Y. Chiang & M. G. Safonov 
  9. % Copyright (c) 1988 by the MathWorks, Inc.
  10. % All Rights Reserved.
  11. clc
  12. clear
  13. disp(' ')
  14. disp('                      << Benchmark Problem >>')
  15. disp('  ');
  16. disp('                           |--> x1        |--> x2 = z (measurement)');
  17. disp('         (control)    ------     k   ------');
  18. disp('           u -------> | M1 | -/\/\/\-| M2 | ------> w (disturbance)');
  19. disp('                      ------         ------');
  20. disp('                       0  0           0  0');
  21. disp('   ===============================================================');
  22. disp('                     [ Design Specifications ]');
  23. disp('    ------------------------------------------------------------');
  24. disp('    | Design | Robustness   | Settling Time | Disturbance w(t) |')
  25. disp('    ---------+--------------+---------------+-------------------');
  26. disp('    |   1    | stable for   | Ts ~= 15 sec  |     impulse      |');
  27. disp('    |        | .5 < k < 2   | for nominal   |                  |');
  28. disp('    ---------+--------------+---------------+-------------------');
  29. disp('    |   2    | maximize MSM | Ts ~= 15 sec  |     impulse      |');
  30. disp('    |        | for nominal  | for nominal   |                  |');
  31. disp('    ---------+--------------+---------------+-------------------');
  32. disp('    |   3    | stable for   | Ts ~= 20 sec  | A sin(.5t + phi) |');
  33. disp('    |        | .5 < k < 2   | for all cases | A, phi unknown   |');
  34. disp('    ------------------------------------------------------------');
  35. disp('   ')
  36. disp(' ')
  37. flag = input('Design # 1, 2 or 3 : ');
  38. Lag = 0;   % Lag = 1; for design # 1, non-min. phase controller
  39. if flag == 1
  40.    Lag = input('  Controller type (0 = min phase, 1 = non-min phase): ')
  41. end;
  42. spring = [1 0.5 2];
  43. %
  44. clc
  45. disp(' ')
  46. disp('  -----------------------------------------------------------')
  47. disp('    1. Set up the nominal system:');
  48. disp(' ');
  49. disp('       k = 1; m1 = 1; m2 = 1;');
  50. disp('       ag = [0      0     1     0');
  51. disp('             0      0     0     1');
  52. disp('          -k/m1   k/m1    0     0');
  53. disp('           k/m2  -k/m2    0     0];');
  54. disp('       bg = [0 0 1/m1 0]`; cg = [0 1 0 0]; dg = 0;');
  55. disp('  -----------------------------------------------------------')
  56. k = 1; m1 = 1; m2 = 1;
  57. ag = [0      0     1     0
  58.       0      0     0     1
  59.    -k/m1   k/m1    0     0
  60.     k/m2  -k/m2    0     0];
  61. bg = [0 0 1/m1 0]'; cg = [0 1 0 0]; dg = 0;
  62. disp(' ');
  63. disp('       Poles of G(s) - -')
  64. lamg = eig(ag)
  65. disp(' ')
  66. disp('                       (strike a key to continue)')
  67. pause
  68. %---------------------------------------
  69. % Augment Plant for H-Inf Design:
  70. %
  71. clc
  72. if Lag == 0
  73.   if flag == 1 | flag == 3 
  74.       disp(' ')
  75.       disp('  -----------------------------------------------------------');
  76.       disp('    2. Set up the augment plant for H-Inf design:');
  77.       disp(' ')
  78.       disp('       Change the spring to "5/4" with scaling factor "Gam"')
  79.       disp('       (the additive uncertainty is now bounded by +- 1)')
  80.       disp('       k0 = 5/4;')
  81.       disp('       A = [   0       0     1     0')
  82.       disp('               0       0     0     1')
  83.       disp('          -k0/m1   k0/m1     0     0')
  84.       disp('           k0/m2  -k0/m2     0     0];')
  85.       disp('       B1 = [0 0 -1/m1  1/m2]`; B2 = -[0 0 1/m1 0]`;')
  86.       disp('       C1 = Gam*[1 -1 0 0];   C2 = [0 1 0 0];')
  87.       disp('       D11 = 0; D12 = 0; D21 = 0; D22 = 0;')
  88.       disp('       where "Gam" is the robustness level of spring')
  89.       disp('       (Gam = 0.75 for 0.5 <= k <= 2).')
  90.       disp('  ----------------------------------------------------------');
  91.       disp('  ')
  92. %      itcha = input('Assign iteration channel: (1=robustness, 2=control, 4=fix no. on each) ');
  93.       itcha = 4;
  94.       if itcha == 4
  95.          disp(' ')
  96.          if flag == 1
  97.             disp(' Try:  Gam = 0.75, Rho = 0.02...')
  98.          end
  99.          if flag == 2
  100.             disp(' Try:  Gam = 0.75, Rho = 0...')
  101.          end
  102.          if flag == 3
  103.             disp(' Try:  Gam = 0.4,  Rho = 0.01...')
  104.          end
  105.          disp(' ')
  106.          Gam = input('Assign the fix robustness level "Gam" of spring constant: ');
  107.          Rho = input('Assign the fix bound "Rho" on control signal: '); 
  108.       end
  109.       if itcha == 2 
  110.          Gam = input('Assign the fix robustness level "Gam" of spring constant: ');
  111.          Rho = 1;
  112.       end
  113.       if itcha == 1 
  114.          Rho = input('Assign the fix bound "Rho" on control signal: ');
  115.          Gam = 1;
  116.       end
  117.       k0 = 1.25;
  118.       A  = [   0       0     1     0
  119.                0       0     0     1
  120.           -k0/m1   k0/m1     0     0
  121.            k0/m2  -k0/m2     0     0]; 
  122.       B1 = [0 0 -1/m1  1/m2]'; B2 = -[0 0 1/m1 0]';
  123.       C1 = Gam*[1 -1 0 0];   C2 = [0 1 0 0];
  124.       D11 = 0; D12 = 0; D21 = 0; D22 = 0;
  125.       if Rho ~= 0
  126.           C1 = [C1;0 0 0 0]; D11 = [D11;0]; D12 = [D12;Rho];
  127.       end
  128.   end
  129. end
  130. if flag == 3
  131.    clc
  132.    disp(' ')
  133.    disp('  ----------------------------------------------------------');
  134.    disp('       Augment the plant with Internal Model:')
  135.    disp('                         (s+1)^2')
  136.    disp('                 I(s) = -----------')
  137.    disp('                        s^2 + 0.5^2')
  138.    disp('  ----------------------------------------------------------');   
  139.    [aint,bint,cint,dint] = tf2ss([1 2 1],[1 0 0.25]);
  140.   if Lag == 0
  141.    [aa,bb,cc,dd] = series(aint,bint,cint,dint,A,B2,[C1;C2],[D12;D22]);
  142.    A = aa; B1 = [zeros(size(aint)*[1;0],1);B1]; B2 = bb;
  143.    C1 = cc(1:2,:); C2 = cc(3,:); D12 = dd(1:2,1); D22 = dd(3,1);
  144.   else
  145.    [ag,bg,cg,dg] = series(aint,bint,cint,dint,ag,bg,cg,dg);
  146.   end
  147.    disp('  ')
  148.    disp('  ')
  149.    disp('                    (strike a key to continue)')
  150.    pause   
  151. end
  152. if flag == 2 
  153.    disp(' ')
  154.    disp('  -----------------------------------------------------------');
  155.    disp('    2. Set up the augment plant for H-Inf design:');
  156.    disp(' ');
  157.    disp('       A = ag; ');
  158.    disp('       B1 = [0 0 0;0 0 0;-1/m1 1 0; 1/m2 0 1];')
  159.    disp('       B2 = -[0 0 1/m1 0]`; ')
  160.    disp('       C1 = [1 -1 0 0;-k k 0 0;k -k 0 0]; C2 = [0 1 0 0];')
  161.    disp('       D11 = [0 0 0;-1 0 0;1 0 0]; D12 = -[0 1 0]`;')
  162.    disp('       D21 = [0 0 0]; D22 = 0;')
  163.    disp('  ----------------------------------------------------------'); 
  164.    A = ag; B1 = [0 0 0;0 0 0;-1/m1 1 0; 1/m2 0 1]; B2 = -[0 0 1/m1 0]'; 
  165.    C1 = [1 -1 0 0;-k k 0 0;k -k 0 0]; C2 = [0 1 0 0];
  166.    D11 = [0 0 0;-1 0 0;1 0 0]; D12 = -[0 1 0]'; 
  167.    D21 = [0 0 0]; D22 = 0;
  168.    itcha = 1:3;   % iterate on 3 robustness channels
  169.    disp('  ')
  170.    disp('  ')
  171.    disp('                    (strike a key to continue)')
  172.    pause
  173. end
  174. %
  175. clc
  176. disp(' ')
  177. disp(' ------------------------------------------------------------------')
  178. disp('    3. Transform the Plant via Pole-Shifting Bilinear Transform:')
  179. disp(' ')
  180. disp('       % select the circle point for mapping ...')
  181. disp('       % pack the 2-port state-space ...')
  182. disp('       B = [B1 B2]; C = [C1;C2]; D = [D11 D12;D21 D22];')
  183. disp('       [aa,bb,cc,dd] = bilin(A,B,C,D,1,`Sft_jw`,cirpt);')
  184. disp('       % split the (aa,..) back to 2-port state-space (A,B1,..);')
  185. disp('       % The H-Inf problem is ready to go ..')
  186. disp(' ------------------------------------------------------------------')
  187. disp(' ');
  188. disp(' ');
  189. if Lag == 0
  190.    disp('                  (strike a key to see an example of bilinear mapping)')
  191.    pause
  192.    bilexp
  193. else
  194.    disp('                  (strike a key to continue..)')
  195.    pause
  196. end
  197. disp(' ');
  198. disp(' ');
  199. disp('A class of H-Inf controllers can be found by iterating the circle point "p1":')
  200. disp(' ')
  201. disp(' ')
  202. if flag == 1
  203.    if Lag == 1
  204.       disp('            (p2 = INF; Try p1 = -0.3)');
  205.    else
  206.       disp('            (Min. Phase Controller: p1 = -0.35)');
  207.    end
  208. end
  209. if flag == 2
  210.       disp('            (Min. Phase Controller: p1 = -0.465)'); 
  211.       disp('            (Non-Min. Phase Controller: p1 = -0.25)');
  212. end
  213. if flag == 3
  214.       disp('            (Min. Phase Controller: p1 = -0.3)');
  215. end
  216. cirpt1 = input('   Input the circle point "p1": ');
  217. cirpt = [-100 cirpt1];
  218. %
  219. if Lag == 1
  220.    sgm = -cirpt1; itcha = 4;
  221.    [tmp1,tmp2]=size(ag);
  222.    ag = ag + sgm*eye(tmp1,tmp2);
  223.    disp(' ')
  224.    disp('Only impose W2 weighting on control signal - -')
  225.    Rho = input('   Assign the fix bound on W2 weight (try 0.1): '); 
  226.    w1 = []; w3 = []; w2 = [Rho;1];
  227.    [A,B1,B2,C1,C2,D11,D12,D21,D22] = augtf(ag,bg,cg,dg,w1,w2,w3);
  228. end
  229. no_u1 = size(B1)*[0;1]; no_u2 = size(B2)*[0;1];
  230. no_y1 = size(C1)*[1;0]; no_y2 = size(C2)*[1;0];
  231. if Lag == 0 
  232.    B = [B1 B2]; C = [C1;C2]; D = [D11 D12;D21 D22];
  233.    [aa,bb,cc,dd] = bilin(A,B,C,D,1,'Sft_jw',cirpt);
  234.    A = aa; B1 = bb(:,1:no_u1); B2 = bb(:,no_u1+1:no_u1+no_u2); 
  235.    C1 = cc(1:no_y1,:); C2 = cc(no_y1+1:no_y1+no_y2,:);
  236.    D11 = dd(1:no_y1,1:no_u1); D12 = dd(1:no_y1,no_u1+1:no_u1+no_u2); 
  237.    D21 = dd(no_y1+1:no_y1+no_y2,1:no_u1); 
  238.    D22 = dd(no_y1+1:no_y1+no_y2,no_u1+1:no_u1+no_u2);
  239. end
  240. %
  241. clc
  242. disp(' ')
  243. disp('  -------------------------------------------------------------')
  244. disp('    4. Starting H-Inf Design: ');
  245. disp(' ')
  246. disp('       TSS_ = mksys(A,B1,B2,C1,C2,D11,D12,D21,D22,`tss`);')
  247. if itcha == 4
  248.       disp('       % Regular H-Inf ')
  249.       disp('       [ss_cp,ss_cl,hinfo]=hinf(TSS_);')
  250. else
  251.       disp('       % H-Inf GAMMA-Iteration for maximum MSM')
  252.       disp('       [gamopt,ss_cp,ss_cl]=hinfopt(TSS_);')
  253. end
  254. disp('       [acp,bcp,ccp,dcp] = branch(ss_cp);')
  255. disp('       [acl,bcl,ccl,dcl] = branch(ss_cl);')
  256. disp('  -------------------------------------------------------------')
  257. disp('  ')
  258. disp('  ')
  259. disp('                    (strike a key to continue)')
  260. pause
  261. format short e
  262. clc
  263. if itcha == 4
  264.       [acp,bcp,ccp,dcp,acl,bcl,ccl,dcl,hinfo] = hinf(...
  265.                               A,B1,B2,C1,C2,D11,D12,D21,D22);
  266. else
  267.       [gamopt,acp,bcp,ccp,dcp,acl,bcl,ccl,dcl] = hinfopt(...
  268.                               A,B1,B2,C1,C2,D11,D12,D21,D22,itcha);
  269. end
  270. disp('  ')
  271. disp('                    (strike a key to continue)')
  272. pause
  273. %
  274. acceva
  275. accplt
  276. %
  277. % ------------ End of ACCDEMO.M % RYC/MGS %