Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[7710] trunk/octave-forge/main/control
|
From: <par...@us...> - 2010-09-13 22:59:34
|
Revision: 7710
http://octave.svn.sourceforge.net/octave/?rev=7710&view=rev
Author: paramaniac
Date: 2010-09-13 22:59:27 +0000 (Mon, 13 Sep 2010)
Log Message:
-----------
control: various small enhancements to ARE solvers
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/care.m
trunk/octave-forge/main/control/inst/dare.m
trunk/octave-forge/main/control/src/slsb02od.cc
Modified: trunk/octave-forge/main/control/inst/care.m
===================================================================
--- trunk/octave-forge/main/control/inst/care.m 2010-09-13 19:32:31 UTC (rev 7709)
+++ trunk/octave-forge/main/control/inst/care.m 2010-09-13 22:59:27 UTC (rev 7710)
@@ -25,25 +25,25 @@
## @strong{Inputs}
## @table @var
## @item a
-## Real matrix.
+## Real matrix (n-by-n).
## @item b
-## Real matrix.
+## Real matrix (n-by-m).
## @item q
-## Real matrix.
+## Real matrix (n-by-n).
## @item r
-## Real matrix.
+## Real matrix (m-by-m).
## @item s
-## Optional real matrix. If @var{s} is not specified, a zero matrix is assumed.
+## Optional real matrix (n-by-m). If @var{s} is not specified, a zero matrix is assumed.
## @end table
##
## @strong{Outputs}
## @table @var
## @item x
-## Unique stabilizing solution of the continuous-time Riccati equation.
+## Unique stabilizing solution of the continuous-time Riccati equation (n-by-n).
## @item l
-## Closed-loop poles.
+## Closed-loop poles (n-by-1).
## @item g
-## Corresponding gain matrix.
+## Corresponding gain matrix (m-by-n).
## @end table
##
## @example
@@ -68,39 +68,41 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: November 2009
-## Version: 0.3
+## Version: 0.4
function [x, l, g] = care (a, b, q, r, s = [])
## TODO: Add SLICOT SG02AD (Solution of continuous- or discrete-time
## algebraic Riccati equations for descriptor systems)
+ ## TODO: extract feedback matrix g from SB02OD (and SG02AD)
+
if (nargin < 4 || nargin > 5)
print_usage ();
endif
- if (! issquare (a))
- error ("care: a is not square");
+ if (! isreal (a) || ! issquare (a))
+ error ("care: a must be real and square");
endif
+
+ if (! isreal (b) || rows (a) != rows (b))
+ error ("care: b must be real and conformal to a");
+ endif
- if (! issquare (q))
- error ("care: q is not square");
+ if (! isreal (q) || ! issquare (q))
+ error ("care: q must be real and square");
endif
- if (! issquare (r))
- error ("care: r is not square");
+ if (! isreal (r) || ! issquare (r))
+ error ("care: r must be real and square");
endif
- if (rows (a) != rows (b))
- error ("care: (a, b) not conformable");
- endif
-
if (columns (r) != columns (b))
error ("care: (b, r) not conformable");
endif
- if (! isempty (s) && any (size (s) != size (b)))
- error ("care: s (%dx%d) must be identically dimensioned with b (%dx%d)",
+ if (! isempty (s) && (! isreal (s) || any (size (s) != size (b))))
+ error ("care: s(%dx%d) must be real and identically dimensioned with b(%dx%d)",
rows (s), columns (s), rows (b), columns (b));
endif
@@ -124,24 +126,15 @@
## solve the riccati equation
if (isempty (s))
- ## unique stabilizing solution
- x = slsb02od (a, b, q, r, b, false, false);
-
- ## corresponding gain matrix
- g = r \ (b.'*x);
+ [x, l] = slsb02od (a, b, q, r, b, false, false);
+
+ g = r \ (b.'*x); # gain matrix
else
- ## unique stabilizing solution
- x = slsb02od (a, b, q, r, s, false, true);
-
- ## corresponding gain matrix
- g = r \ (b.'*x + s.');
+ [x, l] = slsb02od (a, b, q, r, s, false, true);
+
+ g = r \ (b.'*x + s.'); # gain matrix
endif
- ## closed-loop poles
- l = eig (a - b*g);
-
- ## TODO: use alphar, alphai and beta from SB02OD
-
endfunction
Modified: trunk/octave-forge/main/control/inst/dare.m
===================================================================
--- trunk/octave-forge/main/control/inst/dare.m 2010-09-13 19:32:31 UTC (rev 7709)
+++ trunk/octave-forge/main/control/inst/dare.m 2010-09-13 22:59:27 UTC (rev 7710)
@@ -25,25 +25,25 @@
## @strong{Inputs}
## @table @var
## @item a
-## Real matrix.
+## Real matrix (n-by-n).
## @item b
-## Real matrix.
+## Real matrix (n-by-m).
## @item q
-## Real matrix.
+## Real matrix (n-by-n).
## @item r
-## Real matrix.
+## Real matrix (m-by-m).
## @item s
-## Optional real matrix. If @var{s} is not specified, a zero matrix is assumed.
+## Optional real matrix (n-by-m). If @var{s} is not specified, a zero matrix is assumed.
## @end table
##
## @strong{Outputs}
## @table @var
## @item x
-## Unique stabilizing solution of the discrete-time Riccati equation.
+## Unique stabilizing solution of the discrete-time Riccati equation (n-by-n).
## @item l
-## Closed-loop poles.
+## Closed-loop poles (n-by-1).
## @item g
-## Corresponding gain matrix.
+## Corresponding gain matrix (m-by-n).
## @end table
##
## @example
@@ -68,39 +68,41 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: October 2009
-## Version: 0.3
+## Version: 0.4
function [x, l, g] = dare (a, b, q, r, s = [])
## TODO: Add SLICOT SG02AD (Solution of continuous- or discrete-time
## algebraic Riccati equations for descriptor systems)
+ ## TODO: extract feedback matrix g from SB02OD (and SG02AD)
+
if (nargin < 4 || nargin > 5)
print_usage ();
endif
- if (! issquare (a))
- error ("dare: a is not square");
+ if (! isreal (a) || ! issquare (a))
+ error ("dare: a must be real and square");
endif
+
+ if (! isreal (b) || rows (a) != rows (b))
+ error ("dare: b must be real and conformal to a");
+ endif
- if (! issquare (q))
- error ("dare: q is not square");
+ if (! isreal (q) || ! issquare (q))
+ error ("dare: q must be real and square");
endif
- if (! issquare (r))
- error ("dare: r is not square");
+ if (! isreal (r) || ! issquare (r))
+ error ("dare: r must be real and square");
endif
- if (rows (a) != rows (b))
- error ("dare: (a, b) not conformable");
- endif
-
if (columns (r) != columns (b))
error ("dare: (b, r) not conformable");
endif
- if (! isempty (s) && any (size (s) != size (b)))
- error ("dare: s (%dx%d) must be identically dimensioned with b (%dx%d)",
+ if (! isempty (s) && (! isreal (s) || any (size (s) != size (b))))
+ error ("dare: s(%dx%d) must be real and identically dimensioned with b(%dx%d)",
rows (s), columns (s), rows (b), columns (b));
endif
@@ -124,24 +126,15 @@
## solve the riccati equation
if (isempty (s))
- ## unique stabilizing solution
- x = slsb02od (a, b, q, r, b, true, false);
+ [x, l] = slsb02od (a, b, q, r, b, true, false);
- ## corresponding gain matrix
- g = (r + b.'*x*b) \ (b.'*x*a);
+ g = (r + b.'*x*b) \ (b.'*x*a); # gain matrix
else
- ## unique stabilizing solution
- x = slsb02od (a, b, q, r, s, true, true);
-
- ## corresponding gain matrix
- g = (r + b.'*x*b) \ (b.'*x*a + s.');
+ [x, l] = slsb02od (a, b, q, r, s, true, true);
+
+ g = (r + b.'*x*b) \ (b.'*x*a + s.'); # gain matrix
endif
- ## closed-loop poles
- l = eig (a - b*g);
-
- ## TODO: use alphar, alphai and beta from SB02OD
-
endfunction
@@ -167,5 +160,5 @@
%! ge = [ 0.4092 1.7283];
%!
%!assert (x, xe, 1e-4);
-%!assert (l, le, 1e-4);
+%!assert (sort (l), sort (le), 1e-4);
%!assert (g, ge, 1e-4);
\ No newline at end of file
Modified: trunk/octave-forge/main/control/src/slsb02od.cc
===================================================================
--- trunk/octave-forge/main/control/src/slsb02od.cc 2010-09-13 19:32:31 UTC (rev 7709)
+++ trunk/octave-forge/main/control/src/slsb02od.cc 2010-09-13 22:59:27 UTC (rev 7710)
@@ -23,13 +23,14 @@
Author: Lukas Reichlin <luk...@gm...>
Created: February 2010
-Version: 0.2
+Version: 0.2.1
*/
#include <octave/oct.h>
#include <f77-fcn.h>
#include "common.cc"
+#include <complex>
extern "C"
{
@@ -56,12 +57,12 @@
bool* BWORK,
int& INFO);
}
-
+
DEFUN_DLD (slsb02od, args, nargout, "Slicot SB02OD Release 5.0")
{
int nargin = args.length ();
octave_value_list retval;
-
+
if (nargin != 7)
{
print_usage ();
@@ -75,7 +76,7 @@
char uplo = 'U';
char jobl;
char sort = 'S';
-
+
Matrix a = args(0).matrix_value ();
Matrix b = args(1).matrix_value ();
Matrix q = args(2).matrix_value ();
@@ -88,54 +89,55 @@
dico = 'C';
else
dico = 'D';
-
+
if (ijobl == 0)
jobl = 'Z';
else
jobl = 'N';
-
+
int n = a.rows (); // n: number of states
int m = b.columns (); // m: number of inputs
int p = 0; // p: number of outputs, not used because FACT = 'N'
-
+
int lda = max (1, n);
int ldb = max (1, n);
int ldq = max (1, n);
int ldr = max (1, m);
int ldl = max (1, n);
-
+
// arguments out
double rcond;
-
+
int ldx = max (1, n);
Matrix x (ldx, n);
-
+
+ int nu = 2*n;
+ ColumnVector alfar (nu);
+ ColumnVector alfai (nu);
+ ColumnVector beta (nu);
+
+ int lds = max (1, 2*n + m);
+ Matrix s (lds, lds);
+
// unused output arguments
- OCTAVE_LOCAL_BUFFER (double, alfar, 2*n);
- OCTAVE_LOCAL_BUFFER (double, alfai, 2*n);
- OCTAVE_LOCAL_BUFFER (double, beta, 2*n);
-
- int lds = max (1, 2*n + m);
- OCTAVE_LOCAL_BUFFER (double, s, lds*lds);
-
int ldt = max (1, 2*n + m);
OCTAVE_LOCAL_BUFFER (double, t, ldt * 2*n);
-
+
int ldu = max (1, 2*n);
OCTAVE_LOCAL_BUFFER (double, u, ldu * 2*n);
-
+
// tolerance
double tol = 0; // use default value
-
+
// workspace
int liwork = max (1, m, 2*n);
OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
-
+
int ldwork = max (7*(2*n + 1) + 16, 16*n, 2*n + m, 3*m);
OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
-
+
OCTAVE_LOCAL_BUFFER (bool, bwork, 2*n);
-
+
// error indicator
int info;
@@ -153,9 +155,9 @@
l.fortran_vec (), ldl,
rcond,
x.fortran_vec (), ldx,
- alfar, alfai,
- beta,
- s, lds,
+ alfar.fortran_vec (), alfai.fortran_vec (),
+ beta.fortran_vec (),
+ s.fortran_vec (), lds,
t, ldt,
u, ldu,
tol,
@@ -166,13 +168,30 @@
if (f77_exception_encountered)
error ("are: slsb02od: exception in SLICOT subroutine SB02OD");
-
+
if (info != 0)
error ("are: slsb02od: SB02OD returned info = %d", info);
+
+ // assemble complex vector - adapted from DEFUN complex in data.cc
+ alfar.resize (n);
+ alfai.resize (n);
+ beta.resize (n);
+
+ ColumnVector zeror (n);
+ ColumnVector zeroi (n);
+
+ zeror = quotient (alfar, beta);
+ zeroi = quotient (alfai, beta);
+ ComplexColumnVector zero (n, Complex ());
+
+ for (octave_idx_type i = 0; i < n; i++)
+ zero.xelem (i) = Complex (zeror(i), zeroi(i));
+
// return value
retval(0) = x;
+ retval(1) = zero;
}
-
+
return retval;
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
