SB02ND

Optimal state feedback matrix for an optimal control problem

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

Purpose

  To compute the optimal feedback matrix F for the problem of
  optimal control given by

                     -1
       F = (R + B'XB)  (B'XA + L')                           (1)

  in the discrete-time case and

            -1
       F = R  (B'X + L')                                     (2)

  in the continuous-time case, where A, B and L are N-by-N, N-by-M
  and N-by-M matrices respectively; R and X are M-by-M and N-by-N
  symmetric matrices respectively.

  Optionally, matrix R may be specified in a factored form, and L
  may be zero.

Specification
      SUBROUTINE SB02ND( DICO, FACT, UPLO, JOBL, N, M, P, A, LDA, B,
     $                   LDB, R, LDR, IPIV, L, LDL, X, LDX, RNORM, F,
     $                   LDF, OUFACT, IWORK, DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER         DICO, FACT, JOBL, UPLO
      INTEGER           INFO, LDA, LDB, LDF, LDL, LDR, LDWORK, LDX, M,
     $                  N, P
      DOUBLE PRECISION  RNORM
C     .. Array Arguments ..
      INTEGER           IPIV(*), IWORK(*), OUFACT(2)
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), DWORK(*), F(LDF,*),
     $                  L(LDL,*), R(LDR,*), X(LDX,*)

Arguments

Mode Parameters

  DICO    CHARACTER*1
          Specifies the equation from which F is to be determined,
          as follows:
          = 'D':  Equation (1), discrete-time case;
          = 'C':  Equation (2), continuous-time case.

  FACT    CHARACTER*1
          Specifies how the matrix R is given (factored or not), as
          follows:
          = 'N':  Array R contains the matrix R;
          = 'D':  Array R contains a P-by-M matrix D, where R = D'D;
          = 'C':  Array R contains the Cholesky factor of R;
          = 'U':  Array R contains the symmetric indefinite UdU' or
                  LdL' factorization of R. This option is not
                  available for DICO = 'D'.

  UPLO    CHARACTER*1
          Specifies which triangle of the possibly factored matrix R
          (or R + B'XB, on exit) is or should be 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.

Input/Output Parameters
  N       (input) INTEGER
          The order of the matrices A and X.  N >= 0.

  M       (input) INTEGER
          The number of system inputs.  M >= 0.

  P       (input) INTEGER
          The number of system outputs.  P >= 0.
          This parameter must be specified only for FACT = 'D'.

  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
          If DICO = 'D', the leading N-by-N part of this array must
          contain the state matrix A of the system.
          If DICO = 'C', this array is not referenced.

  LDA     INTEGER
          The leading dimension of array A.
          LDA >= MAX(1,N) if DICO = 'D';
          LDA >= 1        if DICO = 'C'.

  B       (input) DOUBLE PRECISION array, dimension (LDB,M)
          The leading N-by-M part of this array must contain the
          input matrix B of the system.
          If DICO = 'D' and FACT = 'D' or 'C', the contents of this
          array is destroyed.
          Otherwise, B is unchanged on exit.

  LDB     INTEGER
          The leading dimension of array B.  LDB >= MAX(1,N).

  R       (input/output) DOUBLE PRECISION array, dimension (LDR,M)
          On entry, if FACT = 'N', 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.
          On entry, if FACT = 'D', the leading P-by-M part of this
          array must contain the direct transmission matrix D of the
          system.
          On entry, if FACT = '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 Cholesky
          factor of the positive definite input weighting matrix R
          (as produced by LAPACK routine DPOTRF).
          On entry, if DICO = 'C' and FACT = 'U', the leading M-by-M
          upper triangular part (if UPLO = 'U') or lower triangular
          part (if UPLO = 'L') of this array must contain the
          factors of the UdU' or LdL' factorization, respectively,
          of the symmetric indefinite input weighting matrix R (as
          produced by LAPACK routine DSYTRF).
          The stricly lower triangular part (if UPLO = 'U') or
          stricly upper triangular part (if UPLO = 'L') of this
          array is used as workspace.
          On exit, if OUFACT(1) = 1, and INFO = 0 (or INFO = M+1),
          the leading M-by-M upper triangular part (if UPLO = 'U')
          or lower triangular part (if UPLO = 'L') of this array
          contains the Cholesky factor of the given input weighting
          matrix (for DICO = 'C'), or that of the matrix R + B'XB
          (for DICO = 'D').
          On exit, if OUFACT(1) = 2, and INFO = 0 (or INFO = M+1),
          the leading M-by-M upper triangular part (if UPLO = 'U')
          or lower triangular part (if UPLO = 'L') of this array
          contains the factors of the UdU' or LdL' factorization,
          respectively, of the given input weighting matrix
          (for DICO = 'C'), or that of the matrix R + B'XB
          (for DICO = 'D').
          On exit R is unchanged if FACT = 'U'.

  LDR     INTEGER.
          The leading dimension of the array R.
          LDR >= MAX(1,M)   if FACT <> 'D';
          LDR >= MAX(1,M,P) if FACT =  'D'.

  IPIV    (input/output) INTEGER array, dimension (M)
          On entry, if FACT = 'U', this array must contain details
          of the interchanges performed and the block structure of
          the d factor in the UdU' or LdL' factorization of matrix R
          (as produced by LAPACK routine DSYTRF).
          On exit, if OUFACT(1) = 2, this array contains details of
          the interchanges performed and the block structure of the
          d factor in the UdU' or LdL' factorization of matrix R (or
          D'D) or R + B'XB (or D'D + B'XB), as produced by LAPACK
          routine DSYTRF.
          This array is not referenced for DICO = 'D' or FACT = 'D',
          or 'C'.

  L       (input) DOUBLE PRECISION array, dimension (LDL,M)
          If JOBL = 'N', the leading N-by-M part of this array must
          contain the cross weighting matrix L.
          If JOBL = 'Z', this array is not referenced.

  LDL     INTEGER
          The leading dimension of array L.
          LDL >= MAX(1,N) if JOBL = 'N';
          LDL >= 1        if JOBL = 'Z'.

  X       (input/output) DOUBLE PRECISION array, dimension (LDX,N)
          On entry, the leading N-by-N part of this array must
          contain the solution matrix X of the algebraic Riccati
          equation as produced by SLICOT Library routines SB02MD or
          SB02OD. Matrix X is assumed non-negative definite.
          On exit, if DICO = 'D', FACT = 'D' or 'C', OUFACT(2) = 1,
          and INFO = 0, the N-by-N upper triangular part of this
          array contains the Cholesky factor of the given matrix X,
          which is found to be positive definite.
          On exit, if DICO = 'D', FACT = 'D' or 'C', OUFACT(2) = 2,
          and INFO = 0, the leading N-by-N part of this array
          contains the matrix of orthonormal eigenvectors of X.
          On exit X is unchanged if DICO = 'C' or FACT = 'N'.

  LDX     INTEGER
          The leading dimension of array X.  LDX >= MAX(1,N).

  RNORM   (input) DOUBLE PRECISION
          If FACT = 'U', this parameter must contain the 1-norm of
          the original matrix R (before factoring it).
          Otherwise, this parameter is not used.

  F       (output) DOUBLE PRECISION array, dimension (LDF,N)
          The leading M-by-N part of this array contains the
          optimal feedback matrix F.

  LDF     INTEGER
          The leading dimension of array F.  LDF >= MAX(1,M).

  OUFACT  (output) INTEGER array, dimension (2)
          Information about the factorization finally used.
          OUFACT(1) = 1:  Cholesky factorization of R (or R + B'XB)
                          has been used;
          OUFACT(1) = 2:  UdU' (if UPLO = 'U') or LdL' (if UPLO =
                          'L') factorization of R (or R + B'XB)
                          has been used;
          OUFACT(2) = 1:  Cholesky factorization of X has been used;
          OUFACT(2) = 2:  Spectral factorization of X has been used.
          The value of OUFACT(2) is not set for DICO = 'C' or for
          DICO = 'D' and FACT = 'N'.

Workspace
  IWORK   INTEGER array, dimension (M)

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if INFO = 0, DWORK(1) returns the optimal value
          of LDWORK, and DWORK(2) contains the reciprocal condition
          number of the matrix R (for DICO = 'C') or of R + B'XB
          (for DICO = 'D').
          If on exit INFO = 0, and OUFACT(2) = 2, then DWORK(3),...,
          DWORK(N+2) contain the eigenvalues of X, in ascending
          order.

  LDWORK  INTEGER
          Dimension of working array DWORK.
          LDWORK >= max(2,3*M)         if FACT = 'N';
          LDWORK >= max(2,2*M)         if FACT = 'U';
          LDWORK >= max(2,3*M)         if FACT = 'C', DICO = 'C';
          LDWORK >= N+3*M+2            if FACT = 'C', DICO = 'D';
          LDWORK >= max(2,min(P,M)+M)  if FACT = 'D', DICO = 'C';
          LDWORK >= max(N+3*M+2,4*N+1) if FACT = 'D', DICO = 'D'.
          For optimum performance LDWORK should be larger.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = i:  if the i-th element of the d factor is exactly zero;
                the UdU' (or LdL') factorization has been completed,
                but the block diagonal matrix d is exactly singular;
          = M+1:  if the matrix R (if DICO = 'C'), or R + B'XB
                (if DICO = 'D') is numerically singular (to working
                precision);
          = M+2:  if one or more of the eigenvalues of X has not
                converged.

Method
  The optimal feedback matrix F is obtained as the solution to the
  system of linear equations

     (R + B'XB) * F = B'XA + L'

  in the discrete-time case and

     R * F = B'X + L'

  in the continuous-time case, with R replaced by D'D if FACT = 'D'.
  The factored form of R, specified by FACT <> 'N', is taken into
  account. If FACT = 'N', Cholesky factorization is tried first, but
  if the coefficient matrix is not positive definite, then UdU' (or
  LdL') factorization is used. The discrete-time case involves
  updating of a triangular factorization of R (or D'D); Cholesky or
  symmetric spectral factorization of X is employed to avoid
  squaring of the condition number of the matrix. When D is given,
  its QR factorization is determined, and the triangular factor is
  used as described above.

Numerical Aspects
  The algorithm consists of numerically stable steps.
                                 3     2
  For DICO = 'C', it requires O(m  + mn ) floating point operations
                        2
  if FACT = 'N' and O(mn ) floating point operations, otherwise.
  For DICO = 'D', the operation counts are similar, but additional
     3
  O(n ) floating point operations may be needed in the worst case.

Further Comments
  None
Example

Program Text

*     SB02ND EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2010 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX, MMAX, PMAX
      PARAMETER        ( NMAX = 20, MMAX = 20, PMAX = 20 )
      INTEGER          NMAX2
      PARAMETER        ( NMAX2 = 2*NMAX )
      INTEGER          LDA, LDB, LDC, LDL, LDR, LDS, LDT, LDU, LDX, LDF
      PARAMETER        ( LDA = NMAX, LDB = NMAX, LDC = PMAX, LDL = NMAX,
     $                   LDR = MAX(MMAX,PMAX), LDS = NMAX2+MMAX,
     $                   LDT = NMAX2+MMAX, LDU = NMAX2, LDX = NMAX,
     $                   LDF = MMAX )
      INTEGER          LIWORK
      PARAMETER        ( LIWORK = MAX( NMAX2,MMAX ) )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = MAX( NMAX+3*MMAX+2, 14*NMAX+23,
     $                   16*NMAX ) )
*     .. Local Scalars ..
      DOUBLE PRECISION TOL, RCOND, RNORM
      INTEGER          I, INFO1, INFO2, J, M, N, P
      CHARACTER*1      DICO, FACT, JOBB, JOBL, SORT, UPLO
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), ALFAI(2*NMAX), ALFAR(2*NMAX),
     $                 B(LDB,MMAX), BETA(2*NMAX), C(LDC,NMAX),
     $                 DWORK(LDWORK), F(LDF,NMAX), L(LDL,MMAX),
     $                 R(LDR,MMAX), S(LDS,NMAX2+MMAX), T(LDT,NMAX2),
     $                 U(LDU,NMAX2), X(LDX,NMAX)
      INTEGER          IPIV(LIWORK), IWORK(LIWORK), OUFACT(2)
      LOGICAL          BWORK(NMAX2)
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         SB02ND, 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, FACT, JOBL, UPLO
      IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99993 ) 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 = 99992 ) M
         ELSE
            READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
            IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
               WRITE ( NOUT, FMT = 99991 ) P
            ELSE
               READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
               IF ( LSAME( FACT, 'D' ) ) THEN
                  READ ( NIN, FMT = * ) ( ( R(I,J), J = 1,M ), I = 1,P )
               ELSE
                  READ ( NIN, FMT = * ) ( ( R(I,J), J = 1,M ), I = 1,M )
               END IF
*              Find the solution matrix X.
               JOBB = 'B'
               SORT = 'S'
               CALL SB02OD( DICO, JOBB, 'Both', UPLO, JOBL, SORT, N, M,
     $                      P, A, LDA, B, LDB, C, LDC, R, LDR, L, LDL,
     $                      RCOND, X, LDX, ALFAR, ALFAI, BETA, S, LDS,
     $                      T, LDT, U, LDU, TOL, IWORK, DWORK, LDWORK,
     $                      BWORK, INFO1 )
*
               IF ( INFO1.NE.0 ) THEN
                  WRITE ( NOUT, FMT = 99998 ) INFO1
               ELSE
                  WRITE ( NOUT, FMT = 99996 )
                  DO 20 I = 1, N
                     WRITE ( NOUT, FMT = 99994 ) ( X(I,J), J = 1,N )
   20             CONTINUE
*                 Compute the optimal feedback matrix F.
                  CALL SB02ND( DICO, FACT, UPLO, JOBL, N, M, P, A, LDA,
     $                         B, LDB, R, LDR, IPIV, L, LDL, X, LDX,
     $                         RNORM, F, LDF, OUFACT, IWORK, DWORK,
     $                         LDWORK, INFO2 )
*
                  IF ( INFO2.NE.0 ) THEN
                     WRITE ( NOUT, FMT = 99997 ) INFO2
                  ELSE
                     WRITE ( NOUT, FMT = 99995 )
                     DO 40 I = 1, M
                        WRITE ( NOUT, FMT = 99994 ) ( F(I,J), J = 1,N )
   40                CONTINUE
                  END IF
               END IF
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' SB02ND EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB02OD = ',I2)
99997 FORMAT (' INFO on exit from SB02ND = ',I2)
99996 FORMAT (' The solution matrix X is ')
99995 FORMAT (/' The optimal feedback matrix F is ')
99994 FORMAT (20(1X,F8.4))
99993 FORMAT (/' N is out of range.',/' N = ',I5)
99992 FORMAT (/' M is out of range.',/' M = ',I5)
99991 FORMAT (/' P is out of range.',/' P = ',I5)
      END
Program Data
 SB02ND EXAMPLE PROGRAM DATA
   2     1     3     0.0     D     N     Z     U
   2.0 -1.0
   1.0  0.0
   1.0
   0.0
   0.0  0.0
   0.0  0.0
   0.0  1.0
   0.0
   0.0
   0.0
Program Results
 SB02ND EXAMPLE PROGRAM RESULTS

 The solution matrix X is 
   1.0000   0.0000
   0.0000   1.0000

 The optimal feedback matrix F is 
   2.0000  -1.0000

Click here to get a compressed (gzip) tar file containing the source code of the routine, the example program, data, documentation, and related files.

Return to index