# HG changeset patch
# User jwe
# Date 1068832126 0
# Node ID c76a32c6f90cda4b078f77306078116232e8c228
# Parent  02d2fcf835fc7c9f587a8937fc8033dd947731a9
[project @ 2003-11-14 17:48:46 by jwe]

diff -r 02d2fcf835fc -r c76a32c6f90c scripts/ChangeLog
--- a/scripts/ChangeLog	Fri Nov 14 17:08:59 2003 +0000
+++ b/scripts/ChangeLog	Fri Nov 14 17:48:46 2003 +0000
@@ -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.
diff -r 02d2fcf835fc -r c76a32c6f90c scripts/control/base/dare.m
--- a/scripts/control/base/dare.m	Fri Nov 14 17:08:59 2003 +0000
+++ b/scripts/control/base/dare.m	Fri Nov 14 17:48:46 2003 +0000
@@ -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
diff -r 02d2fcf835fc -r c76a32c6f90c scripts/control/base/dlqr.m
--- a/scripts/control/base/dlqr.m	Fri Nov 14 17:08:59 2003 +0000
+++ b/scripts/control/base/dlqr.m	Fri Nov 14 17:48:46 2003 +0000
@@ -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
 
diff -r 02d2fcf835fc -r c76a32c6f90c scripts/control/system/is_detectable.m
--- a/scripts/control/system/is_detectable.m	Fri Nov 14 17:08:59 2003 +0000
+++ b/scripts/control/system/is_detectable.m	Fri Nov 14 17:48:46 2003 +0000
@@ -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
 
diff -r 02d2fcf835fc -r c76a32c6f90c scripts/control/system/is_stabilizable.m
--- a/scripts/control/system/is_stabilizable.m	Fri Nov 14 17:08:59 2003 +0000
+++ b/scripts/control/system/is_stabilizable.m	Fri Nov 14 17:48:46 2003 +0000
@@ -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