From bug-octave-request at bevo dot che dot wisc dot edu Thu Nov 6 10:44:36 2003 Subject: bugs and patches for is_stabilizable.m/is_detectable.m From: Gabriele Pannocchia To: bug-octave at bevo dot che dot wisc dot edu Date: Thu, 6 Nov 2003 17:43:48 +0100 Hi John, the functions is_stabilizable.m and is_detectable.m are meant to work for both continuous and discrete time systems. However in the latter case, they may return wrong answers. Here is an example of a non-stabilizable discrete-time system. octave:22> a=diag([-2,-0.5]);b=[0;1];c=[1,0];d=0; octave:23> sys=ss2sys(a,b,c,d,1); octave:24> is_stabilizable (sys) ans = 1 As I anticipated to you, I rewrote these two functions using a different algorithm (Hautus lemma). Patches follow. --- ../octave/scripts/control/system/is_stabilizable.m Fri Aug 9 20:58:14 2002 +++ is_stabilizable.m Thu Nov 6 15:36:42 2003 at @ -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 at @ -17,25 +15,16 @@ ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. ## -*- texinfo -*- -## at deftypefn {Function File} {[@var{retval}, @var{u}] =} is_stabilizable (@var{sys}, @var{tol}) -## at 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 at var{a} matrix. +## at deftypefn {Function File} {@var{retval} =} is_stabilizable (@var{sys}, @var{tol}) +## at 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 at code{retval} = 1 if the system, @var{a}, is stabilizable, -## if the pair ( at var{a}, @var{b}) is stabilizable, or 0 if not. -## at var{u} = orthogonal basis of controllable subspace. +## Test for stabilizability is performed via Hautus Lemma. If at 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. -## at example -## span ([b,a*b,...,a^ b]). -## at end example -## tol is a roundoff paramter, set to 200*eps if omitted. -## at end deftypefn ## See also: size, rows, columns, length, ismatrix, isscalar, isvector ## is_observable, is_stabilizable, is_detectable at @ -44,46 +33,84 @@ ## Created: August 1993 ## Updated by A. S. Hodel (scotte at eng dot auburn dot edu) Aubust, 1995 to use krylovb ## Updated by John Ingram (ingraje at eng dot auburn dot edu) July, 1996 to accept systems +## Modified by Gabriele Pannocchia , November 2003 to use Hautus +## Lemma and to correct the behavior for discrete-time 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; - 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( k > 0 ) - ua = ua(:,1:k); - ## now see if span(ua) is contained in span(U) - retval = (norm(ua - U*U'*ua) < tol); + if(exist("dflg")) + disc = dflg; else - retval = 1; # all poles stable - endif + disc = 0; + end endif -endfunction + if(~exist("tol")) + tol = 200*eps; + end + + + ## 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 + ## 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 Using this path, results are the followings: octave:95> a=diag([-2,-0.5]);b=[0;1];c=[1,0];d=0; octave:96> sys=ss2sys(a,b,c,d,1); octave:97> is_stabilizable (sys) ans = 0 Due to this patch I also modified is_detectable.m: --- ../octave/scripts/control/system/is_detectable.m Fri Aug 9 20:58:14 2002 +++ is_detectable.m Thu Nov 6 15:27:56 2003 at @ -17,18 +17,17 @@ ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. ## -*- texinfo -*- -## at deftypefn {Function File} {[@var{retval}, @var{u}] =} is_detectable (@var{a}, @var{c}, @var{tol}) -## at deftypefnx {Function File} {[@var{retval}, @var{u}] =} is_detectable (@var{sys}, @var{tol}) -## Test for detactability (observability of unstable modes) of -## ( at var{a},@var{c}) dot +## at deftypefn {Function File} {@var{retval} =} is_detectable (@var{a}, @var{c}, @var{tol}, @var{dflg}) +## at deftypefnx {Function File} {@var{retval} =} is_detectable (@var{sys}, @var{tol}) +## Test for detactability (observability of unstable modes) of ( at var{a},@var{c}) dot ## ## Returns 1 if the system at 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) ## ## at strong{See} @code{is_stabilizable} for detailed description of ## arguments and computational method. ## -## Default: tol = 10*norm(a,'fro')*eps ## ## at end deftypefn ## at seealso{is_stabilizable, size, rows, columns, length, ismatrix, at @ -37,28 +36,35 @@ ## Author: A. S. Hodel ## Created: August 1993 ## Updated by John Ingram (ingraje at eng dot auburn dot edu) July 1996. +## Modified by Gabriele Pannocchia , November 2003 to use Hautus +## Lemma and to correct the behavior for discrete-time systems. -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 Gabriele ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------