home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / control / dlyap.m < prev    next >
Encoding:
Text File  |  1999-12-15  |  3.2 KB  |  131 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} {@var{x} = } dlyap (@var{a}, @var{b})
  21. ## Solve the discrete-time Lyapunov equation
  22. ## 
  23. ##   @strong{Inputs}
  24. ##   @table @var
  25. ##     @item a
  26. ##     @var{n} by @var{n} matrix
  27. ##     @item b
  28. ##     Matrix: @var{n} by @var{n}, @var{n} by @var{m}, or @var{p} by @var{n}.
  29. ##   @end table
  30. ## 
  31. ##   @strong{Outputs}
  32. ##   @var{x}: matrix satisfying appropriate discrete time Lyapunov equation.
  33. ##   Options:
  34. ##   @itemize @bullet
  35. ##     @item @var{b} is square: solve @code{a x a' - x + b = 0}
  36. ##     @item @var{b} is not square: @var{x} satisfies either
  37. ##       @example
  38. ##       a x a' - x + b b' = 0
  39. ##       @end example
  40. ##       @noindent
  41. ##       or
  42. ##       @example
  43. ##  a' x a - x + b' b = 0,
  44. ##     @end example
  45. ##     @noindent
  46. ##     whichever is appropriate.
  47. ##   @end itemize
  48. ##   
  49. ## @strong{Method}
  50. ##   Uses Schur decomposition method as in Kitagawa,
  51. ##     @cite{An Algorithm for Solving the Matrix Equation @var{X} =
  52. ##     @var{F}@var{X}@var{F}' + @var{S}},
  53. ##   International Journal of Control, Volume 25, Number 5, pages 745--753
  54. ##   (1977). 
  55. ## 
  56. ## Column-by-column solution method as suggested in
  57. ##   Hammarling, @cite{Numerical Solution of the Stable, Non-Negative
  58. ##   Definite Lyapunov Equation}, IMA Journal of Numerical Analysis, Volume
  59. ##   2, pages 303--323 (1982).
  60. ## 
  61. ## @end deftypefn
  62.  
  63. function x = dlyap (a, b)
  64.  
  65.   ## Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993.
  66.  
  67.   if ((n = is_square (a)) == 0)
  68.     warning ("dlyap: a must be square");
  69.   endif
  70.  
  71.   if ((m = is_square (b)) == 0)
  72.     [n1, m] = size (b);
  73.     if (n1 == n)
  74.       b = b*b';
  75.       m = n1;
  76.     else
  77.       b = b'*b;
  78.       a = a';
  79.     endif
  80.   endif
  81.  
  82.   if (n != m)
  83.     warning ("dlyap: a,b not conformably dimensioned");
  84.   endif
  85.  
  86.   ## Solve the equation column by column.
  87.  
  88.   [u, s] = schur (a);
  89.   b = u'*b*u;
  90.  
  91.   j = n;
  92.   while (j > 0)
  93.     j1 = j;
  94.  
  95.     ## Check for Schur block.
  96.  
  97.     if (j == 1)
  98.       blksiz = 1;
  99.     elseif (s (j, j-1) != 0)
  100.       blksiz = 2;
  101.       j = j - 1;
  102.     else
  103.       blksiz = 1;
  104.     endif
  105.  
  106.     Ajj = kron (s (j:j1, j:j1), s) - eye (blksiz*n);
  107.  
  108.     rhs = reshape (b (:, j:j1), blksiz*n, 1);
  109.  
  110.     if (j1 < n)
  111.       rhs2 = s*(x (:, (j1+1):n) * s (j:j1, (j1+1):n)');
  112.       rhs = rhs + reshape (rhs2, blksiz*n, 1);
  113.     endif
  114.  
  115.     v = - Ajj\rhs;
  116.     x (:, j) = v (1:n);
  117.  
  118.     if(blksiz == 2)
  119.       x (:, j1) = v ((n+1):blksiz*n);
  120.     endif
  121.  
  122.     j = j - 1;
  123.  
  124.   endwhile
  125.  
  126.   ## Back-transform to original coordinates.
  127.  
  128.   x = u*x*u';
  129.  
  130. endfunction
  131.