|
55
|
1 function x = are (a, b, c, opt) |
|
26
|
2 |
|
55
|
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 # |
|
|
15 # see also: balance |
|
|
16 |
|
55
|
17 if (nargin == 3 || nargin == 4) |
|
|
18 if (nargin == 4) |
|
|
19 if (! (strcmp (opt, "N") || strcmp (opt, "P") ... |
|
|
20 || strcmp (opt, "S") || strcmp (opt, "B") ... |
|
|
21 || strcmp (opt, "n") || strcmp (opt, "p") ... |
|
|
22 || strcmp (opt, "s") || strcmp (opt, "b"))) |
|
|
23 printf ("warning: are: opt has an illegal value; setting to B"); |
|
|
24 opt = "B"; |
|
|
25 endif |
|
|
26 else |
|
|
27 opt = "B"; |
|
|
28 endif |
|
|
29 if ((n = is_square(a)) == 0) |
|
|
30 error ("are: a is not square"); |
|
|
31 endif |
|
|
32 |
|
|
33 if (is_controllable(a,b) == 0) |
|
|
34 printf("warning: are: a, b are not controllable"); |
|
26
|
35 endif |
|
55
|
36 if ((m = is_square (b)) == 0) |
|
|
37 b = b * b'; |
|
|
38 m = rows (b); |
|
|
39 endif |
|
|
40 if (is_observable (a, c) == 0) |
|
|
41 printf ("warning: are: a,c are not observable"); |
|
|
42 endif |
|
|
43 if ((p = is_square (c)) == 0) |
|
|
44 c = c' * c; |
|
|
45 p = rows (c); |
|
|
46 endif |
|
|
47 if (n != m || n != p) |
|
|
48 error ("are: a, b, c not conformably dimensioned."); |
|
|
49 endif |
|
|
50 |
|
|
51 # Should check for controllability/observability here |
|
|
52 # use Boley-Golub (Syst. Contr. Letters, 1984) method, not the |
|
|
53 # |
|
|
54 # n-1 |
|
|
55 # rank ([ B A*B ... A^ *B]) method |
|
|
56 |
|
|
57 [u, s] = schur (balance ([a, -b; -c, -a'], opt), "A"); |
|
|
58 n1 = n + 1; |
|
|
59 n2 = 2 * n; |
|
|
60 x = u (n1:n2, 1:n) / u (1:n, 1:n); |
|
26
|
61 else |
|
55
|
62 error("usage: x = are (a, b, c)") |
|
26
|
63 endif |
|
|
64 |
|
|
65 endfunction |
