## SB02OD

### Solution of continuous- or discrete-time algebraic Riccati equations (generalized Schur vectors method)

[Specification] [Arguments] [Method] [References] [Comments] [Example]

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
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 )

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
```   Van Dooren, P.
A Generalized Eigenvalue Approach for Solving Riccati
Equations.
SIAM J. Sci. Stat. Comp., 2, pp. 121-135, 1981.

 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.

 Sima, V.
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.

```
```  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.0
```
Program Results
``` SB02OD EXAMPLE PROGRAM RESULTS

The solution matrix X is
1.7321   1.0000
1.0000   1.7321
```