Mercurial > hg > octave-lyh

comparison scripts/control/base/dlyap.m @ 3431:99ab64f4a09d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2000-01-14 03:53:03 by jwe]
author jwe
date Fri, 14 Jan 2000 04:12:41 +0000
parents
children 7923abdeb4e5
comparison
equal deleted inserted replaced
3430:65b3519ac3a1 3431:99ab64f4a09d
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 ## @end deftypefn
61
62 ## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
63 ## Created: August 1993
64
65 function x = dlyap (a, b)
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