|
55
|
1 function x = are (a, b, c, opt) |
|
26
|
2 |
|
73
|
3 # Usage: x = are (a, b, c {,opt}) |
|
26
|
4 # |
|
55
|
5 # Solves algebraic riccati equation |
|
26
|
6 # |
|
55
|
7 # a' x + x a - x b x + c = 0 |
|
26
|
8 # |
|
55
|
9 # for identically dimensioned square matrices a, b, c. If b (c) is not |
|
|
10 # square, then the function attempts to use b * b' (c' * c) instead. |
|
26
|
11 # |
|
55
|
12 # opt is an option passed to the eigenvalue balancing routine; default |
|
|
13 # is `B'. |
|
26
|
14 # |
|
73
|
15 # See also: balance |
|
|
16 |
|
|
17 # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993. |
|
26
|
18 |
|
55
|
19 if (nargin == 3 || nargin == 4) |
|
|
20 if (nargin == 4) |
|
|
21 if (! (strcmp (opt, "N") || strcmp (opt, "P") ... |
|
|
22 || strcmp (opt, "S") || strcmp (opt, "B") ... |
|
|
23 || strcmp (opt, "n") || strcmp (opt, "p") ... |
|
|
24 || strcmp (opt, "s") || strcmp (opt, "b"))) |
|
|
25 printf ("warning: are: opt has an illegal value; setting to B"); |
|
|
26 opt = "B"; |
|
|
27 endif |
|
|
28 else |
|
|
29 opt = "B"; |
|
|
30 endif |
|
|
31 if ((n = is_square(a)) == 0) |
|
|
32 error ("are: a is not square"); |
|
|
33 endif |
|
|
34 |
|
|
35 if (is_controllable(a,b) == 0) |
|
|
36 printf("warning: are: a, b are not controllable"); |
|
26
|
37 endif |
|
55
|
38 if ((m = is_square (b)) == 0) |
|
|
39 b = b * b'; |
|
|
40 m = rows (b); |
|
|
41 endif |
|
|
42 if (is_observable (a, c) == 0) |
|
|
43 printf ("warning: are: a,c are not observable"); |
|
|
44 endif |
|
|
45 if ((p = is_square (c)) == 0) |
|
|
46 c = c' * c; |
|
|
47 p = rows (c); |
|
|
48 endif |
|
|
49 if (n != m || n != p) |
|
|
50 error ("are: a, b, c not conformably dimensioned."); |
|
|
51 endif |
|
|
52 |
|
|
53 # Should check for controllability/observability here |
|
|
54 # use Boley-Golub (Syst. Contr. Letters, 1984) method, not the |
|
|
55 # |
|
|
56 # n-1 |
|
|
57 # rank ([ B A*B ... A^ *B]) method |
|
|
58 |
|
|
59 [u, s] = schur (balance ([a, -b; -c, -a'], opt), "A"); |
|
|
60 n1 = n + 1; |
|
|
61 n2 = 2 * n; |
|
|
62 x = u (n1:n2, 1:n) / u (1:n, 1:n); |
|
26
|
63 else |
|
55
|
64 error("usage: x = are (a, b, c)") |
|
26
|
65 endif |
|
|
66 |
|
|
67 endfunction |
