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