Mercurial > hg > octave-terminal
changeset 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 | 1d4cfd89ebb6 |
children | c63068b25a63 |
files | scripts/control/are.m scripts/control/is_controllable.m scripts/control/is_observable.m scripts/general/is_square.m |
diffstat | 4 files changed, 146 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/scripts/control/are.m @@ -0,0 +1,60 @@ +function x = are(a,b,c,opt) + +# +# usage: x = are(a,b,c{,opt}) +# +# solves algebraic riccati equation +# +# a' x + x a - x b x + c = 0 +# +# for identically dimensioned square matrices a,b,c. If b (c) is not square, +# then the function attempts to use b*b' (c'*c) instead. +# +# opt is an option passed to the eigenvalue balancing routine; default is 'B' +# +# see also: balance + +if((nargin == 3) || (nargin == 4)) + if(nargin==4) + if((opt ~= 'N') || (opt ~= 'P') || (opt ~= 'S') || (opt ~= 'B')) + printf('warning: are: opt has an illegal value; setting to B'); + opt = 'B'; + endif + else + opt = 'B'; + endif + if( (n = is_square(a) ) == 0) + error('are: a is not square') + endif + if(is_controllable(a,b) == 0) + printf('warning: are: a,b are not controllable') + endif + if( (m = is_square(b)) == 0) + b = b*b'; + m = rows(b); + endif + if(is_observable(a,c) == 0) + printf('warning: are: a,c are not observable') + endif + if( (p = is_square(c)) == 0) + c = c'*c; + p = rows(c); + endif + if( (n ~= m) || (n ~= p)) + error('are: a,b,c not conformably dimensioned.') + endif + + # should check for controllability/observability here + # use Boley-Golub (Syst. Contr. Letters, 1984) method, not the + # n-1 + # rank([ B A*B ... A^ *B]) method + + [u,s] = schur(balance([a,-b;-c,-a'],opt),'A'); + n1=n+1; + n2 = 2*n; + x = u(n1:n2,1:n)/u(1:n,1:n); +else + error('usage: x = are(a,b,c)') +endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/is_controllable.m @@ -0,0 +1,47 @@ +function retval = is_controllable (a,b,tol) + +# usage: is_controllable (a,b{,tol}) +# +# returns 1 the pair(a,b) is controllable, then return value is the +# dimension of x. 0therwise, returns a value of 0 +# +# See also: size, rows, columns, length, is_matrix, is_scalar, is_vector + +# This should really use the method below, but I'm being lazy for now: +# +# Controllability is determined by applying Arnoldi iteration with complete +# re-orthogonalization to obtain an orthogonal basis of the Krylov subspace +# n-1 +# span([b,a*b,...,a^ b]). tol is a roundoff paramter, +# set to 2*eps if omitted + + if ((nargin == 2) || (nargin == 3)) + n = is_square(a); + [nr, nc] = size (b); + if((n == 0) || (n != nr) || (nc == 0)) + retval = 0; + else + + m = b; + tmp = b; + for ii=1:(n-1) + tmp = a*tmp; + m = [m,tmp]; + endfor + + # if n is of any significant size, m will be low rank, so be careful! + if(nargin == 3) + if(is_scalar(tol)) + retval = (rank(m,tol) == n); + else + error('is_controllable: tol must be a scalar') + endif + else + retval = (rank(m) == n); + endif + endif + else + error ("usage: is_controllable (a,b)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/is_observable.m @@ -0,0 +1,18 @@ +function retval = is_observable (a,c,tol) + +# usage: is_observable (a,c{,tol}) +# +# returns 1 the pair(a,c) is observable, then return value is the +# dimension of x. 0therwise, returns a value of 0 +# +# See also: size, rows, columns, length, is_matrix, is_scalar, is_vector + + if (nargin == 2) + retval = is_controllable(a',c'); + elseif (nargin == 3) + retval = iscontrollable(a',c',tol); + else + error ("usage: is_observable (a,c{,tol})"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/is_square.m @@ -0,0 +1,21 @@ +function retval = is_square (x) + +# usage: is_square (x) +# +# If x is square, then return value is the dimension of x. +# otherwise, returns a value of 0 +# +# See also: size, rows, columns, length, is_matrix, is_scalar, is_vector + + if (nargin == 1) + [nr, nc] = size (x); + if( nr == nc) + retval = nr; + else + retval = 0; + endif + else + error ("usage: is_square (x)"); + endif + +endfunction