Purpose
To find the transfer matrix T(s) of a given state-space representation (A,B,C,D). T(s) is expressed as either row or column polynomial vectors over monic least common denominator polynomials.Specification
SUBROUTINE TB04AD( ROWCOL, N, M, P, A, LDA, B, LDB, C, LDC, D,
$ LDD, NR, INDEX, DCOEFF, LDDCOE, UCOEFF, LDUCO1,
$ LDUCO2, TOL1, TOL2, IWORK, DWORK, LDWORK,
$ INFO )
C .. Scalar Arguments ..
CHARACTER ROWCOL
INTEGER INFO, LDA, LDB, LDC, LDD, LDDCOE, LDUCO1,
$ LDUCO2, LDWORK, M, N, NR, P
DOUBLE PRECISION TOL1, TOL2
C .. Array Arguments ..
INTEGER INDEX(*), IWORK(*)
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
$ DCOEFF(LDDCOE,*), DWORK(*),
$ UCOEFF(LDUCO1,LDUCO2,*)
Arguments
Mode Parameters
ROWCOL CHARACTER*1
Indicates whether the transfer matrix T(s) is required
as rows or columns over common denominators as follows:
= 'R': T(s) is required as rows over common denominators;
= 'C': T(s) is required as columns over common
denominators.
Input/Output Parameters
N (input) INTEGER
The order of the state-space representation, i.e. the
order of the original state dynamics matrix A. N >= 0.
M (input) INTEGER
The number of system inputs. M >= 0.
P (input) INTEGER
The number of system outputs. P >= 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 NR-by-NR part of this array contains
the upper block Hessenberg state dynamics matrix A of a
transformed representation for the original system: this
is completely controllable if ROWCOL = 'R', or completely
observable if ROWCOL = 'C'.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,N).
B (input/output) DOUBLE PRECISION array, dimension (LDB,M),
if ROWCOL = 'R', and (LDB,MAX(M,P)) if ROWCOL = 'C'.
On entry, the leading N-by-M part of this array must
contain the original input/state matrix B; if
ROWCOL = 'C', the remainder of the leading N-by-MAX(M,P)
part is used as internal workspace.
On exit, the leading NR-by-M part of this array contains
the transformed input/state matrix B.
LDB INTEGER
The leading dimension of array B. LDB >= MAX(1,N).
C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
On entry, the leading P-by-N part of this array must
contain the original state/output matrix C; if
ROWCOL = 'C', the remainder of the leading MAX(M,P)-by-N
part is used as internal workspace.
On exit, the leading P-by-NR part of this array contains
the transformed state/output matrix C.
LDC INTEGER
The leading dimension of array C.
LDC >= MAX(1,P) if ROWCOL = 'R';
LDC >= MAX(1,M,P) if ROWCOL = 'C'.
D (input) DOUBLE PRECISION array, dimension (LDD,M),
if ROWCOL = 'R', and (LDD,MAX(M,P)) if ROWCOL = 'C'.
The leading P-by-M part of this array must contain the
original direct transmission matrix D; if ROWCOL = 'C',
this array is modified internally, but restored on exit,
and the remainder of the leading MAX(M,P)-by-MAX(M,P)
part is used as internal workspace.
LDD INTEGER
The leading dimension of array D.
LDD >= MAX(1,P) if ROWCOL = 'R';
LDD >= MAX(1,M,P) if ROWCOL = 'C'.
NR (output) INTEGER
The order of the transformed state-space representation.
INDEX (output) INTEGER array, dimension (porm), where porm = P,
if ROWCOL = 'R', and porm = M, if ROWCOL = 'C'.
The degrees of the denominator polynomials.
DCOEFF (output) DOUBLE PRECISION array, dimension (LDDCOE,N+1)
The leading porm-by-kdcoef part of this array contains
the coefficients of each denominator polynomial, where
kdcoef = MAX(INDEX(I)) + 1.
DCOEFF(I,K) is the coefficient in s**(INDEX(I)-K+1) of
the I-th denominator polynomial, where K = 1,2,...,kdcoef.
LDDCOE INTEGER
The leading dimension of array DCOEFF.
LDDCOE >= MAX(1,P) if ROWCOL = 'R';
LDDCOE >= MAX(1,M) if ROWCOL = 'C'.
UCOEFF (output) DOUBLE PRECISION array, dimension
(LDUCO1,LDUCO2,N+1)
If ROWCOL = 'R' then porp = M, otherwise porp = P.
The leading porm-by-porp-by-kdcoef part of this array
contains the coefficients of the numerator matrix U(s).
UCOEFF(I,J,K) is the coefficient in s**(INDEX(iorj)-K+1)
of polynomial (I,J) of U(s), where K = 1,2,...,kdcoef;
if ROWCOL = 'R' then iorj = I, otherwise iorj = J.
Thus for ROWCOL = 'R', U(s) =
diag(s**INDEX(I))*(UCOEFF(.,.,1)+UCOEFF(.,.,2)/s+...).
LDUCO1 INTEGER
The leading dimension of array UCOEFF.
LDUCO1 >= MAX(1,P) if ROWCOL = 'R';
LDUCO1 >= MAX(1,M) if ROWCOL = 'C'.
LDUCO2 INTEGER
The second dimension of array UCOEFF.
LDUCO2 >= MAX(1,M) if ROWCOL = 'R';
LDUCO2 >= MAX(1,P) if ROWCOL = 'C'.
Tolerances
TOL1 DOUBLE PRECISION
The tolerance to be used in determining the i-th row of
T(s), where i = 1,2,...,porm. If the user sets TOL1 > 0,
then the given value of TOL1 is used as an absolute
tolerance; elements with absolute value less than TOL1 are
considered neglijible. If the user sets TOL1 <= 0, then
an implicitly computed, default tolerance, defined in
the SLICOT Library routine TB01ZD, is used instead.
TOL2 DOUBLE PRECISION
The tolerance to be used to separate out a controllable
subsystem of (A,B,C). If the user sets TOL2 > 0, then
the given value of TOL2 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/TOL2 is considered to be of full rank. If the user sets
TOL2 <= 0, then an implicitly computed, default tolerance,
defined in the SLICOT Library routine TB01UD, is used
instead.
Workspace
IWORK INTEGER array, dimension (N+MAX(M,P))
On exit, if INFO = 0, the first nonzero elements of
IWORK(1:N) return the orders of the diagonal blocks of A.
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*(N + 1) + MAX(N*MP + 2*N + MAX(N,MP),
3*MP, PM)),
where MP = M, PM = P, if ROWCOL = 'R';
MP = P, PM = M, if ROWCOL = 'C'.
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
The method for transfer matrices factorized by rows will be described here: T(s) factorized by columns is dealt with by operating on the dual of the original system. Each row of T(s) is simply a single-output relatively left prime polynomial matrix representation, so can be calculated by applying a simplified version of the Orthogonal Structure Theorem to a minimal state-space representation for the corresponding row of the given system. A minimal state-space representation is obtained using the Orthogonal Canonical Form to first separate out a completely controllable one for the overall system and then, for each row in turn, applying it again to the resulting dual SIMO (single-input multi-output) system. Note that the elements of the transformed matrix A so calculated are individually scaled in a way which guarantees a monic denominator polynomial.References
[1] Williams, T.W.C.
An Orthogonal Structure Theorem for Linear Systems.
Control Systems Research Group, Kingston Polytechnic,
Internal Report 82/2, 1982.
Numerical Aspects
3 The algorithm requires 0(N ) operations.Further Comments
NoneExample
Program Text
* TB04AD 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 MAXMP
PARAMETER ( MAXMP = MAX( MMAX, PMAX ) )
INTEGER LDA, LDB, LDC, LDD, LDDCOE, LDUCO1, LDUCO2,
$ NMAXP1
PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = MAXMP,
$ LDD = MAXMP, LDDCOE = MAXMP, LDUCO1 = MAXMP,
$ LDUCO2 = MAXMP, NMAXP1 = NMAX+1 )
INTEGER LIWORK
PARAMETER ( LIWORK = NMAX + MAXMP )
INTEGER LDWORK
PARAMETER ( LDWORK = NMAX*( NMAX + 1 ) +
$ MAX( NMAX*MAXMP + 2*NMAX +
$ MAX( NMAX, MAXMP ), 3*MAXMP ) )
* .. Local Scalars ..
DOUBLE PRECISION TOL1, TOL2
INTEGER I, II, IJ, INDBLK, INFO, J, JJ, KDCOEF, M, N,
$ NR, P, PORM, PORP
CHARACTER*1 ROWCOL
CHARACTER*132 ULINE
LOGICAL LROWCO
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MAXMP), C(LDC,NMAX),
$ D(LDD,MAXMP), DCOEFF(LDDCOE,NMAXP1),
$ DWORK(LDWORK), UCOEFF(LDUCO1,LDUCO2,NMAXP1)
INTEGER INDEX(MAXMP), IWORK(LIWORK)
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL TB04AD
* .. 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, TOL1, TOL2, ROWCOL
LROWCO = LSAME( ROWCOL, 'R' )
ULINE(1:20) = ' '
DO 20 I = 21, 132
ULINE(I:I) = '-'
20 CONTINUE
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99986 ) 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 = 99985 ) M
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), I = 1,N ), J = 1,M )
IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
WRITE ( NOUT, FMT = 99984 ) P
ELSE
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
READ ( NIN, FMT = * ) ( ( D(I,J), J = 1,M ), I = 1,P )
* Find the transfer matrix T(s) of (A,B,C,D).
CALL TB04AD( ROWCOL, N, M, P, A, LDA, B, LDB, C, LDC, D,
$ LDD, NR, INDEX, DCOEFF, LDDCOE, UCOEFF,
$ LDUCO1, LDUCO2, TOL1, TOL2, IWORK, DWORK,
$ LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 ) NR
DO 40 I = 1, NR
WRITE ( NOUT, FMT = 99996 ) ( A(I,J), J = 1,NR )
40 CONTINUE
WRITE ( NOUT, FMT = 99995 )
DO 60 I = 1, NR
WRITE ( NOUT, FMT = 99996 ) ( B(I,J), J = 1,M )
60 CONTINUE
WRITE ( NOUT, FMT = 99994 )
DO 80 I = 1, P
WRITE ( NOUT, FMT = 99996 ) ( C(I,J), J = 1,NR )
80 CONTINUE
INDBLK = 0
DO 100 I = 1, N
IF ( IWORK(I).NE.0 ) INDBLK = INDBLK + 1
100 CONTINUE
IF ( LROWCO ) THEN
PORM = P
PORP = M
WRITE ( NOUT, FMT = 99993 ) INDBLK,
$ ( IWORK(I), I = 1,INDBLK )
ELSE
PORM = M
PORP = P
WRITE ( NOUT, FMT = 99992 ) INDBLK,
$ ( IWORK(I), I = 1,INDBLK )
END IF
WRITE ( NOUT, FMT = 99991 ) ( INDEX(I), I = 1,PORM )
WRITE ( NOUT, FMT = 99990 )
KDCOEF = 0
DO 120 I = 1, PORM
KDCOEF = MAX( KDCOEF, INDEX(I) )
120 CONTINUE
KDCOEF = KDCOEF + 1
DO 160 II = 1, PORM
DO 140 JJ = 1, PORP
WRITE ( NOUT, FMT = 99989 ) II, JJ,
$ ( UCOEFF(II,JJ,IJ), IJ = 1,KDCOEF )
WRITE ( NOUT, FMT = 99988 ) ULINE(1:7*KDCOEF+21)
WRITE ( NOUT, FMT = 99987 )
$ ( DCOEFF(II,IJ), IJ = 1,KDCOEF )
140 CONTINUE
160 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TB04AD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TB04AD = ',I2)
99997 FORMAT (' The order of the transformed state-space representatio',
$ 'n = ',I2,//' The transformed state dynamics matrix A is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' The transformed input/state matrix B is ')
99994 FORMAT (/' The transformed state/output matrix C is ')
99993 FORMAT (/' The controllability index of the transformed state-sp',
$ 'ace representation = ',I2,//' The dimensions of the diag',
$ 'onal blocks of the transformed A are ',/20(I5))
99992 FORMAT (/' The observability index of the transformed state-spac',
$ 'e representation = ',I2,//' The dimensions of the diagon',
$ 'al blocks of the transformed A are ',/20(I5))
99991 FORMAT (/' The degrees of the denominator polynomials are',/20(I5)
$ )
99990 FORMAT (/' The coefficients of polynomials in the transfer matri',
$ 'x T(s) are ')
99989 FORMAT (/' element (',I2,',',I2,') is ',20(1X,F6.2))
99988 FORMAT (1X,A)
99987 FORMAT (20X,20(1X,F6.2))
99986 FORMAT (/' N is out of range.',/' N = ',I5)
99985 FORMAT (/' M is out of range.',/' M = ',I5)
99984 FORMAT (/' P is out of range.',/' P = ',I5)
END
Program Data
TB04AD EXAMPLE PROGRAM DATA 3 2 2 0.0 0.0 R -1.0 0.0 0.0 0.0 -2.0 0.0 0.0 0.0 -3.0 0.0 1.0 -1.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 1.0Program Results
TB04AD EXAMPLE PROGRAM RESULTS
The order of the transformed state-space representation = 3
The transformed state dynamics matrix A is
-2.5000 -0.2887 -0.4082
-0.2887 -1.5000 -0.7071
-0.4082 -0.7071 -2.0000
The transformed input/state matrix B is
-1.4142 -0.7071
0.0000 1.2247
0.0000 0.0000
The transformed state/output matrix C is
0.0000 0.8165 1.1547
0.0000 1.6330 0.5774
The controllability index of the transformed state-space representation = 2
The dimensions of the diagonal blocks of the transformed A are
2 1
The degrees of the denominator polynomials are
2 3
The coefficients of polynomials in the transfer matrix T(s) are
element ( 1, 1) is 1.00 5.00 7.00 0.00
-----------------------------
1.00 5.00 6.00 0.00
element ( 1, 2) is 0.00 1.00 3.00 0.00
-----------------------------
1.00 5.00 6.00 0.00
element ( 2, 1) is 0.00 0.00 1.00 1.00
-----------------------------
1.00 6.00 11.00 6.00
element ( 2, 2) is 1.00 8.00 20.00 15.00
-----------------------------
1.00 6.00 11.00 6.00