Mercurial > hg > octave-nkf
changeset 4611:c76a32c6f90c
[project @ 2003-11-14 17:48:46 by jwe]
author | jwe |
---|---|
date | Fri, 14 Nov 2003 17:48:46 +0000 |
parents | 02d2fcf835fc |
children | d44675070f1a |
files | scripts/ChangeLog scripts/control/base/dare.m scripts/control/base/dlqr.m scripts/control/system/is_detectable.m scripts/control/system/is_stabilizable.m |
diffstat | 5 files changed, 140 insertions(+), 116 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,5 +1,12 @@ 2003-11-14 Gabriele Pannocchia <g.pannocchia@ing.unipi.it> + * control/base/dare.m: + * control/base/dlqr.m: + + * control/system/is_detectable.m: Use Hautus Lemma. + Correct the behavior for discrete-time systems. + * control/system/is_stabilizable.m: Likewise. + * linear-algebra/krylov.m: Return H = [] in Vnrm == 0 case. * linear-algebra/krylovb.m: Fix typo in usage message.
--- a/scripts/control/base/dare.m +++ b/scripts/control/base/dare.m @@ -18,20 +18,20 @@ ## 02111-1307, USA. ## -*- texinfo -*- -## @deftypefn {Function File} {} dare (@var{a}, @var{b}, @var{c}, @var{r}, @var{opt}) +## @deftypefn {Function File} {} dare (@var{a}, @var{b}, @var{q}, @var{r}, @var{opt}) ## ## Return the solution, @var{x} of the discrete-time algebraic Riccati ## equation ## @iftex ## @tex ## $$ -## A^TXA - X + A^TXB (R + B^TXB)^{-1} B^TXA + C = 0 +## A^TXA - X + A^TXB (R + B^TXB)^{-1} B^TXA + Q = 0 ## $$ ## @end tex ## @end iftex ## @ifinfo ## @example -## a' x a - x + a' x b (r + b' x b)^(-1) b' x a + c = 0 +## a' x a - x + a' x b (r + b' x b)^(-1) b' x a + q = 0 ## @end example ## @end ifinfo ## @noindent @@ -44,9 +44,9 @@ ## @item b ## @var{n} by @var{m}. ## -## @item c +## @item q ## @var{n} by @var{n}, symmetric positive semidefinite, or @var{p} by @var{n}. -## In the latter case @math{c:=c'*c} is used. +## In the latter case @math{q:=q'*q} is used. ## ## @item r ## @var{m} by @var{m}, symmetric positive definite (invertible). @@ -74,7 +74,7 @@ ## Created: August 1993 ## Adapted-By: jwe -function x = dare (a, b, c, r, opt) +function x = dare (a, b, q, r, opt) if (nargin == 4 | nargin == 5) if (nargin == 5) @@ -86,28 +86,29 @@ opt = "B"; endif - ## dimension checks are done in is_controllable, is_observable - if (is_controllable (a, b) == 0) - warning ("dare: a,b are not controllable"); - elseif (is_observable (a, c) == 0) - warning ("dare: a,c are not observable"); + + if ((p = issquare (q)) == 0) + q = q'*q; endif - if ((p = issquare (c)) == 0) - c = c'*c; - p = rows (c); - endif + ##Checking positive definiteness + if (isdefinite(r)<=0) + error("dare: r not positive definite"); + end + if (isdefinite(q)<0) + error("dare: q not positive semidefinite"); + end + ## Check r dimensions. - n = rows(a); - m = columns(b); + [n,m] = size(b); if ((m1 = issquare (r)) == 0) - warning ("dare: r is not square"); + error ("dare: r is not square"); elseif (m1 != m) - warning ("b,r are not conformable"); + error ("b,r are not conformable"); endif - s1 = [a, zeros(n) ; -c, eye(n)]; + s1 = [a, zeros(n) ; -q, eye(n)]; s2 = [eye(n), (b/r)*b' ; zeros(n), a']; [c,d,s1,s2] = balance(s1,s2,opt); [aa,bb,u,lam] = qz(s1,s2,"S"); @@ -116,7 +117,7 @@ n2 = 2*n; x = u (n1:n2, 1:n)/u(1:n, 1:n); else - usage ("x = dare (a, b, c, r {,opt})"); + usage ("x = dare (a, b, q, r {,opt})"); endif endfunction
--- a/scripts/control/base/dlqr.m +++ b/scripts/control/base/dlqr.m @@ -101,8 +101,6 @@ ## Created: August 1993 ## Converted to discrete time by R. B. Tenison ## (btenison@eng.auburn.edu) October 1993 -## Modified by Gabriele Pannocchia <pannocchia@ing.unipi.it> -## July 2000 function [k, p, e] = dlqr (a, b, q, r, s) @@ -110,28 +108,10 @@ error ("dlqr: invalid number of arguments"); endif - ## Check a. - if ((n = issquare (a)) == 0) - error ("dlqr: requires 1st parameter(a) to be square"); - endif - - ## Check b. - [n1, m] = size (b); - if (n1 != n) - error ("dlqr: a,b not conformal"); - endif + ## Dimension check is done inside dare.m + [n,m] = size(b); - ## Check q. - if ((n1 = issquare (q)) == 0 || n1 != n) - error ("dlqr: q must be square and conformal with a"); - endif - - ## Check r. - if((m1 = issquare(r)) == 0 || m1 != m) - error ("dlqr: r must be square and conformal with column dimension of b"); - endif - - ## Check if n is there. + ## Check if s is there. if (nargin == 5) [n1, m1] = size (s); if (n1 != n || m1 != m) @@ -139,7 +119,6 @@ endif ## Incorporate cross term into a and q. - ao = a - (b/r)*s'; qo = q - (s/r)*s'; else @@ -148,15 +127,23 @@ qo = q; endif - ## Check that q, (r) are symmetric, positive (semi)definite - if (issymmetric (q) && issymmetric (r) - && all (eig (q) >= 0) && all (eig (r) > 0)) - p = dare (ao, b, qo, r); - k = (r+b'*p*b)\(b'*p*a + s'); - e = eig (a - b*k); - else - error ("dlqr: q (r) must be symmetric positive (semi) definite"); + ## Checking stabilizability and detectability (dimensions are checked inside these calls) + tol = 200*eps; + if (is_stabilizable (ao, b,tol,1) == 0) + error ("dlqr: (a,b) not stabilizable"); endif + dflag = is_detectable (ao, qo, tol,1); + if ( dflag == 0) + warning ("dlqr: (a,q) not detectable"); + elseif ( dflag == -1) + error("dlqr: (a,q) has non minimal modes near unit circle"); + end + + ## Compute the Riccati solution + p = dare (ao, b, qo, r); + k = (r+b'*p*b)\(b'*p*a + s'); + e = eig (a - b*k); + endfunction
--- a/scripts/control/system/is_detectable.m +++ b/scripts/control/system/is_detectable.m @@ -17,18 +17,17 @@ ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{retval}, @var{u}] =} is_detectable (@var{a}, @var{c}, @var{tol}) -## @deftypefnx {Function File} {[@var{retval}, @var{u}] =} is_detectable (@var{sys}, @var{tol}) -## Test for detactability (observability of unstable modes) of -## (@var{a},@var{c}). +## @deftypefn {Function File} {@var{retval} =} is_detectable (@var{a}, @var{c}, @var{tol}, @var{dflg}) +## @deftypefnx {Function File} {@var{retval} =} is_detectable (@var{sys}, @var{tol}) +## Test for detactability (observability of unstable modes) of (@var{a},@var{c}). ## ## Returns 1 if the system @var{a} or the pair (@var{a},@var{c})is -## detectable, 0 if not. +## detectable, 0 if not, and -1 if the system has unobservable modes at the +## imaginary axis (unit circle for discrete-time systems) ## ## @strong{See} @code{is_stabilizable} for detailed description of ## arguments and computational method. ## -## Default: tol = 10*norm(a,'fro')*eps ## ## @end deftypefn ## @seealso{is_stabilizable, size, rows, columns, length, ismatrix, @@ -38,27 +37,32 @@ ## Created: August 1993 ## Updated by John Ingram (ingraje@eng.auburn.edu) July 1996. -function [retval, U] = is_detectable (a, c, tol) +function [retval, U] = is_detectable (a, c, tol, dflg) if( nargin < 1) - usage("[retval,U] = is_detectable(a , c {, tol})"); + usage("retval = is_detectable(a , {c , tol, dlfg})"); elseif(isstruct(a)) ## system form if(nargin == 2) tol = c; elseif(nargin > 2) - usage("[retval,U] = is_detectable(sys {, tol})"); + usage("retval = is_detectable(sys {, tol})"); endif + dflg = is_digital(a); [a,b,c] = sys2ss(a); - elseif(nargin > 3) - usage("[retval,U] = is_detectable(a , c {, tol})"); - endif - if(exist("tol")) - [retval,U] = is_stabilizable (a', c', tol); else - [retval,U] = is_stabilizable (a', c'); - endif + if ((nargin > 4)||(nargin == 1)) + usage("retval = is_detectable(a , {c , tol, dflg})"); + endif + if (~exist("dflg")) + dflg = 0; + end + end + if(~exist("tol")) + tol = 200*eps; + end + retval = is_stabilizable(a',c',tol,dflg); endfunction
--- a/scripts/control/system/is_stabilizable.m +++ b/scripts/control/system/is_stabilizable.m @@ -1,5 +1,3 @@ -## Copyright (C) 1993, 1994, 1995 Auburn University. All rights reserved. -## ## This file is part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it @@ -17,25 +15,16 @@ ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{retval}, @var{u}] =} is_stabilizable (@var{sys}, @var{tol}) -## @deftypefnx {Function File} {[@var{retval}, @var{u}] =} is_stabilizable (@var{a}, @var{b}, @var{tol}) -## Logical check for system stabilizability (i.e., all unstable modes are controllable). -## -## Test for stabilizability is performed via an ordered Schur decomposition -## that reveals the unstable subspace of the system @var{a} matrix. +## @deftypefn {Function File} {@var{retval} =} is_stabilizable (@var{sys}, @var{tol}) +## @deftypefnx {Function File} {@var{retval} =} is_stabilizable (@var{a}, @var{b}, @var{tol}, @var{dflg}) +## Logical check for system stabilizability (i.e., all unstable modes are controllable). +## Returns 1 if the system is stabilizable, 0 if the the system is not stabilizable, -1 +## if the system has non stabilizable modes at the imaginary axis (unit circle for discrete-time +## systems. ## -## Returns @code{retval} = 1 if the system, @var{a}, is stabilizable, -## if the pair (@var{a}, @var{b}) is stabilizable, or 0 if not. -## @var{u} = orthogonal basis of controllable subspace. +## Test for stabilizability is performed via Hautus Lemma. If @var{dflg}!=0 assume that +## discrete-time matrices (a,b) are supplied. ## -## Controllable subspace is determined by applying Arnoldi iteration with -## complete re-orthogonalization to obtain an orthogonal basis of the -## Krylov subspace. -## @example -## span ([b,a*b,...,a^ b]). -## @end example -## tol is a roundoff paramter, set to 200*eps if omitted. -## @end deftypefn ## See also: size, rows, columns, length, ismatrix, isscalar, isvector ## is_observable, is_stabilizable, is_detectable @@ -45,45 +34,81 @@ ## Updated by A. S. Hodel (scotte@eng.auburn.edu) Aubust, 1995 to use krylovb ## Updated by John Ingram (ingraje@eng.auburn.edu) July, 1996 to accept systems -function [retval, U] = is_stabilizable (a, b, tol) +function retval = is_stabilizable (a, b, tol, dflg) - if(nargin < 1) usage("[retval,U] = is_stabilizable(a {, b ,tol})"); + if(nargin < 1) + usage("[retval,U] = is_stabilizable(a {, b ,tol, dflg})"); elseif(isstruct(a)) - ## sustem passed. + ## system passed. if(nargin == 2) tol = b; % get tolerance elseif(nargin > 2) - usage("[retval,U] = is_stabilizable(sys{,tol})"); + usage("retval = is_stabilizable(sys{,tol})"); endif + disc = is_digital(a); [a,b] = sys2ss(a); else ## a,b arguments sent directly. - if(nargin > 3) - usage("[retval,U] = is_stabilizable(a {, b ,tol})"); + if ((nargin > 4)||(nargin == 1)) + usage("retval = is_stabilizable(a {, b ,tol, dflg})"); endif - endif - - if(exist("tol")) - [retval,U] = is_controllable(a,b,tol); - else - [retval,U] = is_controllable(a,b); - tol = 1e2*rows(b)*eps; + if(exist("dflg")) + disc = dflg; + else + disc = 0; + end endif - if( !retval & columns(U) > 0) - ## now use an ordered Schur decomposition to get an orthogonal - ## basis of the unstable subspace... - n = rows(a); - [ua, s] = schur (-(a+eye(n)*tol), "A"); - k = sum( real(eig(a)) >= 0 ); # count unstable poles + if(~exist("tol")) + tol = 200*eps; + end + - if( k > 0 ) - ua = ua(:,1:k); - ## now see if span(ua) is contained in span(U) - retval = (norm(ua - U*U'*ua) < tol); + ## Checking dimensions + n = is_square(a); + if (n==0) + error("is_stabilizable: a must be square"); + end + [nr,m] = size(b); + if (nr!=n) + error("is_stabilizable: (a,b) not conformal"); + end + + ##Computing the eigenvalue of A + L = eig(a); + retval = 1; + specflag = 0; + for i=1:n + if (disc==0) + ## Continuous time case + rL = real(L(i)); + if (rL>=0) + H = [eye(n)*L(i)-a, b]; + f = (rank(H,tol)==n); + if (f==0) + retval = 0; + if (rL==0) + specflag = 1; + end + end + end else - retval = 1; # all poles stable - endif - endif - -endfunction + ## Discrete time case + rL = abs(L(i)); + if (rL>=1) + H = [eye(n)*L(i)-a, b]; + f = (rank(H,tol)==n); + if (f==0) + retval = 0; + if (rL==1) + specflag = 1; + end + end + end + end + end + if (specflag==1) + ## This means that the system has uncontrollable modes at the imaginary axis + ## (or at the unit circle for discrete time systems) + retval = -1; + end