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