Mercurial > octave-nkf
comparison scripts/control/are.m @ 26:e90ea9cbd4de
[project @ 1993-08-10 20:56:55 by jwe]
Initial revision
| author | jwe |
|---|---|
| date | Tue, 10 Aug 1993 20:56:55 +0000 |
| parents | |
| children | b973bf9a9dba |
comparison
equal
deleted
inserted
replaced
| 25:1d4cfd89ebb6 | 26:e90ea9cbd4de |
|---|---|
| 1 function x = are(a,b,c,opt) | |
| 2 | |
| 3 # | |
| 4 # usage: x = are(a,b,c{,opt}) | |
| 5 # | |
| 6 # solves algebraic riccati equation | |
| 7 # | |
| 8 # a' x + x a - x b x + c = 0 | |
| 9 # | |
| 10 # for identically dimensioned square matrices a,b,c. If b (c) is not square, | |
| 11 # then the function attempts to use b*b' (c'*c) instead. | |
| 12 # | |
| 13 # opt is an option passed to the eigenvalue balancing routine; default is 'B' | |
| 14 # | |
| 15 # see also: balance | |
| 16 | |
| 17 if((nargin == 3) || (nargin == 4)) | |
| 18 if(nargin==4) | |
| 19 if((opt ~= 'N') || (opt ~= 'P') || (opt ~= 'S') || (opt ~= 'B')) | |
| 20 printf('warning: are: opt has an illegal value; setting to B'); | |
| 21 opt = 'B'; | |
| 22 endif | |
| 23 else | |
| 24 opt = 'B'; | |
| 25 endif | |
| 26 if( (n = is_square(a) ) == 0) | |
| 27 error('are: a is not square') | |
| 28 endif | |
| 29 if(is_controllable(a,b) == 0) | |
| 30 printf('warning: are: a,b are not controllable') | |
| 31 endif | |
| 32 if( (m = is_square(b)) == 0) | |
| 33 b = b*b'; | |
| 34 m = rows(b); | |
| 35 endif | |
| 36 if(is_observable(a,c) == 0) | |
| 37 printf('warning: are: a,c are not observable') | |
| 38 endif | |
| 39 if( (p = is_square(c)) == 0) | |
| 40 c = c'*c; | |
| 41 p = rows(c); | |
| 42 endif | |
| 43 if( (n ~= m) || (n ~= p)) | |
| 44 error('are: a,b,c not conformably dimensioned.') | |
| 45 endif | |
| 46 | |
| 47 # should check for controllability/observability here | |
| 48 # use Boley-Golub (Syst. Contr. Letters, 1984) method, not the | |
| 49 # n-1 | |
| 50 # rank([ B A*B ... A^ *B]) method | |
| 51 | |
| 52 [u,s] = schur(balance([a,-b;-c,-a'],opt),'A'); | |
| 53 n1=n+1; | |
| 54 n2 = 2*n; | |
| 55 x = u(n1:n2,1:n)/u(1:n,1:n); | |
| 56 else | |
| 57 error('usage: x = are(a,b,c)') | |
| 58 endif | |
| 59 | |
| 60 endfunction |