home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / control / hinfdemo.m < prev    next >
Encoding:
Text File  |  1999-12-18  |  10.3 KB  |  330 lines
  1. ## Copyright (C) 1996,1998 Kai Mueller
  2. ##
  3. ## This file is part of Octave.
  4. ##
  5. ## Octave is free software; you can redistribute it and/or modify it
  6. ## under the terms of the GNU General Public License as published by the
  7. ## Free Software Foundation; either version 2, or (at your option) any
  8. ## later version.
  9. ##
  10. ## Octave is distributed in the hope that it will be useful, but WITHOUT
  11. ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12. ## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13. ## for more details.
  14. ##
  15. ## You should have received a copy of the GNU General Public License
  16. ## along with Octave; see the file COPYING.  If not, write to the Free
  17. ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
  18.  
  19. ## -*- texinfo -*-
  20. ## @deftypefn {Function File} {} hinfdemo ()
  21. ##
  22. ## H_infinity design demos for continuous SISO and MIMO systems and a
  23. ## discrete system.  The SISO system is difficult to control because it
  24. ## is non minimum phase and unstable.  The second design example
  25. ## controls the "jet707" plant, the linearized state space model of a
  26. ## Boeing 707-321 aircraft at v=80m/s (M = 0.26, Ga0 = -3 deg, alpha0 =
  27. ## 4 deg, kappa = 50 deg).  Inputs: (1) thrust and (2) elevator angle
  28. ## outputs: (1) airspeed and (2) pitch angle. The discrete system is a
  29. ## stable and second order.
  30. ##
  31. ## @table @asis
  32. ## @item SISO plant
  33. ## @display
  34. ## @group
  35. ##               s - 2
  36. ##    G(s) = --------------
  37. ##           (s + 2)(s - 1)
  38. ##
  39. ##                             +----+
  40. ##        -------------------->| W1 |---> v1
  41. ##    z   |                    +----+
  42. ##    ----|-------------+                   || T   ||     => min.
  43. ##        |             |                       vz   infty
  44. ##        |    +---+    v   y  +----+
  45. ##      u *--->| G |--->O--*-->| W2 |---> v2
  46. ##        |    +---+       |   +----+
  47. ##        |                |
  48. ##        |    +---+       |
  49. ##        -----| K |<-------
  50. ##             +---+
  51. ## @end group
  52. ## @end display
  53. ##    W1 und W2 are the robustness and performance weighting
  54. ##       functions
  55. ##
  56. ## @item MIMO plant
  57. ## The optimal controller minimizes the H_infinity norm of the
  58. ## augmented plant P (mixed-sensitivity problem):
  59. ## @display
  60. ## @group
  61. ##      w
  62. ##       1 -----------+
  63. ##                    |                   +----+
  64. ##                +---------------------->| W1 |----> z1
  65. ##      w         |   |                   +----+
  66. ##       2 ------------------------+
  67. ##                |   |            |
  68. ##                |   v   +----+   v      +----+
  69. ##             +--*-->o-->| G  |-->o--*-->| W2 |---> z2
  70. ##             |          +----+      |   +----+
  71. ##             |                      |
  72. ##             ^                      v
  73. ##              u (from                 y (to K)
  74. ##                controller
  75. ##                K)
  76. ##
  77. ##
  78. ##                   +    +           +    +
  79. ##                   | z  |           | w  |
  80. ##                   |  1 |           |  1 |
  81. ##                   | z  | = [ P ] * | w  |
  82. ##                   |  2 |           |  2 |
  83. ##                   | y  |           | u  |
  84. ##                   +    +           +    +
  85. ## @end group
  86. ## @end display
  87. ##
  88. ## @item DISCRETE SYSTEM
  89. ##   This is not a true discrete design. The design is carried out
  90. ##   in continuous time while the effect of sampling is described by
  91. ##   a bilinear transformation of the sampled system.
  92. ##   This method works quite well if the sampling period is "small"
  93. ##   compared to the plant time constants.
  94. ##
  95. ## @item The continuous plant
  96. ## @display
  97. ## @group
  98. ##                  1
  99. ##    G (s) = --------------
  100. ##     k      (s + 2)(s + 1)
  101. ##
  102. ## @end group
  103. ## @end display
  104. ## is discretised with a ZOH (Sampling period = Ts = 1 second):
  105. ## @display
  106. ## @group
  107. ## 
  108. ##              0.199788z + 0.073498
  109. ##    G(s) = --------------------------
  110. ##           (z - 0.36788)(z - 0.13534)
  111. ##
  112. ##                             +----+
  113. ##        -------------------->| W1 |---> v1
  114. ##    z   |                    +----+
  115. ##    ----|-------------+                   || T   ||     => min.
  116. ##        |             |                       vz   infty
  117. ##        |    +---+    v      +----+
  118. ##        *--->| G |--->O--*-->| W2 |---> v2
  119. ##        |    +---+       |   +----+
  120. ##        |                |
  121. ##        |    +---+       |
  122. ##        -----| K |<-------
  123. ##             +---+
  124. ## @end group
  125. ## @end display
  126. ##    W1 and W2 are the robustness and performancs weighting
  127. ##       functions
  128. ## @end table
  129. ## @end deftypefn
  130.  
  131. ## Kai P. Mueller 30-APR-1998 <mueller@ifr.ing.tu-bs.de
  132.  
  133. yn = [];
  134. while (length(yn) < 1)
  135.   yn = input(" * [s]iso, [m]imo, or [d]iscrete design? [no default]: ","S");
  136. endwhile
  137. if ((yn(1) == "s") | (yn(1) == 'S'))
  138.   sys_type = 1;
  139. elseif ((yn(1) == "m") | (yn(1) == 'M'))
  140.   sys_type = 2;
  141. elseif ((yn(1) == "d") | (yn(1) == 'D'))
  142.   sys_type = 3;
  143. else
  144.   disp(" *** no system type specified, hinfdemo terminated.");
  145.   return;
  146. endif
  147.  
  148. echo off
  149. switch (sys_type)
  150.  
  151.   case (1)
  152.     ## siso
  153.     disp(" ");
  154.     disp("    ----------------------------------------------");
  155.     disp("    H_infinity optimal control for the SISO plant:");
  156.     disp(" ");
  157.     disp("                    s - 2");
  158.     disp("        G(s) = --------------");
  159.     disp("               (s + 2)(s - 1)");
  160.     disp(" ");
  161.     disp("    ----------------------------------------------");
  162.     disp(" ");
  163.  
  164.     ## weighting on actuator u
  165.     W1 = wgt1o(0.05, 100.0, 425.0);
  166.     ## weighting on controlled variable y
  167.     W2 = wgt1o(10.0, 0.05, 0.001);
  168.     ## plant
  169.     G = tf2sys([1 -2],[1 1 -2]);
  170.  
  171.     ## need One as the pseudo transfer function One = 1
  172.     One = ugain(1);
  173.     disp(" o forming P...");
  174.     psys = buildssic([1 4;2 4;3 1],[3],[2 3 5],[3 4],G,W1,W2,One);
  175.     disp(" ");
  176.     disp(" o controller design...");
  177.     [K, gfin, GW]=hinfsyn(psys, 1, 1, 0.1, 10.0, 0.02);
  178.     disp(" ");
  179.     disp("-- OK ----------------------------------------------");
  180.  
  181.     disp("  Closed loop poles:");
  182.     damp(GW);
  183.     ## disp(" o Testing H_infinity norm: (hinfnorm does not work)");
  184.     ## hinfnorm(GW);
  185.  
  186.     disp(" ");
  187.     yn = input(" * Plot closed loop step response? [n]: ","S");
  188.     if (length(yn) >= 1)
  189.       if ((yn(1) == "y") || (yn(1) == 'Y'))
  190.           disp(" o step responses of T and KS...");
  191.           GW = buildssic([1 2; 2 1], [], [1 2], [-2], G, K);
  192.           figure(1);
  193.           step(GW, 1, 10);
  194.       endif
  195.     endif
  196.  
  197.   case (2)
  198.     ## mimo
  199.     disp(" ");
  200.     disp("    -----------------------------------------------");
  201.     disp("      H_inf optimal control for the jet707 plant");
  202.     disp("    -----------------------------------------------");
  203.     disp(" ");
  204.  
  205.     ## Weighting function on u (robustness weight)
  206.     ww1 = wgt1o(0.01,5,0.9);
  207.     ww2 = wgt1o(0.01,5,2.2);
  208.     W1 = buildssic([1 0;2 0],[],[1 2],[1 2],ww1,ww2);
  209.     ## Weighting function on y (performance weight)
  210.     ww1 = wgt1o(250,0.1,0.0001);
  211.     ww2 = wgt1o(250,0.1,0.0002);
  212.     W2 = buildssic([1 0;2 0],[],[1 2],[1 2],ww1,ww2);
  213.     ## plant (2 x 2 system)
  214.     G = jet707;
  215.  
  216.     disp(" o forming P...");
  217.     One = ugain(2);
  218.     Clst = [1 7; 2 8; 3 7; 4 8; 5 1; 6 2];
  219.     P = buildssic(Clst,[5 6],[3:6 9 10],[1 2 5:8],G,W1,W2,One);
  220.  
  221.     disp(" ");
  222.     disp(" o controller design...");
  223.     K = hinfsyn(P, 2, 2, 0.25, 10.0, 0.005);
  224.  
  225.     disp(" ");
  226.     yn = input(" * Plot closed loop step responses? [n]: ","S");
  227.     if (length(yn) >= 1)
  228.       if ((yn(1) == "y") || (yn(1) == 'Y'))
  229.           disp(" o step responses of T and KS...");
  230.           GW = buildssic([1 3;2 4;3 1;4 2],[],[1 2 3 4],[-3 -4],G,K);
  231.  
  232.           disp(" ");
  233.           disp("  FIGURE 1: speed refence => 1, pitch angle ref. => 0");
  234.           disp("  ===================================================");
  235.           disp("      y1:  speed                      (should be 1)");
  236.           disp("      y2:  pitch            angle (should remain 0)");
  237.           disp("      y3:  thrust      (should be a slow transient)");
  238.           disp("      y6:  elevator  (should be a faster transient)");
  239.           disp(" ");
  240.           disp("  FIGURE 2: speed refence => 0, pitch angle ref. => 1");
  241.           disp("  ===================================================");
  242.           disp("      y1:  speed                  (should remain 0)");
  243.           disp("      y2:  pitch                angle (should be 1)");
  244.           disp("      y3:  thrust      (should be a slow transient)");
  245.           disp("      y6:  elevator  (should be a faster transient)");
  246.           disp(" ");
  247.           figure(1)
  248.           step(GW);
  249.           figure(2)
  250.           step(GW,2);
  251.       endif
  252.     endif
  253.  
  254.   case (3)
  255.     ## discrete
  256.     disp(" ");
  257.     disp("    --------------------------------------------------");
  258.     disp("    Discrete H_infinity optimal control for the plant:");
  259.     disp(" ");
  260.     disp("                       0.199788z + 0.073498");
  261.     disp("            G(s) = --------------------------");
  262.     disp("                   (z - 0.36788)(z - 0.13533)");
  263.     disp("    --------------------------------------------------");
  264.     disp(" ");
  265.  
  266.     ## sampling time
  267.     Ts = 1.0;
  268.     ## weighting on actuator value u
  269.     W1 = wgt1o(0.1, 200.0, 50.0);
  270.     ## weighting on controlled variable y
  271.     W2 = wgt1o(350.0, 0.05, 0.0002);
  272.     ## omega axis
  273.     ww = logspace(-4.99, 3.99, 100);
  274.     if (columns(ww) > 1);  ww = ww';  endif
  275.  
  276.     ## continuous plant
  277.     G = tf2sys(2,[1 3 2]);
  278.     ## discrete plant with zoh
  279.     Gd = c2d(G, Ts);
  280.     ## w-plane (continuous representation of the sampled system)
  281.     Gw = d2c(Gd, "bi");
  282.  
  283.     disp(" ");
  284.     disp(" o building P...");
  285.     ## need One as the pseudo transfer function One = 1
  286.     One = ugain(1);
  287.     psys = buildssic([1 4;2 4;3 1],[3],[2 3 5],[3 4],Gw,W1,W2,One);
  288.     disp(" o controller design...");
  289.     [K, gfin, GWC] = hinfsyn(psys, 1, 1, 0.1, 10.0, 0.02);
  290.  
  291.     disp(" ");
  292.     fig_n = 1;
  293.     yn = input(" * Plot magnitudes of W1KS and W2S? [n]: ","S");
  294.     if (length(yn) >= 1)
  295.       if ((yn(1) == "y") || (yn(1) == 'Y'))
  296.         disp(" o magnitudes of W1KS and W2S...");
  297.         gwx = sysprune(GWC, 1, 1);
  298.         mag1 = bode(gwx, ww);
  299.         if (columns(mag1) > 1);  mag1 = mag1';  endif
  300.         gwx = sysprune(GWC, 2, 1);
  301.         mag2 = bode(gwx, ww);
  302.         if (columns(mag2) > 1);  mag2 = mag2';  endif
  303.         figure(fig_n)
  304.         fig_n = fig_n + 1;
  305.         gset grid
  306.         loglog(ww, [mag1 mag2]);
  307.       endif
  308.     endif
  309.  
  310.     Kd = c2d(K, "bi", Ts);
  311.     GG = buildssic([1 2; 2 1], [], [1 2], [-2], Gd, Kd);
  312.     disp(" o closed loop poles...");
  313.     damp(GG);
  314.  
  315.     disp(" ");
  316.     yn = input(" * Plot closed loop step responses? [n]: ","S");
  317.     if (length(yn) >= 1)
  318.       if ((yn(1) == "y") || (yn(1) == 'Y'))
  319.         disp(" o step responses of T and KS...");
  320.         figure(fig_n)
  321.         step(GG, 1, 10);
  322.       endif
  323.     endif
  324.  
  325. endswitch
  326.  
  327. disp(" o hinfdemo terminated successfully.");
  328.  
  329. ## KPM-hinfdemo/End
  330.