home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / MUDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  9.8 KB  |  280 lines
  1. %
  2. % ---------------------------------------------------------------
  3. %  MUDEMO.M is a script file that demonstrates the design of 
  4. %     1990 ACC benchmark problem using MU-SYNTHESIS.
  5. % ---------------------------------------------------------------
  6. %
  7. % R. Y. Chiang & M. G. Safonov 
  8. % Copyright (c) 1988 by the MathWorks, Inc.
  9. % All Rights Reserved.
  10.  
  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('    |   2    | maximize MSM | Ts ~= 15 sec  |     impulse      |');
  27. disp('    |        | for nominal  | for nominal   |                  |');
  28. disp('    ------------------------------------------------------------');
  29. disp('   ')
  30. disp('                                  (strike a key to continue ...)');
  31. pause
  32. disp(' ')
  33. disp(' ')
  34. spring = [1 0.5 2];
  35. %
  36. clc
  37. disp(' ')
  38. disp('  -----------------------------------------------------------')
  39. disp('    1. Set up the nominal system:');
  40. disp(' ');
  41. disp('       k = 1; m1 = 1; m2 = 1;');
  42. disp('       ag = [0      0     1     0');
  43. disp('             0      0     0     1');
  44. disp('          -k/m1   k/m1    0     0');
  45. disp('           k/m2  -k/m2    0     0];');
  46. disp('       bg = [0 0 1/m1 0]`; cg = [0 1 0 0]; dg = 0;');
  47. disp('  -----------------------------------------------------------')
  48. k = 1; m1 = 1; m2 = 1;
  49. %
  50. ag = [0      0     1     0
  51.       0      0     0     1
  52.    -k/m1   k/m1    0     0
  53.     k/m2  -k/m2    0     0];
  54. bg = [0 0 1/m1 0]';
  55. cg = [0 1 0 0];
  56. dg = 0;
  57. disp(' ');
  58. disp('       Poles of G(s) - -')
  59. lamg = eig(ag)
  60. disp(' ')
  61. disp('                       (strike a key to continue)')
  62. pause
  63. %---------------------------------------
  64. % Augment Plant for H-Inf Design:
  65. %
  66. clc
  67. disp(' ')
  68.    disp('  -----------------------------------------------------------');
  69.    disp('    2. Set up the augment plant for H-Inf design:');
  70.    disp(' ')
  71.    disp('       Two-port state-space:')
  72.    disp('          (A,BB1,BB2,CC1,CC2,DD11,DD12,DD21,DD22)')
  73.    disp(' ')
  74. w2flag = 0;
  75. %   w2flag = input('     Select W2: Dynamic (1); Constant (2): ');
  76. if w2flag ~= 0
  77.    disp('          Port 1: 3 uncertainties + W2 ')
  78.    disp('          Port 2: controller')
  79.    disp('  -----------------------------------------------------------');
  80.    disp(' ')
  81.    rho = input('     Assign the DC of W2 weighting for control energy: ');
  82.    if w2flag == 1
  83.       nuw2 = rho*[1/10 1]; dnw2 = [1/200 1];
  84.       [aw2,bw2,cw2,dw2] = tf2ss(nuw2,dnw2);
  85.    else
  86.       aw2 = []; bw2 = []; cw2 = []; dw2 = rho;
  87.    end
  88.    AA = diagmx(ag,aw2);
  89.    BB1 = [ 0     0    0     0;
  90.            0     0    0     0;
  91.          -1/m1 -1/m1  0    k/m1;
  92.           1/m2   0  -1/m2 -k/m2;
  93.           zeros(size(aw2)*[1;0],4)];
  94.    BB2 = [0;0;1/m1;0;bw2];
  95.    CC1 =  [ 1    -1   0 0;
  96.           -k/m1  k/m1 0 0;
  97.            k/m2 -k/m2 0 0;
  98.             0     0   0 0]; 
  99.    CC1 = [CC1 [zeros(3,size(cw2)*[0;1]);cw2]];
  100.    CC2 = [0 1 0 0 zeros(1,size(cw2)*[0;1])];
  101.    DD11 = [0 0 0 0;-1/m1 -1/m1 0 k/m1;1/m2 0 -1/m2 0; 0 0 0 0];
  102.    DD12 = [0;1/m1;0;dw2]; DD21 = [0 0 0 1]; DD22 = 0;
  103. else
  104.    disp('          Port 1: 3 uncertainties')
  105.    disp('          Port 2: controller')
  106.    disp('  -----------------------------------------------------------');
  107.    AA = ag; BB1 = [0 0 0;0 0 0;-1/m1 -1/m1 0;1/m2 0 -1/m2]; 
  108.    BB2 = [0 0 1/m1 0]'; 
  109.    CC1 = [1 -1 0 0;-k/m1 k/m1 0 0;k/m2 -k/m2 0 0]; CC2 = [0 1 0 0];
  110.    DD11 = [0 0 0;-1/m1 -1/m1 0;1/m2 0 -1/m2]; DD12 = [0 1/m1 0]'; 
  111.    DD21 = [0 0 0]; DD22 = 0;
  112. end   
  113. %
  114. disp(' ');
  115. %
  116. no_u1 = size(BB1)*[0;1]; no_u2 = size(BB2)*[0;1];
  117. no_y1 = size(CC1)*[1;0]; no_y2 = size(CC2)*[1;0];
  118. %
  119. BB = [BB1 BB2]; CC = [CC1;CC2]; DD = [DD11 DD12;DD21 DD22];
  120. %
  121. disp(' ');
  122. disp('       Poles of Augment Plant P(s) - -')
  123. lamp = eig(AA)
  124. disp('                (strike a key to continue)')
  125. pause
  126. % --------------------------------------
  127. % JW-axis shifting from s -----> s~:
  128. %
  129. clg
  130. bilexp
  131. clc
  132. disp(' ')
  133. disp(' --------------------------------------------------------------------')
  134. disp('    3. Transform the Augment Plant to Shifted JW-Axis:')
  135. disp(' ')
  136. disp('       % pack the 2-port state-space ...')
  137. disp('       BB = [BB1 BB2]; CC = [CC1;CC2]; DD = [DD11 DD12;DD21 DD22];')
  138. disp('       % select the circle point for mapping ...')
  139. disp('       [aa,bb,cc,dd] = bilin(AA,BB,CC,DD,1,`Sft_jw`,cirpt);')
  140. disp('       % split the (aa,bb,cc,dd) to 2-port state-space (A,B1,B2,..);')
  141. disp(' --------------------------------------------------------------------') 
  142. disp('      Input the circle point "p1" - - -')       
  143. cirpt1 = input('     (-0.25: non-min. phase controller, -0.3: min. phase controller): ');
  144. cirpt = [-100 cirpt1];
  145. [aa,bb,cc,dd] = bilin(AA,BB,CC,DD,1,'Sft_jw',cirpt);
  146. A = aa; B1 = bb(:,1:no_u1); B2 = bb(:,no_u1+1:no_u1+no_u2); 
  147. C1 = cc(1:no_y1,:); C2 = cc(no_y1+1:no_y1+no_y2,:);
  148. D11 = dd(1:no_y1,1:no_u1); D12 = dd(1:no_y1,no_u1+1:no_u1+no_u2); 
  149. D21 = dd(no_y1+1:no_y1+no_y2,1:no_u1); 
  150. D22 = dd(no_y1+1:no_y1+no_y2,no_u1+1:no_u1+no_u2);
  151. disp(' ')
  152. disp('       Poles of P(s~) - -')
  153. lampp = eig(A)
  154. disp(' ')
  155. disp('                    (strike a key to continue)')
  156. pause
  157. clc
  158. disp(' ')
  159. disp('  ------------------------------------------------------------------')
  160. disp('    4. Starting H-Inf Design in "s~" Domain:');
  161. disp(' ')
  162. disp('       TSS_ = mksys(A,B1,B2,C1,C2,D11,D12,D21,D22,`tss`);')
  163. disp('       % H-Inf GAMMA-Iteration for maximum MSM')
  164. disp('       [gamopt,ss_cp,ss_cl]=hinfopt(TSS_);')
  165. disp('       [acp,bcp,ccp,dcp] = branch(ss_cp);')
  166. disp('       [acl,bcl,ccl,dcl] = branch(ss_cl);')
  167. disp('  ------------------------------------------------------------------')
  168. disp('  ')
  169. disp('  ')
  170. disp('                    (strike a key to continue)')
  171. pause
  172. format short e
  173. [gamopt,acp,bcp,ccp,dcp,acl,bcl,ccl,dcl]=hinfopt(...
  174.                     A,B1,B2,C1,C2,D11,D12,D21,D22);
  175. disp(' ')
  176. disp(' ')
  177. disp('                    (strike a key to continue)')
  178. pause
  179. [nucp,dncp] = ss2tf(acp,bcp,ccp,dcp,1);
  180. w = logspace(-3,3,100);
  181. [acp,bcp,ccp,dcp] = tf2ss(nucp,dncp);
  182. % --------------------------------------------------------------------
  183. % JW-axis shifting on H-Inf cost & controller from s~------>s
  184. %
  185. clc
  186. disp(' ')
  187. disp('   ----------------------------------------------------------------------')
  188. disp('      Transform the cost & controller from "s~" to "s":')
  189. disp(' ')
  190. disp('      [acll,bcll,ccll,dcll] = bilin(acl,bcl,ccl,dcl,-1,`Sft_jw`,cirpt);')
  191. disp('      [af,bf,cf,df] = bilin(acp,bcp,ccp,dcp,-1,`Sft_jw`,cirpt);')
  192. disp('   ----------------------------------------------------------------------')
  193. disp(' ')
  194. disp(' ')
  195. disp('                     (strike a key to continue)')
  196. pause
  197. [acll,bcll,ccll,dcll] = bilin(acl,bcl,ccl,dcl,-1,'Sft_jw',cirpt);
  198. [af,bf,cf,df] = bilin(acp,bcp,ccp,dcp,-1,'Sft_jw',cirpt);   
  199. disp(' ')
  200. disp('     - - Computing the Cost Function Ty1u1(s~) and Ty1u1(s) - -');
  201. [perttw1,diag_d] = ssv(acl,bcl,ccl,dcl,w);
  202. msmw1 = gamopt/max(perttw1)*100;
  203. perttw1 = 20*log10(perttw1);
  204. [pertt1,diag_dl]  = ssv(acll,bcll,ccll,dcll,w);
  205. msm1 = gamopt/max(pertt1)*100;
  206. pertt1  = 20*log10(pertt1);  
  207. clg
  208. semilogx(w,[perttw1;pertt1]);
  209. title('Cost Function in "s~" and "s"');
  210. grid;xlabel('Rad/Sec'); ylabel('PERRON SSV (db)');
  211. text(0.2,0.3,['ADDITIVE MSM IN s~-DOMAIN: +-',num2str(msmw1),' %'],'sc');
  212. text(0.2,0.2,['ADDITIVE MSM IN s-DOMAIN : +-',num2str(msm1),' %'],'sc');
  213. %prtsc
  214. disp(' ')
  215. disp(' ')
  216. disp('                     (strike a key to continue)')
  217. pause
  218. clc
  219. disp(' ')
  220. disp('  ---------------------------------------------------------------')
  221. disp('    5. Start to Curve-Fit the Diagonal Scaling D Matrix:')
  222. disp('       (compute the cost function (acl,bcl,ccl,dcl) again)')
  223. disp('    ')
  224. disp('       w = logspace(-3,3,100);')
  225. disp('       [SS_D] = fitd(diag_d,w,[3 2 2]);')
  226. disp('     % Absorb the D(s) and inv(D(s)) into the augmented plant:')
  227. disp('       [TSS_D] = augd(TSS_,SS_D);')
  228. disp('  ----------------------------------------------------------------')
  229. disp(' ')
  230. disp(' ')
  231. disp('                     (strike a key to continue)')
  232. pause
  233. %
  234. if w2flag == 0
  235.    ord_d = [3 2 2];
  236. else
  237.    ord_d = [2  2 2 3];
  238. end
  239. [AD,BD,CD,DD] = fitd(diag_d,w,ord_d);
  240. [ga1,ph1] = bode(AD,BD(:,1),CD(1,:),DD(1,1),1,w);
  241. [ga2,ph2] = bode(AD,BD(:,2),CD(2,:),DD(2,2),1,w);
  242. [ga3,ph3] = bode(AD,BD(:,3),CD(3,:),DD(3,3),1,w);
  243. if w2flag ~= 0
  244.    [ga4,ph4] = bode(AD,BD(:,4),CD(4,:),DD(4,4));
  245. end
  246. SS_D = mksys(AD,BD,CD,DD);
  247. TSS_ = mksys(A,B1,B2,C1,C2,D11,D12,D21,D22,'tss');
  248. TSS_D = augd(TSS_,SS_D);
  249. clc
  250. disp(' ')
  251. disp('  -------------------------------------------------------------------')
  252. disp('    6. Design the H-Inf Controller for the Diagonally-Scaled Plant')
  253. disp(' ')
  254. disp('       [gamopt,ss_cp,ss_cl] = hinfopt(TSS_D);')
  255. disp('       [acp,bcp,ccp,dcp] = branch(ss_cp);')
  256. disp('       [acl,bcl,ccl,dcl] = branch(ss_cl);')
  257. disp('  -------------------------------------------------------------------')
  258. disp(' ')
  259. disp(' ')
  260. disp('                    (strike a key to continue)')
  261. pause
  262. [gamopt,ss_cp,ss_cl] = hinfopt(TSS_D);
  263. [acp,bcp,ccp,dcp] = branch(ss_cp);
  264. [acl,bcl,ccl,dcl] = branch(ss_cl);
  265. disp(' ')
  266. disp(' ')
  267. disp('                    (strike a key to continue)')
  268. pause
  269. %
  270. mueva
  271. muplt
  272. %
  273. % ------------ End of MUDEMO.M % RYC/MGS %
  274.  
  275.  
  276.  
  277.  
  278.  
  279.  
  280.