home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / scripts / control / dare.m < prev    next >
Text File  |  1996-09-28  |  3KB  |  113 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 = dare (a, b, c, r, opt)
  20.  
  21. # Usage: x = dare (a, b, c, r {,opt})
  22. #
  23. # Solves discrete-time algebraic riccati equation
  24. #
  25. #   a' x a - x + a' x b (r + b' x b)^{-1} b' x a + c = 0
  26. #
  27. # for
  28. #
  29. #   a: nxn
  30. #   b: nxm
  31. #   c: nxn, symmetric positive semidefinite 
  32. #   r: mxm, invertible
  33. #
  34. # If c is not square, then the function attempts to use c'*c instead.
  35. #
  36. # Solution method: Laub's Schur method (IEEE Trans Auto Contr, 1979) applied
  37. # to the appropriate symplectic matrix.
  38. #
  39. # See also: Ran and Rodman, "Stable Hermitian Solutions of Discrete
  40. # Algebraic Riccati Equations," Mathematics of Control, Signals and
  41. # Systems, Vol 5, no 2 (1992)  pp 165-194.
  42. #
  43. # opt is an option passed to the eigenvalue balancing routine default
  44. # is "B". 
  45. #
  46. # See also: balance, are
  47.  
  48. # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993.
  49.  
  50.   if (nargin == 4 || nargin == 5)
  51.     if (nargin == 5)
  52.       if (opt != "N" || opt != "P" || opt != "S" || opt != "B")
  53.     warning ("dare: opt has an invalid value -- setting to B");
  54.     opt = "B";
  55.       endif
  56.     else
  57.       opt = "B";
  58.     endif
  59.  
  60. # Check a matrix dimensions
  61.     if ((n = is_square (a)) == 0)
  62.       error ("dare: a is not square");
  63.     endif
  64.  
  65. # Check a,b compatibility.
  66.  
  67.     [n1, m] = size (b);
  68.  
  69.     if (n1 != n)
  70.       warning ("dare: a,b are not conformable");
  71.     endif
  72.  
  73.     if (is_controllable (a, b) == 0)
  74.       warning ("dare: a,b are not controllable");
  75.     endif
  76.  
  77. # Check a,c compatibility.
  78.  
  79.     if (is_observable (a, c) == 0)
  80.       warning ("dare: a,c are not observable");
  81.     endif
  82.  
  83.     if ((p = is_square (c)) == 0)
  84.       c = c'*c;
  85.       p = rows (c);
  86.     endif
  87.  
  88.     if (n != p)
  89.       error ("dare: a,c are not conformable");
  90.     endif
  91.  
  92. # Check r dimensions.
  93.  
  94.     if ((m1 = is_square (r)) == 0)
  95.       warning ("dare: r is not square");
  96.     elseif (m1 != m)
  97.       warning ("b,r are not conformable");
  98.     endif
  99.  
  100.     brb = (b/r)*b';
  101.     atc = a'\c;
  102.     [d, sy] = balance ([a + brb*atc, -brb/(a'); -atc, inv (a')], opt);
  103.     [u, s] = schur(sy,'D');
  104.     u = d*u;
  105.     n1 = n+1;
  106.     n2 = 2*n;
  107.     x = u (n1:n2, 1:n)/u(1:n, 1:n);
  108.   else
  109.     usage ("x = dare (a, b, c, r {,opt})");
  110.   endif
  111.  
  112. endfunction
  113.