AB01ND

Controllable realization for multi-input systems using orthogonal state and input transformations

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

Purpose

  To find a controllable realization for the linear time-invariant
  multi-input system

          dX/dt = A * X + B * U,

  where A and B are N-by-N and N-by-M matrices, respectively,
  which are reduced by this routine to orthogonal canonical form
  using (and optionally accumulating) orthogonal similarity
  transformations.  Specifically, the pair (A, B) is reduced to
  the pair (Ac, Bc),  Ac = Z' * A * Z,  Bc = Z' * B,  given by

          [ Acont     *    ]         [ Bcont ]
     Ac = [                ],   Bc = [       ],
          [   0    Auncont ]         [   0   ]

     and

             [ A11 A12  . . .  A1,p-1 A1p ]         [ B1 ]
             [ A21 A22  . . .  A2,p-1 A2p ]         [ 0  ]
             [  0  A32  . . .  A3,p-1 A3p ]         [ 0  ]
     Acont = [  .   .   . . .    .     .  ],   Bc = [ .  ],
             [  .   .     . .    .     .  ]         [ .  ]
             [  .   .       .    .     .  ]         [ .  ]
             [  0   0   . . .  Ap,p-1 App ]         [ 0  ]

  where the blocks  B1, A21, ..., Ap,p-1  have full row ranks and
  p is the controllability index of the pair.  The size of the
  block  Auncont is equal to the dimension of the uncontrollable
  subspace of the pair (A, B).

Specification
      SUBROUTINE AB01ND( JOBZ, N, M, A, LDA, B, LDB, NCONT, INDCON,
     $                   NBLK, Z, LDZ, TAU, TOL, IWORK, DWORK, LDWORK,
     $                   INFO )
C     .. Scalar Arguments ..
      CHARACTER         JOBZ
      INTEGER           INDCON, INFO, LDA, LDB, LDWORK, LDZ, M, N, NCONT
      DOUBLE PRECISION  TOL
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), DWORK(*), TAU(*), Z(LDZ,*)
      INTEGER           IWORK(*), NBLK(*)

Arguments

Mode Parameters

  JOBZ    CHARACTER*1
          Indicates whether the user wishes to accumulate in a
          matrix Z the orthogonal similarity transformations for
          reducing the system, as follows:
          = 'N':  Do not form Z and do not store the orthogonal
                  transformations;
          = 'F':  Do not form Z, but store the orthogonal
                  transformations in the factored form;
          = 'I':  Z is initialized to the unit matrix and the
                  orthogonal transformation matrix Z is returned.

Input/Output Parameters
  N       (input) INTEGER
          The order of the original state-space representation,
          i.e. the order of the matrix A.  N >= 0.

  M       (input) INTEGER
          The number of system inputs, or of columns of B.  M >= 0.

  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the leading N-by-N part of this array must
          contain the original state dynamics matrix A.
          On exit, the leading NCONT-by-NCONT part contains the
          upper block Hessenberg state dynamics matrix Acont in Ac,
          given by Z' * A * Z, of a controllable realization for
          the original system. The elements below the first block-
          subdiagonal are set to zero.

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

  B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
          On entry, the leading N-by-M part of this array must
          contain the input matrix B.
          On exit, the leading NCONT-by-M part of this array
          contains the transformed input matrix Bcont in Bc, given
          by Z' * B, with all elements but the first block set to
          zero.

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

  NCONT   (output) INTEGER
          The order of the controllable state-space representation.

  INDCON  (output) INTEGER
          The controllability index of the controllable part of the
          system representation.

  NBLK    (output) INTEGER array, dimension (N)
          The leading INDCON elements of this array contain the
          the orders of the diagonal blocks of Acont.

  Z       (output) DOUBLE PRECISION array, dimension (LDZ,N)
          If JOBZ = 'I', then the leading N-by-N part of this
          array contains the matrix of accumulated orthogonal
          similarity transformations which reduces the given system
          to orthogonal canonical form.
          If JOBZ = 'F', the elements below the diagonal, with the
          array TAU, represent the orthogonal transformation matrix
          as a product of elementary reflectors. The transformation
          matrix can then be obtained by calling the LAPACK Library
          routine DORGQR.
          If JOBZ = 'N', the array Z is not referenced and can be
          supplied as a dummy array (i.e. set parameter LDZ = 1 and
          declare this array to be Z(1,1) in the calling program).

  LDZ     INTEGER
          The leading dimension of array Z. If JOBZ = 'I' or
          JOBZ = 'F', LDZ >= MAX(1,N); if JOBZ = 'N', LDZ >= 1.

  TAU     (output) DOUBLE PRECISION array, dimension (N)
          The elements of TAU contain the scalar factors of the
          elementary reflectors used in the reduction of B and A.

Tolerances
  TOL     DOUBLE PRECISION
          The tolerance to be used in rank determination when
          transforming (A, B). If the user sets TOL > 0, then
          the given value of TOL is used as a lower bound for the
          reciprocal condition number (see the description of the
          argument RCOND in the SLICOT routine MB03OD);  a
          (sub)matrix whose estimated condition number is less than
          1/TOL is considered to be of full rank.  If the user sets
          TOL <= 0, then an implicitly computed, default tolerance,
          defined by  TOLDEF = N*N*EPS,  is used instead, where EPS
          is the machine precision (see LAPACK Library routine
          DLAMCH).

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.

  LDWORK  INTEGER
          The length of the array DWORK.
          LDWORK >= MAX(1, N, 3*M).
          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.

