From bug-octave-request at bevo dot che dot wisc dot edu Thu Nov 6 10:51:15 2003 Subject: changes to dlqr.m dare.m From: Gabriele Pannocchia To: bug-octave at bevo dot che dot wisc dot edu Date: Thu, 6 Nov 2003 17:50:24 +0100 John, changes to dlqr.m and dare.m follows. These changes are made based on the discussions I had with Jim about the behavior of these two functions. Notice that these modified functions rely on the patches I just sent for is_stabilizable.m and is_detectable.m and on a new function, isdefinite.m, also attached below. --- ../octave/scripts/control/base/dare.m Fri Aug 9 20:58:13 2002 +++ dare.m Thu Nov 6 15:45:57 2003 at @ -18,20 +18,20 @@ ## 02111-1307, USA. ## -*- texinfo -*- -## at deftypefn {Function File} {} dare (@var{a}, @var{b}, @var{c}, @var{r}, @var{opt}) +## at deftypefn {Function File} {} dare (@var{a}, @var{b}, @var{q}, @var{r}, @var{opt}) ## ## Return the solution, at var{x} of the discrete-time algebraic Riccati ## equation ## at iftex ## at 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 ## $$ ## at end tex ## at end iftex ## at ifinfo ## at 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 ## at end example ## at end ifinfo ## at noindent at @ -44,7 +44,7 @@ ## at item b ## at var{n} by @var{m}. ## -## at item c +## at item q ## at var{n} by @var{n}, symmetric positive semidefinite, or @var{p} by @var{n}. ## In the latter case at math{c:=c'*c} is used. ## at @ -73,8 +73,9 @@ ## Author: A. S. Hodel ## Created: August 1993 ## Adapted-By: jwe +## Modified by: Gabriele Pannocchia , November 2003 -function x = dare (a, b, c, r, opt) +function x = dare (a, b, q, r, opt) if (nargin == 4 | nargin == 5) if (nargin == 5) at @ -86,28 +87,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"); at @ -116,7 +118,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 --- ../octave/scripts/control/base/dlqr.m Fri Aug 9 20:58:13 2002 +++ dlqr.m Thu Nov 6 15:56:02 2003 at @ -101,8 +101,8 @@ ## Created: August 1993 ## Converted to discrete time by R. B. Tenison ## (btenison at eng dot auburn dot edu) October 1993 -## Modified by Gabriele Pannocchia -## July 2000 +## Modified by Gabriele Pannocchia +## July 2000, November 2003 function [k, p, e] = dlqr (a, b, q, r, s) at @ -110,28 +110,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) at @ -139,7 +121,6 @@ endif ## Incorporate cross term into a and q. - ao = a - (b/r)*s'; qo = q - (s/r)*s'; else at @ -148,15 +129,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"); - endif + ## 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 NEW FUNCTIONS isdefinite.m ## -*- texinfo -*- ## at deftypefn {Function File} {} isdefinite (@var{x},@var{tol}) ## If at var{x} is symmetric positive definite within the tolerance specified by @var{tol}, ## then return 1. If at var{x} is symmetric positive semidefinite, then return 0. ## Otherwise, return -1. If at var{tol} is omitted, use a tolerance equal to 100 times ## the machine precision. ## at end deftypefn ## at seealso{size, rows, columns, length, ismatrix, issymmetric, ## issquare, and isvector} ## Created by: Gabriele Pannocchia , November 2003 function retval = isdefinite (x,tol) if ((nargin == 1)||(nargin == 2)) if (nargin==1) tol = 100*eps; endif sym = issymmetric(x,tol); if (sym>0) ##Matrix is symmetric. Checking eigenvalues mineig = min(eig(x)); if (mineig>tol) retval = 1; elseif (mineig>-tol) retval = 0; else retval = -1; end else error("isdefinite: matrix x must be symmetric"); endif else usage ("isdefinite (x,tol)"); endif 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 -------------------------------------------------------------