### Solution of continuous- or discrete-time algebraic Riccati equations for descriptor systems

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

Purpose

```  To solve for X either the continuous-time algebraic Riccati
equation
-1
Q + A'XE + E'XA - (L+E'XB)R  (L+E'XB)' = 0 ,              (1)

or the discrete-time algebraic Riccati equation
-1
E'XE = A'XA - (L+A'XB)(R + B'XB)  (L+A'XB)' + Q ,         (2)

where A, E, B, Q, R, and L are N-by-N, 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 pencil (A - BF,E), where F is
the optimal gain matrix,
-1
F = R  (L+E'XB)' ,        for (1),

and
-1
F = (R+B'XB)  (L+A'XB)' , for (2).
-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.

It is assumed that E is nonsingular, but this condition is not
checked. Note that the definition (1) of the continuous-time
algebraic Riccati equation, and the formula for the corresponding
optimal gain matrix, require R to be nonsingular, but the
associated linear quadratic optimal problem could have a unique
solution even when matrix R is singular, under mild assumptions
(see METHOD). The routine SG02AD works accordingly in this case.

```
Specification
```      SUBROUTINE SG02AD( DICO, JOBB, FACT, UPLO, JOBL, SCAL, SORT, ACC,
\$                   N, M, P, A, LDA, E, LDE, B, LDB, Q, LDQ, R,
\$                   LDR, L, LDL, RCONDU, X, LDX, ALFAR, ALFAI,
\$                   BETA, S, LDS, T, LDT, U, LDU, TOL, IWORK,
\$                   DWORK, LDWORK, BWORK, IWARN, INFO )
C     .. Scalar Arguments ..
CHARACTER         ACC, DICO, FACT, JOBB, JOBL, SCAL, SORT, UPLO
INTEGER           INFO, IWARN, LDA, LDB, LDE, LDL, LDQ, LDR, LDS,
\$                  LDT, LDU, LDWORK, LDX, M, N, P
DOUBLE PRECISION  RCONDU, TOL
C     .. Array Arguments ..
LOGICAL           BWORK(*)
INTEGER           IWORK(*)
DOUBLE PRECISION  A(LDA,*), ALFAI(*), ALFAR(*), B(LDB,*), BETA(*),
\$                  DWORK(*), E(LDE,*), 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, or Q and R, 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
SG02AD, for obtaining the results when JOBB = 'G' and
JOBL = 'N'.

SCAL    CHARACTER*1
If JOBB = 'B', specifies whether or not a scaling strategy
should be used to scale Q, R, and L, as follows:
= 'G':  General scaling should be used;
= 'N':  No scaling should be used.
SCAL is not used if JOBB = 'G'.

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.

ACC     CHARACTER*1
Specifies whether or not iterative refinement should be
used to solve the system of algebraic equations giving
the solution matrix X, as follows:
= 'R':  Use iterative refinement;
= 'N':  Do not use iterative refinement.

```
Input/Output Parameters
```  N       (input) INTEGER
The actual state dimension, i.e., the order of the
matrices A, E, 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 descriptor system.

LDA     INTEGER
The leading dimension of array A.  LDA >= MAX(1,N).

E       (input) DOUBLE PRECISION array, dimension (LDE,N)
The leading N-by-N part of this array must contain the
matrix E of the descriptor system.

LDE     INTEGER
The leading dimension of array E.  LDE >= 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 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' and SCAL = 'G', then Q 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,*)
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.
If FACT = 'D' or 'B', the leading P-by-M part of this
array must contain the direct transmission matrix D of the
system.
If JOBB = 'B' and SCAL = 'G', then R 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,*)
If JOBL = 'N' and JOBB = 'B', the leading N-by-M part of
this array must contain the cross weighting matrix L.
If JOBB = 'B' and SCAL = 'G', then L 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'.

RCONDU  (output) DOUBLE PRECISION
If N > 0 and INFO = 0 or INFO = 7, 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)
If INFO = 0, 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, or INFO >= 5).
For instance, if SORT = 'S', the leading N elements of
these arrays contain the closed-loop spectrum of the
system. Specifically,
lambda(k) = [ALFAR(k)+j*ALFAI(k)]/BETA(k) for
k = 1,2,...,N.

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,
corresponding to the scaled Q, R, and L, if JOBB = 'B'
and SCAL = '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)
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, corresponding to the scaled Q, R, and L, if
JOBB = 'B' and SCAL = 'G'. That is,

(T   T  )
( 11  12)
T = (       ),
(0   T  )
(     22)

where T  , T   and T   are N-by-N matrices.
11   12      22

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'.

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).
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
If JOBB = 'B' and SCAL = 'G', then U corresponds to the
scaled pencil. If a basis for the stable deflating
subspace of the original problem is needed, then the
submatrix U   must be multiplied by the scaling factor
21
contained in DWORK(4).

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
M-by-M 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 bottom
right lower triangular matrix obtained while compressing
the matrix pencil of order 2N+M to obtain a pencil of
order 2N. If ACC = 'R', and INFO = 0 or INFO = 7, DWORK(3)
returns the reciprocal pivot growth factor (see SLICOT
Library routine MB02PD) for the LU factorization of the
coefficient matrix of the system of algebraic equations
giving the solution matrix X; if DWORK(3) is much
less than 1, then the computed X and RCONDU could be
unreliable. If INFO = 0 or INFO = 7, DWORK(4) returns the
scaling factor used to scale Q, R, and L. DWORK(4) is set
to 1 if JOBB = 'G' or SCAL = 'N'.

LDWORK  INTEGER
The length of the array DWORK.
LDWORK >= MAX(7*(2*N+1)+16,16*N),           if JOBB = 'G';
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)

```
Warning Indicator
```  IWARN   INTEGER
= 0:  no warning;
= 1:  the computed solution may be inaccurate due to poor
scaling or eigenvalues too close to the boundary of
the stability domain (the imaginary axis, if
DICO = 'C', or the unit circle, if DICO = 'D').

```
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 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 the spectrum is too close to the boundary of
the stability domain;
= 7:  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 E is nonsingular, the triple (E,A,B) is
strongly stabilizable and detectable (see ); if, in addition,

