Purpose
To solve for X either the continuous-time algebraic Riccati
equation
-1
Q + A'X + XA - (L+XB)R (L+XB)' = 0 (1)
or the discrete-time algebraic Riccati equation
-1
X = A'XA - (L+A'XB)(R + B'XB) (L+A'XB)' + Q (2)
where A, B, Q, R, and L are N-by-N, N-by-M, N-by-N, M-by-M and
N-by-M matrices, respectively, such that Q = C'C, R = D'D and
L = C'D; X is an N-by-N symmetric matrix.
The routine also returns the computed values of the closed-loop
spectrum of the system, i.e., the stable eigenvalues lambda(1),
..., lambda(N) of the corresponding Hamiltonian or symplectic
pencil, in the continuous-time case or discrete-time case,
respectively.
-1
Optionally, matrix G = BR B' may be given instead of B and R.
Other options include the case with Q and/or R given in a
factored form, Q = C'C, R = D'D, and with L a zero matrix.
The routine uses the method of deflating subspaces, based on
reordering the eigenvalues in a generalized Schur matrix pair.
A standard eigenproblem is solved in the continuous-time case
if G is given.
Specification
SUBROUTINE SB02OD( DICO, JOBB, FACT, UPLO, JOBL, SORT, N, M, P, A,
$ LDA, B, LDB, Q, LDQ, R, LDR, L, LDL, RCOND, X,
$ LDX, ALFAR, ALFAI, BETA, S, LDS, T, LDT, U,
$ LDU, TOL, IWORK, DWORK, LDWORK, BWORK, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, FACT, JOBB, JOBL, SORT, UPLO
INTEGER INFO, LDA, LDB, LDL, LDQ, LDR, LDS, LDT, LDU,
$ LDWORK, LDX, M, N, P
DOUBLE PRECISION RCOND, TOL
C .. Array Arguments ..
LOGICAL BWORK(*)
INTEGER IWORK(*)
DOUBLE PRECISION A(LDA,*), ALFAI(*), ALFAR(*), B(LDB,*), BETA(*),
$ DWORK(*), L(LDL,*), Q(LDQ,*), R(LDR,*),
$ S(LDS,*), T(LDT,*), U(LDU,*), X(LDX,*)
Arguments
Mode Parameters
DICO CHARACTER*1
Specifies the type of Riccati equation to be solved as
follows:
= 'C': Equation (1), continuous-time case;
= 'D': Equation (2), discrete-time case.
JOBB CHARACTER*1
Specifies whether or not the matrix G is given, instead
of the matrices B and R, as follows:
= 'B': B and R are given;
= 'G': G is given.
FACT CHARACTER*1
Specifies whether or not the matrices Q and/or R (if
JOBB = 'B') are factored, as follows:
= 'N': Not factored, Q and R are given;
= 'C': C is given, and Q = C'C;
= 'D': D is given, and R = D'D;
= 'B': Both factors C and D are given, Q = C'C, R = D'D.
UPLO CHARACTER*1
If JOBB = 'G', or FACT = 'N', specifies which triangle of
the matrices G and Q (if FACT = 'N'), or Q and R (if
JOBB = 'B'), is stored, as follows:
= 'U': Upper triangle is stored;
= 'L': Lower triangle is stored.
JOBL CHARACTER*1
Specifies whether or not the matrix L is zero, as follows:
= 'Z': L is zero;
= 'N': L is nonzero.
JOBL is not used if JOBB = 'G' and JOBL = 'Z' is assumed.
SLICOT Library routine SB02MT should be called just before
SB02OD, for obtaining the results when JOBB = 'G' and
JOBL = 'N'.
SORT CHARACTER*1
Specifies which eigenvalues should be obtained in the top
of the generalized Schur form, as follows:
= 'S': Stable eigenvalues come first;
= 'U': Unstable eigenvalues come first.
Input/Output Parameters
N (input) INTEGER
The actual state dimension, i.e. the order of the matrices
A, Q, and X, and the number of rows of the matrices B
and L. N >= 0.
M (input) INTEGER
The number of system inputs. If JOBB = 'B', M is the
order of the matrix R, and the number of columns of the
matrix B. M >= 0.
M is not used if JOBB = 'G'.
P (input) INTEGER
The number of system outputs. If FACT = 'C' or 'D' or 'B',
P is the number of rows of the matrices C and/or D.
P >= 0.
Otherwise, P is not used.
A (input) DOUBLE PRECISION array, dimension (LDA,N)
The leading N-by-N part of this array must contain the
state matrix A of the system.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,N).
B (input) DOUBLE PRECISION array, dimension (LDB,*)
If JOBB = 'B', the leading N-by-M part of this array must
contain the input matrix B of the system.
If JOBB = 'G', the leading N-by-N upper triangular part
(if UPLO = 'U') or lower triangular part (if UPLO = 'L')
of this array must contain the upper triangular part or
lower triangular part, respectively, of the matrix
-1
G = BR B'. The stricly lower triangular part (if
UPLO = 'U') or stricly upper triangular part (if
UPLO = 'L') is not referenced.
LDB INTEGER
The leading dimension of array B. LDB >= MAX(1,N).
Q (input) DOUBLE PRECISION array, dimension (LDQ,N)
If FACT = 'N' or 'D', the leading N-by-N upper triangular
part (if UPLO = 'U') or lower triangular part (if UPLO =
'L') of this array must contain the upper triangular part
or lower triangular part, respectively, of the symmetric
state weighting matrix Q. The stricly lower triangular
part (if UPLO = 'U') or stricly upper triangular part (if
UPLO = 'L') is not referenced.
If JOBB = 'B', the triangular part of this array defined
by UPLO is modified internally, but is restored on exit.
If FACT = 'C' or 'B', the leading P-by-N part of this
array must contain the output matrix C of the system.
If JOBB = 'B', this part is modified internally, but is
restored on exit.
LDQ INTEGER
The leading dimension of array Q.
LDQ >= MAX(1,N) if FACT = 'N' or 'D',
LDQ >= MAX(1,P) if FACT = 'C' or 'B'.
R (input) DOUBLE PRECISION array, dimension (LDR,M)
If FACT = 'N' or 'C', the leading M-by-M upper triangular
part (if UPLO = 'U') or lower triangular part (if UPLO =
'L') of this array must contain the upper triangular part
or lower triangular part, respectively, of the symmetric
input weighting matrix R. The stricly lower triangular
part (if UPLO = 'U') or stricly upper triangular part (if
UPLO = 'L') is not referenced.
The triangular part of this array defined by UPLO is
modified internally, but is restored on exit.
If FACT = 'D' or 'B', the leading P-by-M part of this
array must contain the direct transmission matrix D of the
system. This part is modified internally, but is restored
on exit.
If JOBB = 'G', this array is not referenced.
LDR INTEGER
The leading dimension of array R.
LDR >= MAX(1,M) if JOBB = 'B' and FACT = 'N' or 'C';
LDR >= MAX(1,P) if JOBB = 'B' and FACT = 'D' or 'B';
LDR >= 1 if JOBB = 'G'.
L (input) DOUBLE PRECISION array, dimension (LDL,M)
If JOBL = 'N' (and JOBB = 'B'), the leading N-by-M part of
this array must contain the cross weighting matrix L.
This part is modified internally, but is restored on exit.
If JOBL = 'Z' or JOBB = 'G', this array is not referenced.
LDL INTEGER
The leading dimension of array L.
LDL >= MAX(1,N) if JOBL = 'N' and JOBB = 'B';
LDL >= 1 if JOBL = 'Z' or JOBB = 'G'.
RCOND (output) DOUBLE PRECISION
An estimate of the reciprocal of the condition number (in
the 1-norm) of the N-th order system of algebraic
equations from which the solution matrix X is obtained.
X (output) DOUBLE PRECISION array, dimension (LDX,N)
The leading N-by-N part of this array contains the
solution matrix X of the problem.
LDX INTEGER
The leading dimension of array X. LDX >= MAX(1,N).
ALFAR (output) DOUBLE PRECISION array, dimension (2*N)
ALFAI (output) DOUBLE PRECISION array, dimension (2*N)
BETA (output) DOUBLE PRECISION array, dimension (2*N)
The generalized eigenvalues of the 2N-by-2N matrix pair,
ordered as specified by SORT (if INFO = 0). For instance,
if SORT = 'S', the leading N elements of these arrays
contain the closed-loop spectrum of the system matrix
A - BF, where F is the optimal feedback matrix computed
based on the solution matrix X. Specifically,
lambda(k) = [ALFAR(k)+j*ALFAI(k)]/BETA(k) for
k = 1,2,...,N.
If DICO = 'C' and JOBB = 'G', the elements of BETA are
set to 1.
S (output) DOUBLE PRECISION array, dimension (LDS,*)
The leading 2N-by-2N part of this array contains the
ordered real Schur form S of the first matrix in the
reduced matrix pencil associated to the optimal problem,
or of the corresponding Hamiltonian matrix, if DICO = 'C'
and JOBB = 'G'. That is,
(S S )
( 11 12)
S = ( ),
(0 S )
( 22)
where S , S and S are N-by-N matrices.
11 12 22
Array S must have 2*N+M columns if JOBB = 'B', and 2*N
columns, otherwise.
LDS INTEGER
The leading dimension of array S.
LDS >= MAX(1,2*N+M) if JOBB = 'B',
LDS >= MAX(1,2*N) if JOBB = 'G'.
T (output) DOUBLE PRECISION array, dimension (LDT,2*N)
If DICO = 'D' or JOBB = 'B', the leading 2N-by-2N part of
this array contains the ordered upper triangular form T of
the second matrix in the reduced matrix pencil associated
to the optimal problem. That is,
(T T )
( 11 12)
T = ( ),
(0 T )
( 22)
where T , T and T are N-by-N matrices.
11 12 22
If DICO = 'C' and JOBB = 'G' this array is not referenced.
LDT INTEGER
The leading dimension of array T.
LDT >= MAX(1,2*N+M) if JOBB = 'B',
LDT >= MAX(1,2*N) if JOBB = 'G' and DICO = 'D',
LDT >= 1 if JOBB = 'G' and DICO = 'C'.
U (output) DOUBLE PRECISION array, dimension (LDU,2*N)
The leading 2N-by-2N part of this array contains the right
transformation matrix U which reduces the 2N-by-2N matrix
pencil to the ordered generalized real Schur form (S,T),
or the Hamiltonian matrix to the ordered real Schur
form S, if DICO = 'C' and JOBB = 'G'. That is,
(U U )
( 11 12)
U = ( ),
(U U )
( 21 22)
where U , U , U and U are N-by-N matrices.
11 12 21 22
LDU INTEGER
The leading dimension of array U. LDU >= MAX(1,2*N).
Tolerances
TOL DOUBLE PRECISION
The tolerance to be used to test for near singularity of
the original matrix pencil, specifically of the triangular
factor obtained during the reduction process. If the user
sets TOL > 0, then the given value of TOL is used as a
lower bound for the reciprocal condition number of that
matrix; a matrix whose estimated condition number is less
than 1/TOL is considered to be nonsingular. If the user
sets TOL <= 0, then a default tolerance, defined by
TOLDEF = EPS, is used instead, where EPS is the machine
precision (see LAPACK Library routine DLAMCH).
This parameter is not referenced if JOBB = 'G'.
Workspace
IWORK INTEGER array, dimension (LIWORK)
LIWORK >= MAX(1,M,2*N) if JOBB = 'B',
LIWORK >= MAX(1,2*N) if JOBB = 'G'.
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK. If JOBB = 'B' and N > 0, DWORK(2) returns the
reciprocal of the condition number of the M-by-M lower
triangular matrix obtained after compressing the matrix
pencil of order 2N+M to obtain a pencil of order 2N.
If INFO = 0 or INFO = 6, DWORK(3) returns the scaling
factor used internally, which should multiply the
submatrix Y2 to recover X from the first N columns of U
(see METHOD).
LDWORK INTEGER
The length of the array DWORK.
LDWORK >= MAX(3,6*N), if JOBB = 'G',
DICO = 'C';
LDWORK >= MAX(7*(2*N+1)+16,16*N), if JOBB = 'G',
DICO = 'D';
LDWORK >= MAX(7*(2*N+1)+16,16*N,2*N+M,3*M), if JOBB = 'B'.
For optimum performance LDWORK should be larger.
BWORK LOGICAL array, dimension (2*N)
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: if the computed extended matrix pencil is singular,
possibly due to rounding errors;
= 2: if the QZ (or QR) algorithm failed;
= 3: if reordering of the (generalized) eigenvalues
failed;
= 4: if after reordering, roundoff changed values of
some complex eigenvalues so that leading eigenvalues
in the (generalized) Schur form no longer satisfy
the stability condition; this could also be caused
due to scaling;
= 5: if the computed dimension of the solution does not
equal N;
= 6: if a singular matrix was encountered during the
computation of the solution matrix X.
Method
The routine uses a variant of the method of deflating subspaces
proposed by van Dooren [1]. See also [2], [3].
It is assumed that (A,B) is stabilizable and (C,A) is detectable.
Under these assumptions the algebraic Riccati equation is known to
have a unique non-negative definite solution.
The first step in the method of deflating subspaces is to form the
extended Hamiltonian matrices, dimension 2N + M given by
discrete-time continuous-time
|A 0 B| |I 0 0| |A 0 B| |I 0 0|
|Q -I L| - z |0 -A' 0|, |Q A' L| - s |0 -I 0|.
|L' 0 R| |0 -B' 0| |L' B' R| |0 0 0|
Next, these pencils are compressed to a form (see [1])
lambda x A - B .
f f
This generalized eigenvalue problem is then solved using the QZ
algorithm and the stable deflating subspace Ys is determined.
If [Y1'|Y2']' is a basis for Ys, then the required solution is
-1
X = Y2 x Y1 .
A standard eigenvalue problem is solved using the QR algorithm in
the continuous-time case when G is given (DICO = 'C', JOBB = 'G').
References
[1] Van Dooren, P.
A Generalized Eigenvalue Approach for Solving Riccati
Equations.
SIAM J. Sci. Stat. Comp., 2, pp. 121-135, 1981.
[2] Mehrmann, V.
The Autonomous Linear Quadratic Control Problem. Theory and
Numerical Solution.
Lect. Notes in Control and Information Sciences, vol. 163,
Springer-Verlag, Berlin, 1991.
[3] Sima, V.
Algorithms for Linear-Quadratic Optimization.
Pure and Applied Mathematics: A Series of Monographs and
Textbooks, vol. 200, Marcel Dekker, Inc., New York, 1996.
Numerical Aspects
This routine is particularly suited for systems where the matrix R is ill-conditioned. Internal scaling is used.Further Comments
To obtain a stabilizing solution of the algebraic Riccati equations set SORT = 'S'. The routine can also compute the anti-stabilizing solutions of the algebraic Riccati equations, by specifying SORT = 'U'.Example
Program Text
* SB02OD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX, MMAX, PMAX
PARAMETER ( NMAX = 20, MMAX = 20, PMAX = 20 )
INTEGER NMAX2M, NMAX2
PARAMETER ( NMAX2M = 2*NMAX+MMAX, NMAX2 = 2*NMAX )
INTEGER LDA, LDB, LDL, LDQ, LDR, LDS, LDT, LDU, LDX
PARAMETER ( LDA = NMAX, LDB = NMAX, LDL = NMAX,
$ LDQ = MAX(NMAX,PMAX), LDR = MAX(MMAX,PMAX),
$ LDS = NMAX2M, LDT = NMAX2M, LDU = NMAX2,
$ LDX = NMAX )
INTEGER LIWORK
PARAMETER ( LIWORK = MAX(MMAX,NMAX2) )
INTEGER LDWORK
PARAMETER ( LDWORK = MAX(14*NMAX+23,16*NMAX) )
INTEGER LBWORK
PARAMETER ( LBWORK = NMAX2 )
* .. Local Scalars ..
DOUBLE PRECISION RCOND, TOL
INTEGER I, INFO, J, M, N, P
CHARACTER*1 DICO, FACT, JOBB, JOBL, SORT, UPLO
LOGICAL LJOBB
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), ALFAI(NMAX2), ALFAR(NMAX2),
$ B(LDB,MMAX), BETA(NMAX2), DWORK(LDWORK),
$ L(LDL,MMAX), Q(LDQ,NMAX), R(LDR,MMAX),
$ S(LDS,NMAX2M), T(LDT,NMAX2), U(LDU,NMAX2),
$ X(LDX,NMAX)
INTEGER IWORK(LIWORK)
LOGICAL BWORK(LBWORK)
C .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL SB02OD
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, M, P, TOL, DICO, JOBB, FACT, UPLO, JOBL,
$ SORT
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99995 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99994 ) M
ELSE
LJOBB = LSAME( JOBB, 'B' )
IF ( LJOBB ) THEN
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,N )
END IF
IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
WRITE ( NOUT, FMT = 99993 ) P
ELSE
IF ( LSAME( FACT, 'N' ) .OR. LSAME( FACT, 'D' ) ) THEN
READ ( NIN, FMT = * )
$ ( ( Q(I,J), J = 1,N ), I = 1,N )
ELSE
READ ( NIN, FMT = * )
$ ( ( Q(I,J), J = 1,N ), I = 1,P )
END IF
IF ( LJOBB ) THEN
IF ( LSAME( FACT, 'N' ) .OR. LSAME( FACT, 'C' ) ) THEN
READ ( NIN, FMT = * )
$ ( ( R(I,J), J = 1,M ), I = 1,M )
ELSE
READ ( NIN, FMT = * )
$ ( ( R(I,J), J = 1,M ), I = 1,P )
END IF
IF ( LSAME( JOBL, 'N' ) )
$ READ ( NIN, FMT = * )
$ ( ( L(I,J), J = 1,M ), I = 1,N )
END IF
* Find the solution matrix X.
CALL SB02OD( DICO, JOBB, FACT, UPLO, JOBL, SORT, N, M, P,
$ A, LDA, B, LDB, Q, LDQ, R, LDR, L, LDL,
$ RCOND, X, LDX, ALFAR, ALFAI, BETA, S, LDS,
$ T, LDT, U, LDU, TOL, IWORK, DWORK, LDWORK,
$ BWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 )
DO 20 I = 1, N
WRITE ( NOUT, FMT = 99996 ) ( X(I,J), J = 1,N )
20 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' SB02OD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB02OD = ',I2)
99997 FORMAT (' The solution matrix X is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' N is out of range.',/' N = ',I5)
99994 FORMAT (/' M is out of range.',/' M = ',I5)
99993 FORMAT (/' P is out of range.',/' P = ',I5)
END
Program Data
SB02OD EXAMPLE PROGRAM DATA 2 1 3 0.0 C B B U Z S 0.0 1.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0Program Results
SB02OD EXAMPLE PROGRAM RESULTS The solution matrix X is 1.7321 1.0000 1.0000 1.7321