home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-bin / lib / octave / 1.1.1 / m / control / are.m < prev    next >
Encoding:
Text File  |  1996-10-12  |  2.6 KB  |  90 lines
  1. # Copyright (C) 1993, 1994, 1995 John W. Eaton
  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 x = are (a, b, c, opt)
  20.  
  21. # Usage: x = are (a, b, c {,opt})
  22. #
  23. # Solves algebraic riccati equation
  24. #
  25. #   a' x + x a - x b x + c = 0
  26. #
  27. # for identically dimensioned square matrices a, b, c.  If b (c) is not
  28. # square, then the function attempts to use b * b' (c' * c) instead.
  29. #
  30. # Solution method: apply Laub's Schur method (IEEE Trans. Auto. Contr,
  31. # 1979) to the appropriate Hamiltonian matrix.
  32. #
  33. # opt is an option passed to the eigenvalue balancing routine default is "B".
  34. #
  35. # See also: balance
  36.  
  37. # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993.
  38.  
  39.   if (nargin == 3 || nargin == 4)
  40.     if (nargin == 4)
  41.       if (! (strcmp (opt, "N") || strcmp (opt, "P") ...
  42.          || strcmp (opt, "S") || strcmp (opt, "B") ...
  43.          || strcmp (opt, "n") || strcmp (opt, "p") ...
  44.          || strcmp (opt, "s") || strcmp (opt, "b")))
  45.     warning ("are: opt has an invalid value; setting to B");
  46.     opt = "B";
  47.       endif
  48.     else
  49.       opt = "B";
  50.     endif
  51.     if ((n = is_square(a)) == 0)
  52.       error ("are: a is not square");
  53.     endif
  54.  
  55.     if (is_controllable(a,b) == 0)
  56.       warning ("are: a, b are not controllable");
  57.     endif
  58.     if ((m = is_square (b)) == 0)
  59.       b = b * b';
  60.       m = rows (b);
  61.     endif
  62.     if (is_observable (a, c) == 0)
  63.       warning ("are: a,c are not observable");
  64.     endif
  65.     if ((p = is_square (c)) == 0)
  66.       c = c' * c;
  67.       p = rows (c);
  68.     endif
  69.     if (n != m || n != p)
  70.       error ("are: a, b, c not conformably dimensioned.");
  71.     endif
  72.  
  73. # Should check for controllability/observability here
  74. # use Boley-Golub (Syst. Contr. Letters, 1984) method, not the
  75. #
  76. #                     n-1
  77. # rank ([ B A*B ... A^   *B]) method 
  78.  
  79.     [d, h] = balance ([a, -b; -c, -a'], opt);
  80.     [u, s] = schur (h, "A");
  81.     u = d * u;
  82.     n1 = n + 1;
  83.     n2 = 2 * n;
  84.     x = u (n1:n2, 1:n) / u (1:n, 1:n);
  85.   else
  86.     usage ("x = are (a, b, c)")
  87.   endif
  88.  
  89. endfunction
  90.