home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / DINTDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  5.9 KB  |  148 lines
  1. % DINTDEMO.M --- A demo of hinf design for a double integrator plant
  2. %
  3.  
  4. % R. Y. Chiang & M. G. Safonov 
  5. % Copyright (c) 1988 by the MathWorks, Inc.
  6. % All Rights Reserved.
  7.  
  8. clc
  9. disp(' ')
  10. disp(' ----------------------------------------------------------------')
  11. disp('         << H-INF DESIGN FOR A DOUBLE INTEGRATOR PLANT >>')
  12. disp(' ')
  13. disp('   Model: a double integrator with moment inertia J (1/J/s/s)')
  14. disp(' ')
  15. disp('   For example: rigid body dynamics of a spacecraft (J = 5700)')
  16. disp(' ')
  17. disp('                           * * * * *')
  18. disp('      * * * * * * * *      * space *      * * * * * * * *')
  19. disp('      * solar panel * ---- * craft * ---- * solar panel *')
  20. disp('      * * * * * * * *      *1/j/s/s*      * * * * * * * *')
  21. disp('                           * * * * *')
  22. disp('          ')
  23. disp('   Objective: Find a stabilizing Mixed-Sensitivity H-Inf')
  24. disp('              controller.')
  25. disp(' ----------------------------------------------------------------')
  26. disp('   ')
  27. disp('                                  (strike a key to continue ...)');
  28. pause
  29. den = conv([1 0],[1 0]);
  30. [ag,bg,cg,dg] = tf2ss(1/5700,den);
  31. clc
  32. disp(' ' )
  33. disp(' ----------------------------------------------------------------')
  34. disp('                 << JW-AXIS POLES/ZEROS >>')
  35. disp(' ')
  36. disp('  One special feature of state-space MIXED-SENSITIVITY approach:')
  37. disp('   ')
  38. disp('     PLANT DOES NOT ALLOWED TO HAVE JW-AXIS POLE/ZERO !!')
  39. disp('   ')
  40. disp('  But this does not mean we can not deal with the situation.')
  41. disp('  Simply use the jw-axis shifting technique to shift the poles:')
  42. disp(' ')
  43. disp('    Step 1:  Shift the plant (Ag <----- Ag + SIGMA * eye(Ag))')
  44. disp('    Step 2:  Do a standard mixed-sensitivity H-Inf design')
  45. disp('    Step 3:  Shift back the controller ')
  46. disp('             (Acp <------ Acp - SIGMA * eye(Acp))')
  47. disp('                       Done !!')
  48. disp('  where SIGMA can be a small positive number (e.g., 0.1)')
  49. disp(' ---------------------------------------------------------------')
  50. disp('   ')
  51. disp('                                  (strike a key to continue ...)');
  52. pause
  53. SIGMA = 0.1;
  54. ag0 = ag + SIGMA*eye(2);
  55. disp('   ')
  56. disp(' Poles of the plant/shifted plant:')
  57. [eig(ag) eig(ag0)]
  58. disp('   ')
  59. disp('                                  (strike a key to continue ...)');
  60. pause
  61. clc
  62. disp(' ')
  63. disp(' --------------------------------------------------------------------')
  64. disp('               << W3 WEIGHTING & JW-AXIS ZEROS >>')
  65. disp('')
  66. disp('  We have delt with the jw-axis plant POLES by shifting the axis.')
  67. disp('  But there are still two plant ZEROS at infinity, which are also on ')
  68. disp('  the jw-axis ....')
  69. disp('  We can solve this problem via a clever W3 weighting:')
  70. disp('  ')
  71. disp('                        W3 = s^2 / 100')
  72. disp('  ')
  73. disp('  where the double differentiator makes the plant full rank at  ')
  74. disp('  infinity, but also serves as the compelmentary sensitivity')
  75. disp('  wighting function to control the system bandwidth.')
  76. disp('  In our example, system bandwidth is set to be 10 rad/sec.')
  77. disp(' ---------------------------------------------------------------------')
  78. disp('   ')
  79. disp('                                  (strike a key to continue ...)');
  80. pause
  81. clc
  82. disp(' ')
  83. disp(' ---------------------------------------------------------------------')
  84. disp('       << SENSITIVITY WEIGHTING (W1), THE H-INF DESIGN KNOB >>')
  85. disp(' ')
  86. disp('  By tuning W1 filter, we can find an H-Inf controller:')
  87. disp('                          2                                2')
  88. disp('            beta * [alfa*S + 2*zeta1*w1c*sqrt(alfa)*S + w1c ]')
  89. disp('      W1 = ----------------------------------------------------')
  90. disp('                          2                                2')
  91. disp('                   [beta*S + 2*zeta2*w1c*sqrt(beta)*S + w1c ]')
  92. disp('   ')
  93. disp('  where ')
  94. disp('      beta: DC gain of the filter (controls the disturbance rejection)')
  95. disp('      alfa: high frequency gain (controls the peak overshoot)')
  96. disp('      w1c: filter cross-over frequency')
  97. disp('      zeta1, zeta2: damping ratios of the corner frequencies.')
  98. disp(' ')
  99. disp('  Inverse of W1 is the desired shape of the sensitivity function.')
  100. disp('  Here, we choose 2nd order W1 filter for a uniform loop-shaping on ')
  101. disp('  L(s) = G(s)*F(s).')
  102. disp(' --------------------------------------------------------------------')
  103. disp('   ')
  104. disp('                                  (strike a key to continue ...)');
  105. pause
  106. w1c = input('Input the sensitivity cross-over frequqncy (try 3): ');
  107. beta = 100;
  108. alfa = input('Input the sensitivity upper bound "ALFA" (try 1.5): ');
  109. alfa = 1/alfa;
  110. zeta1 = 0.7; zeta2 = 0.7;
  111. w1 = [beta*[alfa 2*zeta1*w1c*sqrt(alfa) w1c*w1c];...
  112.            [beta 2*zeta2*w1c*sqrt(beta) w1c*w1c]];
  113. w2 = [];
  114. w3 = [1 0 0;0 0 100];
  115. clc
  116. disp(' -----------------------------------------------------------')
  117. disp('                   << H-INF DESIGN >>')
  118. disp(' ')
  119. disp('   >> w1 = [beta*[alfa 2*zeta1*w1c*sqrt(alfa) w1c*w1c];..')
  120. disp('                 [beta 2*zeta2*w1c*sqrt(beta) w1c*w1c]];')
  121. disp('   >> w2 = [];')
  122. disp('   >> w3 = [1 0 0;0 0 100];')
  123. disp('   >> ss_g = mksys(ag0,bg,cg,dg);')
  124. disp('   >> TSS_ = augtf(ss_g,w1,w2,w3);')
  125. disp('   >> [ss_cp,ss_cl,hinfo] = hinf(TSS_);')
  126. disp('   >> [acp,bcp,ccp,dcp] = branch(ss_cp);')
  127. disp('   >> [acl,bcl,ccl,dcl] = branch(ss_cl);')
  128. disp(' -----------------------------------------------------------')
  129. disp('   ')
  130. disp('                                  (strike a key to continue ...)');
  131. pause
  132. disp(' ')
  133. disp(' - - - Working Plant Augmentation - - - Wait - - -')
  134. ss_g = mksys(ag0,bg,cg,dg);
  135. TSS_ = augtf(ss_g,w1,w2,w3); 
  136. clc
  137. [ss_cp,ss_cl,hinfo] = hinf(TSS_);
  138. [acp,bcp,ccp,dcp] = branch(ss_cp);
  139. [acl,bcl,ccl,dcl] = branch(ss_cl);
  140. acp = acp-SIGMA*eye(acp);
  141. if max(real(eig(acl))) < 0
  142.   dinteva
  143.   disp(' ')
  144.   dintplt
  145. end
  146. %
  147. % ---------- End of DINTDEMO.M % RYC/MGS
  148.