home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / HIMAT_X1.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  11.7 KB  |  332 lines
  1.  echo on
  2. %
  3. %    himat_x1
  4. %
  5. %  This MATLAB script file is an example of how to design H-infinity 
  6. %   and mu-synthesis control laws, as applied to pitch axis control 
  7. %   for HIMAT. The purpose of this MATLAB script is to familiarize the
  8. %   user with the MATLAB function files developed for control design
  9. %   in this toolbox.
  10. %
  11. %   These control designs will be compared with a loop-shaping 
  12. %   control design using singular values, mu-analysis tools, the
  13. %   output sensitivity, input complementary sensitivity
  14. %   transfer functions, and time domain responses to step functions.
  15. %   
  16. %
  17. %   The files are set up as follows:
  18.  
  19.  pause    % strike any key to continue
  20. %                                           
  21. %                                              
  22. %   himat_x1    setup file for H-infinity control design, first 
  23. %               mu-synthesis iteration.
  24. %   himat_x2    second iteration of mu-synthesis using d-scales 
  25. %               from himat_x1.
  26. %   himat_x3    third iteration of mu-synthesis using d-scales 
  27. %               from himat_x2.
  28. %   himat_x4     fourth and final iteration of mu-synthesis
  29. %   himat_x5    comparing mu-synthesis control law with loop-shaping 
  30. %               design using singular values and mu.
  31. %   himat_x6    time domain analysis of both control laws.
  32.  
  33.  pause   % strike any key to continue
  34. clc
  35.  
  36. % First load the state space description of the HIMAT plane.
  37. %  This matrix is described in the manual. The states of the
  38. %  plant model are: forward speed (v), angle-of-attack (alpha),
  39. %  pitch rate (q) and pitch angle (theta). The inputs 
  40. %  are the elevon position and the canard position. The outputs 
  41. %  that are to be kept small are angle-of-attack (alpha) and pitch 
  42. %  angle (theta).
  43.  
  44.  pause   % strike any key to continue
  45.  mkhimat;
  46.  minfo(himat)
  47.  
  48. %    Next construct the weights associated with the multiplicative
  49. %      input uncertainty, using the command ND2SYS
  50. % 
  51. %                             50(s+100)
  52. %                   wdel  =   ---------
  53. %                             (s+10000)
  54.  
  55. wdel = nd2sys([50,5000],[1,10000]);
  56.  
  57. pause % strike any key to continue
  58. clc
  59.  
  60. %
  61. %    and the output error weight (performance),
  62. % 
  63. %             .5*(s+3)
  64. %       wp  =   --------
  65. %           (s+.03)
  66.  
  67. wp = nd2sys([0.5,1.5],[1,0.03]);
  68.  
  69. pause % strike any key to continue
  70. clc
  71.  
  72. %  Do a frequency_response of each weight and plot them on
  73. %   the same graph for comparison purposes. Evaluate each transfer
  74. %   function from 1e-4 to 1e+4 for a total of 40 points.
  75. %   On the graph the solid line is the uncertainty weight and
  76. %   the dashed line is the performance weight.
  77.  
  78.  om1 = logspace(-3,3,40);     % frequency vector
  79.  wdelg = frsp(wdel,om1);      % frsp of disturbance weight
  80.  wpg = frsp(wp,om1);          % frsp of performance weight
  81.  wghts = sbs(wdelg,wpg);      % place responses "side-by-side"
  82.  clg;
  83.  vplot('liv,lm',wghts);
  84.  echo off;
  85.  title('Uncertainty Weight and Performance Weight for HIMAT');
  86.  xlabel('Frequency (rad/s)');
  87.  ylabel('Magnitude');
  88.  pause   % strike any key to continue
  89.  clear wdelg wpg wghts
  90.  echo on;
  91.  pause   % strike any key to continue
  92.  
  93. %  Construct a two input/two-output weight for both the 
  94. %   uncertainty weight, WDEL and the output error weight,
  95. %   WP. These weights are associated with each channel 
  96. %   of the multivariable system. Each weight has 2 states, 
  97. %   2 inputs and 2 outputs as one can see using MINFO.
  98.  
  99.  wdel = daug(wdel,wdel);
  100.  wp = daug(wp,wp);
  101.  minfo(wdel)
  102.  minfo(wp)
  103.  pause   % strike any key to continue
  104.  
  105. %  Now that the plant and the weights have been defined, 
  106. %   construct the interconnection structure, HIMAT_IC, 
  107. %   based on the diagram provided. A MATLAB function file, 
  108. %   HIMATIC, has been written to do this for you.
  109. %   SYSIC is used to formulate the interconnection structure.
  110. %   HIMAT_IC has 8 states, 6 inputs and 6 outputs. The 
  111. %   first two inputs and outputs [1 2] correspond to the 
  112. %   multiplicative input uncertainty block, 
  113. %   the second two inputs and outputs [3 4] correspond to
  114. %   the disturbance inputs and error outputs. Inputs and 
  115. %   outputs [5 6] are the measurements and controls being 
  116. %   sent to and received from the controller.
  117.  
  118.  himatic
  119.  minfo(himat_ic)
  120.  clear wp wdel 
  121.  pause   % strike any key to continue
  122.  
  123. % The next step is to design a H-infinity control law for HIMAT.
  124. %  The function HINF designs a H-infinity control law based
  125. %  on the interconnection structure provided. HINF requires 
  126. %  the system interconnection structure (p), number of 
  127. %  measurements (nmeas), number of controls (ncon), a lower bound on 
  128. %  gamma (gmin), an upper bound on gamma (gmax), a tolerance for 
  129. %  the  gamma iteration (tol)  and which method to use in solving 
  130. %  the Riccati equations (ricmethd). Optional input arguments are the
  131. %  epsilon for Hamiltionian jw-axis eigenvalues (epr) and the epsilon
  132. %  for the Riccati solution positive definite tests (epp).  The two 
  133. %  methods are via eigenvalue or Schur  decomposition. HINF returns 
  134. %  the control law, k, the closed loop system, g, and the gamma 
  135. %  value achieved, gf.
  136. %
  137. %      [k,g,gf]=hinfsyn(p,nmeas,ncon,gmin,gmax,tol,ricmethd,epr,epp)
  138.  
  139.  pause   % strike any key to continue
  140.  
  141. % In this example, the system interconnection structure is 
  142. %  HIMAT_IC, with 2 measurements, 2 controls, a gamma low of 0.8, 
  143. %  gamma upper of 6, a tolerance on the gamma iteration of .1, 
  144. %  and we'll use the Schur decomposition method (2) to solve the 
  145. %  Riccati equations. The default values of epr (1e-10) and epp (1e-6)
  146. %  will be used for the epsilon tests.
  147.  
  148.  [k1,g1,gf1]=hinfsyn(himat_ic,2,2,0.8,6,0.05,2);
  149.  
  150. % A H-infinity control law has been designed which achieves an
  151. %  infinity norm of 1.6938 for the interconnection structure
  152. %  provided. First, we will examine aspects of the controller
  153. %  that was just designed, starting with the controller POLES.
  154.  
  155. rifd(spoles(k1))
  156. pause % strike any key to continue 
  157. clc
  158.  
  159. %  Next, a magnitude plot of the frequency response of k1
  160.  
  161. om2=logspace(0,4,16);
  162. k1_g = frsp(k1,om2);
  163. vplot('liv,lm',k1_g)
  164. echo off
  165. title('Magnitude plot of controller, k1')
  166. xlabel('strike any key to continue')
  167. pause
  168. echo on
  169. pause   % strike any key to continue
  170. format short e
  171. clc
  172.  
  173. %  Onto the closed loop, first checking that it is stable.
  174. %   Use the command SPOLES and RIFD.
  175.  
  176.  rifd(spoles(g1))
  177.  pause   % strike any key to continue
  178.  
  179. %  Now calculate a closed loop frequency response. This results in
  180. %  a 4x4 varying matrix G1G
  181.  
  182.  g1g=frsp(g1,om2);
  183.  
  184. %  Next, find the peak norm of the frequency response, using
  185. %   PKVNORM. This corresponds to the maximum singular of the
  186. %   closed loop system.
  187.  
  188.  g1pk=pkvnorm(g1g);
  189.  echo off
  190.  fprintf( '\n final gamma %g, \n max singular value %g \n \n',gf1,g1pk)
  191.  echo on
  192.  g1gs=vsvd(g1g);
  193.  vplot('liv,m',g1gs)
  194.  echo off
  195.  title('Singular value plot of the H_infinity closed loop');
  196.  xlabel('Frequency (rad/s)');
  197.  ylabel('Magnitude');
  198.  grid;
  199.  pause
  200.  echo on
  201.  pause   % strike any key to continue
  202.  
  203. % The H-infinity control law can be analyzed using mu-analysis.
  204. %  The closed-loop system, G1, has 4 inputs and 4 outputs. 
  205. %  The first two inputs and outputs correspond to the 
  206. %  uncertainty block, and the second two correspond to the 
  207. %  disturbance rejection block, or performance block.
  208. %  Therefore, we can define the uncertainty inputs and outputs
  209. %  as a full 2x2 uncertainty block and the disturbance 
  210. %  rejection inputs and outputs as a full 2x2 performance 
  211. %  block. If the mu value for this control design is 1, then
  212. %  we are able to achieve robust performance 
  213. %  for the set of defined weights and this control design.
  214.  
  215.  pause   % strike any key to continue
  216.  
  217. %   mu  -  usage: [bnds,row_d,sens,row_p] = mu(matin,blk,'opt')
  218.  
  219. %  The mu-analysis program, MU, requires a frequency response
  220. %  or constant matrix and a block structure. In this example,
  221. %  the frequency response is G1G and the block structure is
  222. %  two, 2x2 full blocks. MU returns the upper and lower 
  223. %  bounds for mu in BNDS1, the frequency
  224. %  varying d-scales DV1, the rational perturbation which 
  225. %  destabilizes the closed-loop plant RP1 and the 
  226. %  sensitivity of mu to the d-scales SENS1.
  227.  
  228.  blk=[2 2;2 2];
  229.  [bnds1,dv1,sens1,rp1]=mu(g1g,blk);
  230.  
  231. % Plot the maximum singular value and mu on the same plot
  232. %   solid line - maximum singular value
  233. %   dashed line - mu
  234.  
  235.  both=sbs(sel(g1gs,1,1),bnds1);
  236.  vplot('liv,m',both)
  237.  echo off
  238.  title('Max. singular value and mu for the  H-infinity Control Law');
  239.  xlabel('Frequency (rad/s)');
  240.  ylabel('Magnitude');
  241.  grid;
  242.  pause   % strike any key to continue
  243.  echo on
  244.  
  245. % clear variables not used and pack the memory space 
  246. %  due to the limitations on the Mac and PC's.
  247.  
  248.  clear g1g g1gs both 
  249.  pack
  250.  
  251. % This concludes the H-infinity control design portion of this
  252. %   demonstration program. If all you wanted was the H-infinity 
  253. %   control law, you can stop here just type in the number 1.
  254. %   Otherwise, we'll precede with the first iteration of 
  255. %   mu-synthesis.
  256.  
  257. echo off
  258.  fini = input(' Enter the number 1 to stop this program ');
  259.  if fini == 1
  260.    return
  261.  end
  262. echo on
  263.  
  264. % Designing an H-infinity control law is the first step in 
  265. %   mu-synthesis. The next step is to fit the d-scales which
  266. %   are output from the MU program. The MUSYNFIT function
  267. %   is designed to fit the d-scales for each individual block. 
  268. %   If there are n full blocks in the mu-problem, n-1 d-scales 
  269. %   are fit, since the d-scales are normalized to make 
  270. %   the last one 1.
  271.  
  272.  pause   % strike any key to continue
  273.  
  274.  musynfit 
  275.  
  276. %  MUSYNFIT requires the left d-scale from the previous 
  277. %   mu-synthesis iteration (PRE_DSYSL), the d-scales from the 
  278. %   mu-analysis problem (DDATA), and the sensitivity of 
  279. %   the d-scales (SENS), along with the block
  280. %   structure (BLK), number of controller inputs (NMEAS) 
  281. %   and outputs (NCNTRL). Since this is the first mu-synthesis 
  282. %   iteration, the left d-scale is input as 'first'.
  283. %   This creates an identity matrix of the correct size for dsysl.
  284. %   Output from MUSYNFIT is the left d-scale (DSYSL) and  
  285. %   right d-scale (DSYSR) to be multiplied into the 
  286. %   interconnection structure.
  287.  
  288.  pause    % strike any key to continue
  289.  
  290. %   The previous d-scale is required since if new d-scales
  291. %   are continuously generated of a given order, the number
  292. %   of states in the interconnection structure will grow 
  293. %   dramatically with each iteration. Therefore, the
  294. %   previous d-scale is normalized output of the new 
  295. %   d-scale data, and the new d-scale is wrapped around
  296. %   the original interconnection structure.
  297. %   For the first iteration, there are no d-scales in the problem.
  298.  
  299. %   We recommend a a 3rd order transfer function. This increases 
  300. %   the number of states in the interconnection structure to
  301. %   3*(size of full block)*2. If one selects a different order
  302. %   for the d-scale weight, the gamma iteration will converge
  303. %   to a different value.
  304.  
  305.  [dsysl,dsysr]=musynfit('first',dv1,sens1,blk,2,2);
  306.  
  307. % Now, wrap the new d-scales around the original 
  308. %  interconnection structure.
  309. %
  310. %          himat_ic2 = dsysl*himat_ic*(dsysr)^-1
  311.  
  312.  pause    % strike any key to continue
  313.  himat_ic2=mmult(dsysl,mmult(himat_ic,minv(dsysr)));
  314.  
  315.  minfo(himat_ic2)
  316. echo off
  317.  clear a b c dsysr fini g1 g1pk
  318.  pack;
  319. echo on
  320.  
  321. % Thus ends the first iteration of mu-synthesis. Notice that 
  322. %   the new interconnection structure has 20 states, 12 more 
  323. %   states than the original interconnection structure. These 
  324. %   12  states are due the frequency varying d-scales which 
  325. %   were fit. (Note, this is if you fit the d-scales using a 
  326. %   3rd order system.)
  327.  
  328. %  Type in 'himat_x2' to continue the mu-synthesis iteration
  329. %    procedure.
  330. %
  331. % Copyright MUSYN INC 1991,  All Rights Reserved
  332.