home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / dgkfdemo.m < prev    next >
Text File  |  1999-12-24  |  12KB  |  353 lines
  1. ## Copyright (C) 1996 Auburn University.  All Rights Reserved
  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 } { } dgkfdemo ( ) 
  21. ## Octave Controls toolbox demo: H2/Hinfinity options demos
  22. ##@end deftypefn
  23.  
  24. function dgkfdemo ()
  25. ## Written by A. S. Hodel June 1995
  26.  
  27.   save_val = page_screen_output;
  28.   page_screen_output = 1;
  29.   while (1)
  30.     clc
  31.     menuopt=0;
  32.     while(menuopt > 10 || menuopt < 1)
  33.       menuopt = menu('Octave H2/Hinfinity options demo', ...
  34.         'LQ regulator', ...
  35.         'LG state estimator', ...
  36.         'LQG optimal control design', ...
  37.         'H2 gain of a system', ...
  38.         'H2 optimal controller of a system', ...
  39.         'Hinf gain of a system', ...
  40.         'Hinf optimal controller of a SISO system', ...
  41.         'Hinf optimal controller of a MIMO system', ...
  42.         'Discrete-time Hinf optimal control by bilinear transform', ...
  43.         'Return to main demo menu');
  44.     endwhile
  45.     if (menuopt == 1)
  46.       disp('Linear/Quadratic regulator design:')
  47.       disp('Compute optimal state feedback via the lqr command...')
  48.       help lqr
  49.       disp(' ')
  50.       disp('Example:')
  51.       A = [0, 1; -2, -1]
  52.       B = [0; 1]
  53.       Q = [1, 0; 0, 0]
  54.       R = 1
  55.       disp("Q = state penalty matrix; R = input penalty matrix")
  56.       prompt
  57.       disp('Compute state feedback gain k, ARE solution P, and closed-loop')
  58.       disp('poles as follows:');
  59.       cmd = "[k, p, e] = lqr(A,B,Q,R)";
  60.       run_cmd
  61.       prompt
  62.       disp("A similar approach can be used for LTI discrete-time systems")
  63.       disp("by using the dlqr command in place of lqr (see LQG example).")
  64.     elseif (menuopt == 2)
  65.       disp('Linear/Gaussian estimator design:')
  66.       disp('Compute optimal state estimator via the lqe command...')
  67.       help lqe
  68.       disp(' ')
  69.       disp('Example:')
  70.       A = [0, 1; -2, -1]
  71.       disp("disturbance entry matrix G")
  72.       G = eye(2)
  73.       disp("Output measurement matrix C")
  74.       C = [0, 1]
  75.       SigW = [1, 0; 0, 1]
  76.       SigV = 1
  77.       disp("SigW = input disturbance intensity matrix;")
  78.       disp("SigV = measurement noise intensity matrix")
  79.       prompt
  80.       disp('Compute estimator feedback gain k, ARE solution P, and estimator')
  81.       disp('poles via the command: ')
  82.       cmd = "[k, p, e] = lqe(A,G,C,SigW,SigV)";
  83.       run_cmd
  84.       disp("A similar approach can be used for LTI discrete-time systems")
  85.       disp("by using the dlqe command in place of lqe (see LQG example).")
  86.     elseif (menuopt == 3)
  87.       disp('LQG optimal controller of a system:')
  88.       disp('Input accepted as either A,B,C matrices or in system data structure form')
  89.       disp('in both discrete and continuous time.')
  90.       disp("Example 1: continuous time design:")
  91.       prompt
  92.       help lqg
  93.       disp("Example system")
  94.       A = [0, 1; .5, .5];
  95.       B = [0; 2];
  96.       G = eye(2)
  97.       C = [1, 1];
  98.       sys = ss2sys(A, [B, G], C);
  99.       sys = syssetsg(sys,"in", ...
  100.                ["control input"; "disturbance 1"; "disturbance 2"]);
  101.       sysout(sys)
  102.       prompt
  103.       disp("Filtering/estimator parameters:")
  104.       SigW = eye(2)
  105.       SigV = 1
  106.       prompt
  107.       disp("State space (LQR) parameters Q and R are:")
  108.       Q = eye(2)
  109.       R = 1
  110.       cmd = "[K,Q1,P1,Ee,Er] = lqg(sys,SigW,SigV,Q,R,1);";
  111.       run_cmd
  112.       disp("Check: closed loop system A-matrix is")
  113.       disp(" [A,      B*Cc]")
  114.       disp(" [Bc*C,   Ac  ]")
  115.       cmd = "[Ac, Bc, Cc] = sys2ss(K);";
  116.       run_cmd
  117.       cmd = "Acl = [A, B*Cc; Bc*C, Ac]";
  118.       run_cmd
  119.       disp("Check: poles of Acl:")
  120.       Acl_poles = sortcom(eig(Acl))
  121.       disp("Predicted poles from design = union(Er,Ee)")
  122.       cmd = "pred_poles = sortcom([Er; Ee])";
  123.       run_cmd
  124.       disp("Example 2: discrete-time example")
  125.       cmd1 = "Dsys = ss2sys(A, [G, B], C, [0, 0, 0], 1);";
  126.       cmd2 = "[K,Q1,P1,Ee,Er] = lqg(Dsys,SigW, SigV,Q,R);";
  127.       disp("Run commands:")
  128.       cmd = cmd1;
  129.       run_cmd
  130.       cmd = cmd2;
  131.       run_cmd
  132.       prompt
  133.       disp("Check: closed loop system A-matrix is")
  134.       disp(" [A,      B*Cc]")
  135.       disp(" [Bc*C,   Ac  ]")
  136.       [Ac,Bc,Cc] = sys2ss(K);
  137.       Acl = [A, B*Cc; Bc*C, Ac]
  138.       prompt
  139.       disp("Check: poles of Acl:")
  140.       Acl_poles = sortcom(eig(Acl))
  141.       disp("Predicted poles from design = union(Er,Ee)")
  142.       pred_poles = sortcom([Er;Ee])
  143.     elseif (menuopt == 4)
  144.       disp('H2 gain of a system: (Energy in impulse response)')
  145.       disp('Example 1: Stable plant:')
  146.       cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
  147.       run_cmd
  148.       disp("Put into Packed system form:")
  149.       cmd = "Asys = ss2sys(A,B,C);";
  150.       run_cmd
  151.       disp("Evaluate system 2-norm (impulse response energy):");
  152.       cmd = "AsysH2 = h2norm(Asys)";
  153.       run_cmd
  154.       disp("Compare with a plot of the system impulse response:")
  155.       tt = 0:0.1:20;
  156.       for ii=1:length(tt)
  157.         ht(ii) = C*expm(A*tt(ii))*B;
  158.       endfor
  159.       plot(tt,ht)
  160.       title("impulse response of example plant")
  161.       prompt
  162.       disp('Example 2: unstable plant')
  163.       cmd = "A = [0, 1; 2, 1]";
  164.       eval(cmd);
  165.       cmd = "B = [0; 1]";
  166.       eval(cmd);
  167.       cmd = "C = [1, 0]";
  168.       eval(cmd);
  169.       cmd = "sys_poles = eig(A)";
  170.       run_cmd
  171.       prompt
  172.       disp('Put into system data structure form:')
  173.       cmd="Bsys = ss2sys(A,B,C);";
  174.       run_cmd
  175.       disp('Evaluate 2-norm:')
  176.       cmd = "BsysH2 = h2norm(Bsys)";
  177.       run_cmd
  178.       disp(' ')
  179.       prompt('NOTICE: program returns a value without an error signal.')
  180.       disp('')
  181.  
  182.     elseif (menuopt == 5)
  183.       disp('H2 optimal controller of a system: command = h2syn:')
  184.       prompt
  185.       help h2syn
  186.       prompt
  187.       disp("Example system: double integrator with output noise and")
  188.       disp("input disturbance:")
  189.       disp(" ");
  190.       disp("       -------------------->y2");
  191.       disp("       |   _________");
  192.       disp("u(t)-->o-->| 1/s^2 |-->o-> y1");
  193.       disp("       ^   ---------   ^");
  194.       disp("       |               |");
  195.       disp("      w1(t)           w2(t)");
  196.       disp(" ")
  197.       disp("w enters the system through B1, u through B2")
  198.       disp("z = [y1; y2] is obtained through C1, y=y1 through C2");
  199.       disp(" ")
  200.       cmd = "A = [0, 1; 0, 0];  B1 = [0, 0; 1, 0]; B2 = [0; 1];";
  201.       disp(cmd)
  202.       eval(cmd);
  203.       cmd = "C1 = [1, 0; 0, 0]; C2 = [1, 0];    D11 = zeros(2);";
  204.       disp(cmd)
  205.       eval(cmd);
  206.       cmd = "D12 = [0; 1];  D21 = [0, 1];  D22 = 0; D = [D11, D12; D21, D22];";
  207.       disp(cmd)
  208.       eval(cmd);
  209.       disp("Design objective: compute U(s)=K(s)Y1(s) to minimize the closed")
  210.       disp("loop impulse response from w(t) =[w1; w2] to z(t) = [y1; y2]");
  211.       prompt
  212.       disp("First: pack system:")
  213.       cmd="Asys = ss2sys(A, [B1, B2], [C1; C2], D);";
  214.       run_cmd
  215.       disp("Open loop multivariable Bode plot: (will take a moment)")
  216.       cmd="bode(Asys);";
  217.       run_cmd
  218.       prompt("Press a key to close plot and continue");
  219.       closeplot
  220.       disp("Controller design command: (only need 1st two output arguments)")
  221.       cmd="[K,gain, Kc, Kf, Pc,  Pf] = h2syn(Asys,1,1);";
  222.       run_cmd
  223.       disp("Controller is:")
  224.       cmd = "sysout(K)";
  225.       run_cmd
  226.       disp(["returned gain value is: ",num2str(gain)]);
  227.       disp("Check: close the loop and then compute h2norm:")
  228.       prompt
  229.       cmd="K_loop = sysgroup(Asys,K);";
  230.       run_cmd
  231.       cmd = "Kcl = syscnct(K_loop,[3,4],[4,3]);";
  232.       run_cmd
  233.       cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);";
  234.       run_cmd
  235.       cmd="gain_Kcl = h2norm(Kcl)";
  236.       run_cmd
  237.       cmd="gain_err = gain_Kcl - gain";
  238.       run_cmd
  239.       disp("Check: multivarible bode plot:")
  240.       cmd="bode(Kcl);";
  241.       run_cmd
  242.       prompt
  243.       disp("Related functions: is_dgkf, is_contr, is_stabi,")
  244.       disp("                is_obsrv, is_detec")
  245.     elseif (menuopt == 6)
  246.       disp('Hinfinity gain of a system: (max gain over all j-omega)')
  247.       disp('Example 1: Stable plant:')
  248.       cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
  249.       run_cmd
  250.       disp('Pack into system format:')
  251.       cmd = "Asys = ss2sys(A,B,C);";
  252.       run_cmd
  253.       disp('The infinity norm must be computed iteratively by')
  254.       disp('binary search.  For this example, we select tolerance tol = 0.01, ')
  255.       disp('min gain gmin = 1e-2, max gain gmax=1e4.')
  256.       disp('Search quits when upper bound <= (1+tol)*lower bound.')
  257.       cmd = "tol = 0.01; gmin = 1e-2; gmax = 1e+4;";
  258.       run_cmd
  259.       cmd = "[AsysHinf,gmin,gmax] = hinfnorm(Asys,tol,gmin,gmax)"
  260.       run_cmd
  261.       disp("Check: look at max value of magntude Bode plot of Asys:");
  262.       [M,P,w] = bode(Asys);
  263.       xlabel('Omega')
  264.       ylabel('|Asys(j omega)| ')
  265.       grid();
  266.       semilogx(w,M);
  267.       disp(["Max magnitude is ",num2str(max(M)), ...
  268.     ", compared with gmin=",num2str(gmin)," and gmax=", ...
  269.         num2str(gmax),"."])
  270.       prompt
  271.       disp('Example 2: unstable plant')
  272.       cmd = "A = [0, 1; 2, 1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
  273.       run_cmd
  274.       disp("Pack into system format:")
  275.       cmd = "Bsys = ss2sys(A,B,C);";
  276.       run_cmd
  277.       disp('Evaluate with BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)')
  278.       BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)
  279.       disp(' ')
  280.       disp('NOTICE: program returns a value without an error signal.')
  281.       disp('')
  282.  
  283.     elseif (menuopt == 7)
  284.       disp('Hinfinity optimal controller of a system: command = hinfsyn:')
  285.       prompt
  286.       help hinfsyn
  287.       prompt
  288.       disp("Example system: double integrator with output noise and")
  289.       disp("input disturbance:")
  290.       A = [0, 1; 0, 0]
  291.       B1 = [0, 0; 1, 0]
  292.       B2 = [0; 1]
  293.       C1 = [1, 0; 0, 0]
  294.       C2 = [1, 0]
  295.       D11 = zeros(2);
  296.       D12 = [0; 1];
  297.       D21 = [0, 1];
  298.       D22 = 0;
  299.       D = [D11, D12; D21, D22]
  300.       prompt
  301.       disp("First: pack system:")
  302.       cmd="Asys = ss2sys(A, [B1, B2], [C1; C2], D);";
  303.       run_cmd
  304.       prompt
  305.       disp("Open loop multivariable Bode plot: (will take a moment)")
  306.       cmd="bode(Asys);";
  307.       run_cmd
  308.       prompt
  309.       disp("Controller design command: (only need 1st two output arguments)")
  310.       gmax = 1000
  311.       gmin = 0.1
  312.       gtol = 0.01
  313.       cmd="[K,gain] = hinfsyn(Asys,1,1,gmin,gmax,gtol);";
  314.       run_cmd
  315.       disp("Check: close the loop and then compute h2norm:")
  316.       prompt
  317.       cmd="K_loop = sysgroup(Asys,K);";
  318.       run_cmd
  319.       cmd = "Kcl = syscnct(K_loop,[3,4],[4,3]);";
  320.       run_cmd
  321.       cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);";
  322.       run_cmd
  323.       cmd="gain_Kcl = hinfnorm(Kcl)";
  324.       run_cmd
  325.       cmd="gain_err = gain_Kcl - gain";
  326.       run_cmd
  327.       disp("Check: multivarible bode plot:")
  328.       cmd="bode(Kcl);";
  329.       run_cmd
  330.       prompt
  331.       disp("Related functions: is_dgkf, is_contr, is_stabi,")
  332.       disp("                   is_obsrv, is_detec, buildssc")
  333.     elseif (menuopt == 8)
  334.       disp('Hinfinity optimal controller of MIMO system: command = hinfsyn:')
  335.       prompt
  336.       help hinfsyn
  337.       prompt
  338.       disp("Example system: Boeing 707-321 airspeed/pitch angle control")
  339.       disp(" ")
  340.       hinfdemo
  341.     elseif (menuopt == 9)
  342.       disp("Discrete time H-infinity control via bilinear transform");
  343.       prompt
  344.       dhinfdm
  345.     elseif (menuopt == 10)
  346.       return
  347.     endif
  348.     prompt
  349.   endwhile  
  350.   page_screen_output = save_val;
  351. endfunction
  352.  
  353.