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.0Program 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