home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / dare.m < prev    next >
Text File  |  1999-12-24  |  3KB  |  124 lines
  1. ## Copyright (C) 1996, 1997 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
  7. ## the Free Software Foundation; either version 2, or (at your option)
  8. ## any later version.
  9. ##
  10. ## Octave is distributed in the hope that it will be useful, but
  11. ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13. ## General Public License 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
  18. ## 02111-1307, USA.
  19.  
  20. ## -*- texinfo -*-
  21. ## @deftypefn {Function File} {} dare (@var{a}, @var{b}, @var{c}, @var{r}, @var{opt})
  22. ## 
  23. ## Return the solution, @var{x} of the discrete-time algebraic Riccati
  24. ## equation
  25. ## @iftex
  26. ## @tex
  27. ## $$
  28. ## A^TXA - X + A^TXB (R + B^TXB)^{-1} B^TXA + C = 0
  29. ## $$
  30. ## @end tex
  31. ## @end iftex
  32. ## @ifinfo
  33. ## @example
  34. ## a' x a - x + a' x b (r + b' x b)^(-1) b' x a + c = 0
  35. ## @end example
  36. ## @end ifinfo
  37. ## @noindent
  38. ## 
  39. ## @strong{Inputs}
  40. ## @table @var
  41. ## @item a
  42. ## @var{n} by @var{n}.
  43. ## 
  44. ## @item b
  45. ## @var{n} by @var{m}.
  46. ## 
  47. ## @item c
  48. ## @var{n} by @var{n}, symmetric positive semidefinite, or @var{p} by @var{n}.
  49. ## In the latter case @math{c:=c'*c} is used.
  50. ## 
  51. ## @item r
  52. ## @var{m} by @var{m}, symmetric positive definite (invertible).
  53. ## 
  54. ## @item opt
  55. ## (optional argument; default = @code{"B"}):
  56. ## String option passed to @code{balance} prior to ordered @var{QZ} decomposition.
  57. ## @end table
  58. ## 
  59. ## @strong{Outputs}
  60. ## @var{x} solution of DARE.
  61. ## 
  62. ## @strong{Method}
  63. ## Generalized eigenvalue approach (Van Dooren; SIAM J.
  64. ##  Sci. Stat. Comput., Vol 2) applied  to the appropriate symplectic pencil.
  65. ## 
  66. ##  See also: Ran and Rodman, "Stable Hermitian Solutions of Discrete
  67. ##  Algebraic Riccati Equations," Mathematics of Control, Signals and
  68. ##  Systems, Vol 5, no 2 (1992)  pp 165-194.
  69. ## 
  70. ## @end deftypefn
  71.  
  72. ## See also: balance, are
  73.  
  74. ## Author: A. S. Hodel <scotte@eng.auburn.edu>
  75. ## Created: August 1993
  76. ## Adapted-By: jwe
  77.  
  78. function x = dare (a, b, c, r, opt)
  79.  
  80.   if (nargin == 4 | nargin == 5)
  81.     if (nargin == 5)
  82.       if (opt != "N" || opt != "P" || opt != "S" || opt != "B")
  83.     warning ("dare: opt has an invalid value -- setting to B");
  84.     opt = "B";
  85.       endif
  86.     else
  87.       opt = "B";
  88.     endif
  89.  
  90.     ## dimension checks are done in is_contr, is_obsrv
  91.     if (is_contr (a, b) == 0)
  92.       warning ("dare: a,b are not controllable");
  93.     elseif (is_obsrv (a, c) == 0)
  94.       warning ("dare: a,c are not observable");
  95.     endif
  96.  
  97.     if ((p = is_sqr (c)) == 0)
  98.       c = c'*c;
  99.       p = rows (c);
  100.     endif
  101.  
  102.     ## Check r dimensions.
  103.     n = rows(a);
  104.     m = columns(b);
  105.     if ((m1 = is_sqr (r)) == 0)
  106.       warning ("dare: r is not square");
  107.     elseif (m1 != m)
  108.       warning ("b,r are not conformable");
  109.     endif
  110.  
  111.     s1 = [a, zeros(n) ; -c, eye(n)];
  112.     s2 = [eye(n), (b/r)*b' ; zeros(n), a'];
  113.     [c,d,s1,s2] = balance(s1,s2,opt);
  114.     [aa,bb,u,lam] = qz(s1,s2,"S");
  115.     u = d*u;
  116.     n1 = n+1;
  117.     n2 = 2*n;
  118.     x = u (n1:n2, 1:n)/u(1:n, 1:n);
  119.   else
  120.     usage ("x = dare (a, b, c, r {,opt})");
  121.   endif
  122.  
  123. endfunction
  124.