Method
  Matrix B is first QR-decomposed and the appropriate orthogonal
  similarity transformation applied to the matrix A. Leaving the
  first rank(B) states unchanged, the remaining lower left block
  of A is then QR-decomposed and the new orthogonal matrix, Q1,
  is also applied to the right of A to complete the similarity
  transformation. By continuing in this manner, a completely
  controllable state-space pair (Acont, Bcont) is found for the
  given (A, B), where Acont is upper block Hessenberg with each
  subdiagonal block of full row rank, and Bcont is zero apart from
  its (independent) first rank(B) rows.
  NOTE that the system controllability indices are easily
  calculated from the dimensions of the blocks of Acont.

References
  [1] Konstantinov, M.M., Petkov, P.Hr. and Christov, N.D.
      Orthogonal Invariants and Canonical Forms for Linear
      Controllable Systems.
      Proc. 8th IFAC World Congress, Kyoto, 1, pp. 49-54, 1981.

  [2] Paige, C.C.
      Properties of numerical algorithms related to computing
      controllablity.
      IEEE Trans. Auto. Contr., AC-26, pp. 130-138, 1981.

  [3] Petkov, P.Hr., Konstantinov, M.M., Gu, D.W. and
      Postlethwaite, I.
      Optimal Pole Assignment Design of Linear Multi-Input Systems.
      Leicester University, Report 99-11, May 1996.

Numerical Aspects
                            3
  The algorithm requires 0(N ) operations and is backward stable.

Further Comments
  If the system matrices A and B are badly scaled, it would be
  useful to scale them with SLICOT routine TB01ID, before calling
  the routine.

Example

Program Text

*     AB01ND EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX, MMAX
      PARAMETER        ( NMAX = 20, MMAX = 20 )
      INTEGER          LDA, LDB, LDZ
      PARAMETER        ( LDA = NMAX, LDB = NMAX, LDZ = NMAX )
      INTEGER          LIWORK
      PARAMETER        ( LIWORK = MMAX )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = MAX( NMAX, 3*MMAX ) )
*     .. Local Scalars ..
      DOUBLE PRECISION TOL
      INTEGER          I, INFO, INDCON, J, M, N, NCONT
      CHARACTER*1      JOBZ
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), DWORK(LDWORK),
     $                 TAU(NMAX), Z(LDZ,NMAX)
      INTEGER          IWORK(LIWORK), NBLK(NMAX)
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         AB01ND, DORGQR
*     .. 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, TOL, JOBZ
      IF ( N.LE.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99990 ) N
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
         IF ( M.LE.0 .OR. M.GT.MMAX ) THEN
            WRITE ( NOUT, FMT = 99989 ) M
         ELSE
            READ ( NIN, FMT = * ) ( ( B(I,J), I = 1,N ), J = 1,M )
*           Find a controllable ssr for the given system.
            CALL AB01ND( JOBZ, N, M, A, LDA, B, LDB, NCONT, INDCON,
     $                   NBLK, Z, LDZ, TAU, TOL, IWORK, DWORK, LDWORK,
     $                   INFO )
*
            IF ( INFO.NE.0 ) THEN
               WRITE ( NOUT, FMT = 99998 ) INFO
            ELSE
               WRITE ( NOUT, FMT = 99997 ) NCONT
               WRITE ( NOUT, FMT = 99996 )
               DO 20 I = 1, NCONT
                  WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,NCONT )
   20          CONTINUE
               WRITE ( NOUT, FMT = 99994 ) ( NBLK(I), I = 1,INDCON )
               WRITE ( NOUT, FMT = 99993 )
               DO 40 I = 1, NCONT
                  WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,M )
   40          CONTINUE
               WRITE ( NOUT, FMT = 99992 ) INDCON
               IF ( LSAME( JOBZ, 'F' ) )
     $            CALL DORGQR( N, N, N, Z, LDZ, TAU, DWORK, LDWORK,
     $                         INFO )
               IF ( LSAME( JOBZ, 'F' ).OR.LSAME( JOBZ, 'I' ) ) THEN
                  WRITE ( NOUT, FMT = 99991 )
                  DO 60 I = 1, N
                     WRITE ( NOUT, FMT = 99995 ) ( Z(I,J), J = 1,N )
   60             CONTINUE
               END IF
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' AB01ND EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from AB01ND = ',I2)
99997 FORMAT (' The order of the controllable state-space representati',
     $       'on = ',I2)
99996 FORMAT (/' The transformed state dynamics matrix of a controllab',
     $       'le realization is ')
99995 FORMAT (20(1X,F8.4))
99994 FORMAT (/' and the dimensions of its diagonal blocks are ',
     $       /20(1X,I2))
99993 FORMAT (/' The transformed input/state matrix B of a controllabl',
     $       'e realization is ')
99992 FORMAT (/' The controllability index of the transformed system r',
     $       'epresentation = ',I2)
99991 FORMAT (/' The similarity transformation matrix Z is ')
99990 FORMAT (/' N is out of range.',/' N = ',I5)
99989 FORMAT (/' M is out of range.',/' M = ',I5)
      END
Program Data
 AB01ND EXAMPLE PROGRAM DATA
   3     2     0.0     I
  -1.0   0.0   0.0
  -2.0  -2.0  -2.0
  -1.0   0.0  -3.0
   1.0   0.0   0.0
   0.0   2.0   1.0
Program Results
 AB01ND EXAMPLE PROGRAM RESULTS

 The order of the controllable state-space representation =  2

 The transformed state dynamics matrix of a controllable realization is 
  -3.0000   2.2361
   0.0000  -1.0000

 and the dimensions of its diagonal blocks are 
  2

 The transformed input/state matrix B of a controllable realization is 
   0.0000  -2.2361
   1.0000   0.0000

 The controllability index of the transformed system representation =  1

 The similarity transformation matrix Z is 
   0.0000   1.0000   0.0000
  -0.8944   0.0000  -0.4472
  -0.4472   0.0000   0.8944

Return to index