home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / control / is_dgkf.m < prev    next >
Text File  |  1999-03-05  |  7KB  |  229 lines
  1. # Copyright (C) 1996,1998 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 [retval,dgkf_struct] = is_dgkf(Asys,nu,ny,tol)
  20.   # function [retval,dgkf_struct] = is_dgkf(Asys,nu,ny{,tol})
  21.   # Determine whether a continuous time state space system meets
  22.   # assumptions of DGKF algorithm.  
  23.   # Partitions system into: [z] = [A  | Bw  Bu  ][w]
  24.   #                         [y]   [Cz | Dzw Dzu ][u]
  25.   #                               [Cy | Dyw Dyu ]
  26.   #If necessary, orthogonal transformations Qw, Qz and nonsingular
  27.   # transformations Ru, Ry are applied to respective vectors w, z, u, y 
  28.   # in order to satisfy DGKF assumptions.  Loop shifting is used if Dyu block
  29.   # is nonzero.
  30.   #
  31.   # inputs: 
  32.   #        Asys: structured system
  33.   #       nu: number of controlled inputs
  34.   #        ny: number of measured outputs
  35.   #        tol: threshhold for 0.  Default: 200*eps
  36.   # outputs: 
  37.   #    retval: true(1) if system passes check, false(0) otherwise
  38.   #    dgkf_struct: data structure of is_dgkf results.  Entries:
  39.   #      nw, nz: dimensions of w, z
  40.   #      A: system A matrix
  41.   #      Bw: (n x nw) Qw-transformed disturbance input matrix
  42.   #      Bu: (n x nu) Ru-transformed controlled input matrix;
  43.   #          Note: B = [Bw Bu] 
  44.   #      Cz: (nz x n) Qz-transformed error output matrix
  45.   #      Cy: (ny x n) Ry-transformed measured output matrix 
  46.   #          Note: C = [Cz; Cy] 
  47.   #      Dzu, Dyw: off-diagonal blocks of transformed D matrix that enter 
  48.   #          z, y from u, w respectively
  49.   #      Ru: controlled input transformation matrix 
  50.   #      Ry: observed output transformation matrix
  51.   #      Dyu_nz: nonzero if the Dyu block is nonzero.
  52.   #      Dyu: untransformed Dyu block
  53.   #      dflg: nonzero if the system is discrete-time
  54.   #   
  55.   #    is_dgkf exits with an error if the system is mixed discrete/continuous
  56.   #
  57.   # [1] Doyle, Glover, Khargonekar, Francis, "State Space Solutions
  58.   #     to Standard H2 and Hinf Control Problems," IEEE TAC August 1989
  59.   #
  60.   # [2] Maciejowksi, J.M.: "Multivariable feedback design,"
  61.   #     Addison-Wesley, 1989, ISBN 0-201-18243-2
  62.   #
  63.   # [3] Keith Glover and John C. Doyle: "State-space formulae for all
  64.   #     stabilizing controllers that satisfy and h-infinity-norm bound
  65.   #     and relations to risk sensitivity,"
  66.   #     Systems & Control Letters 11, Oct. 1988, pp 167-172.
  67.   # [4] P. A. Iglesias and K. Glover, "State-space approach to discrete-time
  68.   #     H-infinity control."  Int. J. Control, 1991, V. 54, #5, 1031-1073.
  69.   #
  70.   
  71.   #  Written by A. S. Hodel
  72.   #  Updated by John Ingram July 1996 to accept structured systems
  73.   #
  74.   # Revised by Kai P Mueller April 1998 to solve the general H_infinity
  75.   # problem using unitary transformations Q (on w and z)
  76.   # and non-singular transformations R (on u and y) such
  77.   # that the Dzu and Dyw matrices of the transformed plant
  78.   #
  79.   #    ~
  80.   #    P  (the variable Asys here)
  81.   #
  82.   # become
  83.   #
  84.   #    ~            -1         T
  85.   #    D  = Q   D   R   = [ 0 I ]  or [ I ],
  86.   #     12   12  12  12
  87.   #
  88.   #    ~            T
  89.   #    D  = R   D   Q   = [ 0 I ] or [ I ].
  90.   #     21   21  21  21
  91.   #
  92.   # This transformation together with the algorithm in [1] solves
  93.   # the general problem (see [2] for example). 
  94.  
  95.   if (nargin < 3) | (nargin > 4)
  96.     usage("[retval,dgkf_struct] = is_dgkf(Asys,nu,ny{,tol})");
  97.   elseif (! is_scalar(nu) | ! is_scalar(ny) )
  98.     error("is_dgkf: arguments 2 and 3 must be scalars")
  99.   elseif (! is_struct(Asys) )
  100.     error("Argument 1 must be a system data structure");
  101.   endif
  102.   if(nargin < 4)
  103.     tol = 200*eps;
  104.   elseif( !is_sample(tol) )
  105.     error("is_dgkf: tol must be a positive scalar")
  106.   endif
  107.  
  108.   retval = 1;        # assume passes test
  109.  
  110.   dflg = is_digital(Asys);
  111.   [Anc, Anz, nin, nout ] = sysdimensions(Asys);
  112.  
  113.   if( Anz == 0 & Anc == 0 )
  114.     error("is_dgkf: no system states");
  115.   elseif( nu >= nin )
  116.     error("is_dgkf: insufficient number of disturbance inputs");
  117.   elseif( ny >= nout )
  118.     error("is_dgkf: insufficient number of regulated outputs");
  119.   endif
  120.  
  121.   nw = nin - nu;           nw1 = nw + 1;
  122.   nz = nout - ny;          nz1 = nz + 1;
  123.  
  124.   [A,B,C,D] = sys2ss(Asys);
  125.   # scale input/output for numerical reasons
  126.   if(norm(C,'fro')*norm(B,'fro') == 0)
  127.     error("||C||*||B|| = 0; no dynamic connnection from inputs to outputs");
  128.   endif
  129.   xx = sqrt(norm(B, Inf) / norm(C, Inf));
  130.   B = B / xx;  C = C * xx;
  131.  
  132.   # partition matrices
  133.               Bw = B(:,1:nw);        Bu = B(:,nw1:nin);
  134.   Cz = C(1:nz,:);    Dzw = D(1:nz,1:nw);    Dzu = D(1:nz,nw1:nin);
  135.   Cy = C(nz1:nout,:);    Dyw = D(nz1:nout,1:nw);    Dyu = D(nz1:nout,nw1:nin);
  136.  
  137.   # Check for loopo shifting
  138.   Dyu_nz = (norm(Dyu,Inf) != 0);
  139.   if (Dyu_nz)
  140.     warning("is_dgkf: D22 nonzero; performing loop shifting");
  141.   endif
  142.  
  143.   # 12 - rank condition at w = 0
  144.   xx =[A, Bu; Cz, Dzu];
  145.   [nr, nc] = size(xx);
  146.   irank = rank(xx);
  147.   if (irank != nc)
  148.     retval = 0;
  149.     warning(sprintf("rank([A Bu; Cz Dzu]) = %d, need %d; n=%d, nz=%d, nu=%d", ...
  150.     irank,nc,(Anc+Anz),nz,nu));
  151.     warning(" *** 12-rank condition violated at w = 0.");
  152.   endif
  153.  
  154.   # 21 - rank condition at w = 0
  155.   xx =[A, Bw; Cy, Dyw];
  156.   [nr, nc] = size(xx);
  157.   irank = rank(xx);
  158.   if (irank != nr)
  159.     retval = 0;
  160.     warning(sprintf("rank([A Bw; Cy Dyw]) = %d, need %d; n=%d, ny=%d, nw=%d", ...
  161.     irank,nr,(Anc+Anz),ny,nw));
  162.     warning(" *** 21-rank condition violated at w = 0.");
  163.   endif
  164.  
  165.   # can Dzu be transformed to become [0 I]' or [I]?
  166.   # This ensures a normalized weight
  167.   [Qz, Ru] = qr(Dzu);
  168.   irank = rank(Ru);
  169.   if (irank != nu)
  170.     retval = 0;
  171.     warning(sprintf("*** rank(Dzu(%d x %d) = %d", nz, nu, irank));
  172.     warning(" *** Dzu does not have full column rank.");
  173.   endif
  174.   if (nu >= nz)
  175.     Qz = Qz(:,1:nu)';
  176.   else
  177.     Qz = [Qz(:,(nu+1):nz), Qz(:,1:nu)]';
  178.   endif
  179.   Ru = Ru(1:nu,:);
  180.  
  181.   # can Dyw be transformed to become [0 I] or [I]?
  182.   # This ensures a normalized weight
  183.   [Qw, Ry] = qr(Dyw');
  184.   irank = rank(Ry);
  185.   if (irank != ny)
  186.     retval = 0;
  187.     warning(sprintf("*** rank(Dyw(%d x %d) = %d", ny, nw, irank));
  188.     warning(" *** Dyw does not have full row rank.");
  189.   endif
  190.  
  191.   if (ny >= nw)
  192.     Qw = Qw(:,1:ny);
  193.   else
  194.     Qw = [Qw(:,(ny+1):nw), Qw(:,1:ny)];
  195.   endif
  196.   Ry = Ry(1:ny,:)';
  197.  
  198.   # transform P by Qz/Ru and Qw/Ry
  199.   Bw  = Bw*Qw;
  200.   Bu  = Bu/Ru;
  201.   B   = [Bw, Bu];
  202.   Cz  = Qz*Cz;
  203.   Cy  = Ry\Cy;
  204.   C   = [Cz; Cy];
  205.   Dzw = Qz*Dzw*Qw;
  206.   Dzu = Qz*Dzu/Ru;
  207.   Dyw = Ry\Dyw*Qw;
  208.  
  209.   # pack the return structure
  210.   dgkf_struct.nw    = nw;
  211.   dgkf_struct.nu    = nu;
  212.   dgkf_struct.nz    = nz;
  213.   dgkf_struct.ny    = ny;
  214.   dgkf_struct.A        = A;
  215.   dgkf_struct.Bw    = Bw;
  216.   dgkf_struct.Bu    = Bu;
  217.   dgkf_struct.Cz    = Cz;
  218.   dgkf_struct.Cy    = Cy;
  219.   dgkf_struct.Dzw    = Dzw;
  220.   dgkf_struct.Dzu    = Dzu;
  221.   dgkf_struct.Dyw    = Dyw;
  222.   dgkf_struct.Dyu    = Dyu;
  223.   dgkf_struct.Ru    = Ru;
  224.   dgkf_struct.Ry    = Ry;
  225.   dgkf_struct.Dyu_nz    = Dyu_nz;
  226.   dgkf_struct.dflg    = dflg;
  227.  
  228. endfunction
  229.