home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / are.m < prev    next >
Text File  |  1999-12-24  |  3KB  |  118 lines
  1. ## Copyright (C) 1993, 1994, 1995 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} {} are (@var{a}, @var{b}, @var{c}, @var{opt})
  21. ## Solve the algebraic Riccati equation
  22. ## @iftex
  23. ## @tex
  24. ## $$
  25. ## A^TX + XA - XBX + C = 0
  26. ## $$
  27. ## @end tex
  28. ## @end iftex
  29. ## @ifinfo
  30. ## @example
  31. ## a' * x + x * a - x * b * x + c = 0
  32. ## @end example
  33. ## @end ifinfo
  34. ## 
  35. ## @strong{Inputs}
  36. ## @noindent
  37. ## for identically dimensioned square matrices 
  38. ## @table @var
  39. ## @item a
  40. ## @var{n}x@var{n} matrix.
  41. ## @item b
  42. ##   @var{n}x@var{n} matrix or @var{n}x@var{m} matrix; in the latter case
  43. ##   @var{b} is replaced by @math{b:=b*b'}.
  44. ## @item c
  45. ##   @var{n}x@var{n} matrix or @var{p}x@var{m} matrix; in the latter case
  46. ##   @var{c} is replaced by @math{c:=c'*c}.
  47. ## @item opt
  48. ## (optional argument; default = @code{"B"}):
  49. ## String option passed to @code{balance} prior to ordered Schur decomposition.
  50. ## @end table
  51. ## 
  52. ## @strong{Outputs}
  53. ## @var{x}: solution of the ARE.
  54. ## 
  55. ## @strong{Method}
  56. ## Laub's Schur method (IEEE Transactions on
  57. ## Automatic Control, 1979) is applied to the appropriate Hamiltonian
  58. ## matrix.
  59. ## 
  60. ## @end deftypefn
  61.  
  62. ## See also: balance, dare
  63.  
  64. function x = are (a, b, c, opt)
  65. ## Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993.
  66.  
  67.   if (nargin == 3 || nargin == 4)
  68.     if (nargin == 4)
  69.       if (! (strcmp (opt, "N") || strcmp (opt, "P") ...
  70.          || strcmp (opt, "S") || strcmp (opt, "B") ...
  71.          || strcmp (opt, "n") || strcmp (opt, "p") ...
  72.          || strcmp (opt, "s") || strcmp (opt, "b")))
  73.     warning ("are: opt has an invalid value; setting to B");
  74.     opt = "B";
  75.       endif
  76.     else
  77.       opt = "B";
  78.     endif
  79.     if ((n = is_sqr(a)) == 0)
  80.       error ("are: a is not square");
  81.     endif
  82.  
  83.     if (is_contr(a,b) == 0)
  84.       warning ("are: a, b are not controllable");
  85.     endif
  86.     if ((m = is_sqr (b)) == 0)
  87.       b = b * b';
  88.       m = rows (b);
  89.     endif
  90.     if (is_obsrv (a, c) == 0)
  91.       warning ("are: a,c are not observable");
  92.     endif
  93.     if ((p = is_sqr (c)) == 0)
  94.       c = c' * c;
  95.       p = rows (c);
  96.     endif
  97.     if (n != m || n != p)
  98.       error ("are: a, b, c not conformably dimensioned.");
  99.     endif
  100.  
  101. ## Should check for controllability/observability here
  102. ## use Boley-Golub (Syst. Contr. Letters, 1984) method, not the
  103. ##
  104. ##                     n-1
  105. ## rank ([ B A*B ... A^   *B]) method 
  106.  
  107.     [d, h] = balance ([a, -b; -c, -a'], opt);
  108.     [u, s] = schur (h, "A");
  109.     u = d * u;
  110.     n1 = n + 1;
  111.     n2 = 2 * n;
  112.     x = u (n1:n2, 1:n) / u (1:n, 1:n);
  113.   else
  114.     usage ("x = are (a, b, c)")
  115.   endif
  116.  
  117. endfunction
  118.