-    [ Q   L ]
R := [       ] >= 0 ,
[ L'  R ]

then the pencils

discrete-time                   continuous-time

|A   0   B|     |E   0   0|    |A   0   B|     |E   0   0|
|Q  -E'  L| - z |0  -A'  0| ,  |Q   A'  L| - s |0  -E'  0| ,   (3)
|L'  0   R|     |0  -B'  0|    |L'  B'  R|     |0   0   0|

are dichotomic, i.e., they have no eigenvalues on the boundary of
the stability domain. The above conditions are sufficient for
regularity of these pencils. A necessary condition is that
rank([ B'  L'  R']') = m.

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 matrices in (3), of order 2N + M. Next, these pencils are
compressed to a form of order 2N (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  .

```
References
```   Van Dooren, P.
A Generalized Eigenvalue Approach for Solving Riccati
Equations.
SIAM J. Sci. Stat. Comp., 2, pp. 121-135, 1981.

 Arnold, III, W.F. and Laub, A.J.
Generalized Eigenproblem Algorithms and Software for
Algebraic Riccati Equations.
Proc. IEEE, 72, 1746-1754, 1984.

 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, or even singular.

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

```*     SG02AD 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, NMMAX
PARAMETER        ( NMAX2M = 2*NMAX+MMAX, NMAX2 = 2*NMAX,
\$                   NMMAX  = MAX(NMAX,MMAX) )
INTEGER          LDA, LDB, LDE, LDL, LDQ, LDR, LDS, LDT, LDU, LDX
PARAMETER        ( LDA = NMAX, LDB = NMAX, LDE = 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,2*NMAX+MMAX,
\$                                3*MMAX) )
INTEGER          LBWORK
PARAMETER        ( LBWORK = NMAX2 )
*     .. Local Scalars ..
DOUBLE PRECISION RCONDU, TOL
INTEGER          I, INFO, IWARN, J, M, N, P
CHARACTER*1      ACC, DICO, FACT, JOBB, JOBL, SCAL, SORT, UPLO
LOGICAL          LJOBB
*     .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX),  ALFAI(NMAX2),  ALFAR(NMAX2),
\$                 B(LDB,NMMAX), BETA(NMAX2),   DWORK(LDWORK),
\$                 E(LDE,NMAX),  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)
*     .. External Functions ..
LOGICAL          LSAME
EXTERNAL         LSAME
*     .. External Subroutines ..
*     .. 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,
\$                      SCAL, SORT, ACC
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 )
READ ( NIN, FMT = * ) ( ( E(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 SG02AD( DICO, JOBB, FACT, UPLO, JOBL, SCAL, SORT,
\$                      ACC, N, M, P, A, LDA, E, LDE, B, LDB, Q,
\$                      LDQ, R, LDR, L, LDL, RCONDU, X, LDX, ALFAR,
\$                      ALFAI, BETA, S, LDS, T, LDT, U, LDU, TOL,
\$                      IWORK, DWORK, LDWORK, BWORK, IWARN, 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 (' SG02AD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SG02AD = ',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
``` SG02AD EXAMPLE PROGRAM DATA
2     1     3     0.0     C     B     B     U     Z     N     S     N
0.0  1.0
0.0  0.0
1.0  0.0
0.0  1.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
``` SG02AD EXAMPLE PROGRAM RESULTS

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