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