home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / MUEVA.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  5.3 KB  |  149 lines
  1. %
  2. % ---------------------------------------------------------------
  3. %  MUEVA.M is a script file that evaluates the performance of
  4. %     the 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.  
  12. clc
  13. disp(' ');
  14. disp('          << Start to Evaluate the Robust Performance >>');
  15. %
  16. %w = [w [0.1:0.1:2]];
  17. %[w,wind] = sort(w);
  18. %
  19. % --------------------------------------------------------------------
  20. % JW-axis shifting on H-Inf cost & controller from s~------>s
  21. %
  22. disp(' ')
  23. disp('   ----------------------------------------------------------------------')
  24. disp('      Transform the cost & controller from "s~" to "s":')
  25. disp(' ')
  26. disp('      [acll,bcll,ccll,dcll] = bilin(acl,bcl,ccl,dcl,-1,`Sft_jw`,cirpt);')
  27. disp('      [aff,bff,cff,dff] = bilin(acp,bcp,ccp,dcp,-1,`Sft_jw`,cirpt);')
  28. disp('      Balancing the controller for better numerical property:')
  29. disp('      [af,bf,cf,hsv] = obalreal(aff,bff,cff); df = dff;')
  30. disp('   ----------------------------------------------------------------------')
  31. disp(' ')
  32. disp(' ')
  33. disp('                     (strike a key to continue)')
  34. pause
  35. [acll,bcll,ccll,dcll] = bilin(acl,bcl,ccl,dcl,-1,'Sft_jw',cirpt);
  36. [aff,bff,cff,dff] = bilin(acp,bcp,ccp,dcp,-1,'Sft_jw',cirpt);   
  37. [af,bf,cf,hsv] = obalreal(aff,bff,cff); df = dff;
  38. clc
  39. disp(' ')
  40. disp('     - - Computing the cost function Ty1u1(s~) and Ty1u1(s) - -');
  41. [perttw,diag_dl] = ssv(acl,bcl,ccl,dcl,w); 
  42. msmw = gamopt/max(perttw)*100;
  43. perttw = 20*log10(perttw);
  44. [pertt,diag_dl] = ssv(acll,bcll,ccll,dcll,w); 
  45. msm = gamopt/max(pertt)*100;
  46. pertt = 20*log10(pertt);  
  47. disp(' ')
  48. disp('       - - Evaluating Controller - - ')
  49. disp(' ')
  50. lamf = eig(af);
  51. stabflag = size(find(real(lamf)>0));
  52. if stabflag*[1;0] ~= 0
  53.         disp('                Controller is unstable !!');
  54. else
  55.         disp('                Controller is stable - - -');
  56. end
  57. disp(' ')
  58. tzf = tzero(af,bf,cf,df);
  59. stabflag = size(find(real(tzf)>0));
  60. if stabflag*[1;0] ~= 0
  61.            disp('                Controller is non-minimum phase !!');
  62. else
  63.            disp('                Controller is minimum phase - - -');
  64. end
  65. disp(' ')
  66. disp('                             (strike a key to continue)')
  67. pause
  68. clc
  69. disp(' ')
  70. disp(' ')
  71. disp('      - - Computing the gain/phase of F(s) & L(s) - -');
  72. [gf,pf] = bode(af,bf,-cf,-df,1,w); gf = 20*log10(gf);
  73. % --------------------------------------------------------------
  74. % Evaluate the disturbance response from u,w ---> z 
  75. %
  76. disp(' ')
  77. disp('      - - Evaluating the disturbance response from w to z - -')
  78. BB1 = BB1(:,1); DD21 = 0; DD22 = 0; 
  79. % ---------------------------------------------------------------
  80. % Including the control energy, disturbance at M1 or M2:
  81. %
  82. BB1a = [0 0 0 1/m2]';         % disturbance at M2
  83. BB1a = [BB1a,[0 0 0 0]'];     % add V(s) 
  84. BB1b = [0 0 1/m1 0]';         % disturbance at M1
  85. BB1b = [BB1b,[0 0 0 0]'];     % add V(s)
  86. CC2a = [1 0 0 0;0 1 0 0;0 0 0 0]; 
  87. DD21a = [0 0;0 0;0 0]; DD22a = [0;0;1];
  88. DD21  = [0 1]; DD22 = 0;
  89. t = 0:0.1:30;
  90. dist_w = sin(0.5*t)';
  91. disp(' ')
  92. disp('      - - Inject the sensor noise: v(t) = 0.001*sin(100*t) - -');
  93. disp(' ')
  94. nos_v  = 0.001*sin(100*t)';
  95. [tmp1, tmp2]=size(t);
  96. ipp = zeros(tmp1, tmp2)'; ipp(1) = 2/(t(2)-t(1));
  97. U1 = [ipp nos_v];          % impulse + sensor noise
  98. % -----------------------------------------------------------------------
  99. %
  100. no_spr = size(spring)*[0;1];
  101. for k0no = 1:no_spr
  102.       k0 = spring(k0no)
  103.       AA = [   0       0     1     0
  104.                0       0     0     1
  105.           -k0/m1   k0/m1     0     0
  106.            k0/m2  -k0/m2     0     0];
  107.       % dist. at M2
  108.       [aw2zu,bw2zu,cw2zu,dw2zu] = lftf(...
  109.             AA,BB1a,BB2,CC2a,CC2,DD21a,DD22a,DD21,DD22,af,bf,cf,df);
  110.       lamw2zu = eig(aw2zu);
  111.       stabflag = size(find(real(lamw2zu)>0));
  112.       if stabflag*[1;0] ~= 0
  113.          disp('      Closed-Loop unstable - - -');
  114.       else
  115.          disp('      Closed-Loop stable - - -');
  116.       end
  117.       % dist. at M1
  118.       [aw1zu,bw1zu,cw1zu,dw1zu] = lftf(...
  119.             AA,BB1b,BB2,CC2a,CC2,DD21a,DD22a,DD21,DD22,af,bf,cf,df);
  120.       imp_w2 = lsim(aw2zu,bw2zu,cw2zu,dw2zu,U1,t);
  121.       x1_ipw2(:,k0no) = imp_w2(:,1);
  122.       x2_ipw2(:,k0no) = imp_w2(:,2);
  123.       u_ipw2(:,k0no)  = imp_w2(:,3);
  124.       imp_w1 = lsim(aw1zu,bw1zu,cw1zu,dw1zu,U1,t);
  125.       x1_ipw1(:,k0no) = imp_w1(:,1);
  126.       x2_ipw1(:,k0no) = imp_w1(:,2);
  127.       u_ipw1(:,k0no)  = imp_w1(:,3);
  128.       svw2z(:,k0no) = bode(aw2zu,bw2zu,cw2zu(2,:),dw2zu(2,:),1,w);
  129.       svw2z(:,k0no) = 20*log10(svw2z(:,k0no));
  130.       disp(' ')
  131.       ag = [0      0     1     0
  132.             0      0     0     1
  133.        -k0/m1   k0/m1    0     0
  134.         k0/m2  -k0/m2    0     0]; 
  135.       bg = [0 0 1/m1 0]'; 
  136.       cg = [0 1 0 0]; dg = 0;
  137.       [al,bl,cl,dl] = series(af,bf,-cf,-df,ag,bg,cg,dg);
  138.       [nul(k0no,:),dnl(k0no,:)] = ss2tf(al,bl,cl,dl,1);
  139.       [gl(:,k0no),pl(:,k0no)] = bode(al,bl,cl,dl,1,w); 
  140.       [gm(k0no,1),pm(k0no,1),wg(k0no,1),wp(k0no,1)] = ...
  141.                      margin(gl(:,k0no),pl(:,k0no),w);
  142.       gm = 20*log10(gm); 
  143.       gl(:,k0no) = 20*log10(gl(:,k0no));
  144. end     % of spring constant loop
  145. %
  146. gmin = min(gm); pmin = min(pm);
  147. gmax = max(gm); pmax = max(pm);
  148. %
  149.