Purpose
To compute a reduced order model (Ar,Br,Cr,Dr) for a stable original state-space representation (A,B,C,D) by using the optimal Hankel-norm approximation method in conjunction with square-root balancing.Specification
SUBROUTINE AB09CD( DICO, EQUIL, ORDSEL, N, M, P, NR, A, LDA, B,
$ LDB, C, LDC, D, LDD, HSV, TOL1, TOL2, IWORK,
$ DWORK, LDWORK, IWARN, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, EQUIL, ORDSEL
INTEGER INFO, IWARN, LDA, LDB, LDC, LDD, LDWORK,
$ M, N, NR, P
DOUBLE PRECISION TOL1, TOL2
C .. Array Arguments ..
INTEGER IWORK(*)
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
$ DWORK(*), HSV(*)
Arguments
Mode Parameters
DICO CHARACTER*1
Specifies the type of the original system as follows:
= 'C': continuous-time system;
= 'D': discrete-time system.
EQUIL CHARACTER*1
Specifies whether the user wishes to preliminarily
equilibrate the triplet (A,B,C) as follows:
= 'S': perform equilibration (scaling);
= 'N': do not perform equilibration.
ORDSEL CHARACTER*1
Specifies the order selection method as follows:
= 'F': the resulting order NR is fixed;
= 'A': the resulting order NR is automatically determined
on basis of the given tolerance TOL1.
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. M >= 0.
P (input) INTEGER
The number of system outputs. P >= 0.
NR (input/output) INTEGER
On entry with ORDSEL = 'F', NR is the desired order of
the resulting reduced order system. 0 <= NR <= N.
On exit, if INFO = 0, NR is the order of the resulting
reduced order model. NR is set as follows:
if ORDSEL = 'F', NR is equal to MIN(MAX(0,NR-KR+1),NMIN),
where KR is the multiplicity of the Hankel singular value
HSV(NR+1), NR is the desired order on entry, and NMIN is
the order of a minimal realization of the given system;
NMIN is determined as the number of Hankel singular values
greater than N*EPS*HNORM(A,B,C), where EPS is the machine
precision (see LAPACK Library Routine DLAMCH) and
HNORM(A,B,C) is the Hankel norm of the system (computed
in HSV(1));
if ORDSEL = 'A', NR is equal to the number of Hankel
singular values greater than MAX(TOL1,N*EPS*HNORM(A,B,C)).
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the state dynamics matrix A.
On exit, if INFO = 0, the leading NR-by-NR part of this
array contains the state dynamics matrix Ar of the
reduced order system in a real Schur form.
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 original input/state matrix B.
On exit, if INFO = 0, the leading NR-by-M part of this
array contains the input/state matrix Br of the reduced
order system.
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.
On exit, if INFO = 0, the leading P-by-NR part of this
array contains the state/output matrix Cr of the reduced
order system.
LDC INTEGER
The leading dimension of array C. LDC >= MAX(1,P).
D (input/output) DOUBLE PRECISION array, dimension (LDD,M)
On entry, the leading P-by-M part of this array must
contain the original input/output matrix D.
On exit, if INFO = 0, the leading P-by-M part of this
array contains the input/output matrix Dr of the reduced
order system.
LDD INTEGER
The leading dimension of array D. LDD >= MAX(1,P).
HSV (output) DOUBLE PRECISION array, dimension (N)
If INFO = 0, it contains the Hankel singular values of
the original system ordered decreasingly. HSV(1) is the
Hankel norm of the system.
Tolerances
TOL1 DOUBLE PRECISION
If ORDSEL = 'A', TOL1 contains the tolerance for
determining the order of reduced system.
For model reduction, the recommended value is
TOL1 = c*HNORM(A,B,C), where c is a constant in the
interval [0.00001,0.001], and HNORM(A,B,C) is the
Hankel-norm of the given system (computed in HSV(1)).
For computing a minimal realization, the recommended
value is TOL1 = N*EPS*HNORM(A,B,C), where EPS is the
machine precision (see LAPACK Library Routine DLAMCH).
This value is used by default if TOL1 <= 0 on entry.
If ORDSEL = 'F', the value of TOL1 is ignored.
TOL2 DOUBLE PRECISION
The tolerance for determining the order of a minimal
realization of the given system. The recommended value is
TOL2 = N*EPS*HNORM(A,B,C). This value is used by default
if TOL2 <= 0 on entry.
If TOL2 > 0, then TOL2 <= TOL1.
Workspace
IWORK INTEGER array, dimension (LIWORK)
LIWORK = MAX(1,M), if DICO = 'C';
LIWORK = MAX(1,N,M), if DICO = 'D'.
On exit, if INFO = 0, IWORK(1) contains NMIN, the order of
the computed minimal realization.
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( LDW1, LDW2 ), where
LDW1 = N*(2*N+MAX(N,M,P)+5) + N*(N+1)/2,
LDW2 = N*(M+P+2) + 2*M*P + MIN(N,M) +
MAX( 3*M+1, MIN(N,M)+P ).
For optimum performance LDWORK should be larger.
Warning Indicator
IWARN INTEGER
= 0: no warning;
= 1: with ORDSEL = 'F', the selected order NR is greater
than the order of a minimal realization of the
given system. In this case, the resulting NR is set
automatically to a value corresponding to the order
of a minimal realization of the system.
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: the reduction of A to the real Schur form failed;
= 2: the state matrix A is not stable (if DICO = 'C')
or not convergent (if DICO = 'D');
= 3: the computation of Hankel singular values failed;
= 4: the computation of stable projection failed;
= 5: the order of computed stable projection differs
from the order of Hankel-norm approximation.
Method
Let be the stable linear system
d[x(t)] = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t) (1)
where d[x(t)] is dx(t)/dt for a continuous-time system and x(t+1)
for a discrete-time system. The subroutine AB09CD determines for
the given system (1), the matrices of a reduced order system
d[z(t)] = Ar*z(t) + Br*u(t)
yr(t) = Cr*z(t) + Dr*u(t) (2)
such that
HSV(NR) <= INFNORM(G-Gr) <= 2*[HSV(NR+1) + ... + HSV(N)],
where G and Gr are transfer-function matrices of the systems
(A,B,C,D) and (Ar,Br,Cr,Dr), respectively, and INFNORM(G) is the
infinity-norm of G.
The optimal Hankel-norm approximation method of [1], based on the
square-root balancing projection formulas of [2], is employed.
References
[1] Glover, K.
All optimal Hankel norm approximation of linear
multivariable systems and their L-infinity error bounds.
Int. J. Control, Vol. 36, pp. 1145-1193, 1984.
[2] Tombs M.S. and Postlethwaite I.
Truncated balanced realization of stable, non-minimal
state-space systems.
Int. J. Control, Vol. 46, pp. 1319-1330, 1987.
Numerical Aspects
The implemented methods rely on an accuracy enhancing square-root
technique.
3
The algorithms require less than 30N floating point operations.
Further Comments
NoneExample
Program Text
* AB09CD 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 LDA, LDB, LDC, LDD
PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX,
$ LDD = PMAX )
INTEGER LIWORK
PARAMETER ( LIWORK = MAX( NMAX, MMAX ) )
INTEGER LDWORK
PARAMETER ( LDWORK = MAX( NMAX*( 2*NMAX +
$ MAX( NMAX, MMAX, PMAX ) + 5 ) +
$ ( NMAX*( NMAX + 1 ) )/2,
$ NMAX*( MMAX + PMAX + 2 ) +
$ 2*MMAX*PMAX +
$ MIN( NMAX, MMAX ) + MAX( 3*MMAX + 1,
$ MIN( NMAX, MMAX ) + PMAX ) ) )
* .. Local Scalars ..
DOUBLE PRECISION TOL1, TOL2
INTEGER I, INFO, IWARN, J, M, N, NR, P
CHARACTER*1 DICO, EQUIL, ORDSEL
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ D(LDD,MMAX), DWORK(LDWORK), HSV(NMAX)
INTEGER IWORK(LIWORK)
* .. External Subroutines ..
EXTERNAL AB09CD
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* .. 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, NR, TOL1, TOL2, DICO, EQUIL, ORDSEL
IF ( N.LT.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.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99989 ) 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 = 99988 ) 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 a reduced ssr for (A,B,C).
CALL AB09CD( DICO, EQUIL, ORDSEL, N, M, P, NR,
$ A, LDA, B, LDB, C, LDC, D, LDD, HSV, TOL1,
$ TOL2, IWORK, DWORK, LDWORK, IWARN, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 ) NR
WRITE ( NOUT, FMT = 99987 )
WRITE ( NOUT, FMT = 99995 ) ( HSV(J), J = 1,N )
IF( NR.GT.0 ) WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, NR
WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,NR )
20 CONTINUE
IF( NR.GT.0 ) WRITE ( NOUT, FMT = 99993 )
DO 40 I = 1, NR
WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,M )
40 CONTINUE
IF( NR.GT.0 ) WRITE ( NOUT, FMT = 99992 )
DO 60 I = 1, P
WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,NR )
60 CONTINUE
WRITE ( NOUT, FMT = 99991 )
DO 70 I = 1, P
WRITE ( NOUT, FMT = 99995 ) ( D(I,J), J = 1,M )
70 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' AB09CD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from AB09CD = ',I2)
99997 FORMAT (' The order of reduced model = ',I2)
99996 FORMAT (/' The reduced state dynamics matrix Ar is ')
99995 FORMAT (20(1X,F8.4))
99993 FORMAT (/' The reduced input/state matrix Br is ')
99992 FORMAT (/' The reduced state/output matrix Cr is ')
99991 FORMAT (/' The reduced input/output matrix Dr is ')
99990 FORMAT (/' N is out of range.',/' N = ',I5)
99989 FORMAT (/' M is out of range.',/' M = ',I5)
99988 FORMAT (/' P is out of range.',/' P = ',I5)
99987 FORMAT (/' The Hankel singular values are')
END
Program Data
AB09CD EXAMPLE PROGRAM DATA (Continuous system) 7 2 3 0 1.E-1 1.E-14 C N A -0.04165 0.0000 4.9200 -4.9200 0.0000 0.0000 0.0000 -5.2100 -12.500 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.3300 -3.3300 0.0000 0.0000 0.0000 0.0000 0.5450 0.0000 0.0000 0.0000 -0.5450 0.0000 0.0000 0.0000 0.0000 0.0000 4.9200 -0.04165 0.0000 4.9200 0.0000 0.0000 0.0000 0.0000 -5.2100 -12.500 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.3300 -3.3300 0.0000 0.0000 12.500 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 12.500 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000Program Results
AB09CD EXAMPLE PROGRAM RESULTS The order of reduced model = 5 The Hankel singular values are 2.5139 2.0846 1.9178 0.7666 0.5473 0.0253 0.0246 The reduced state dynamics matrix Ar is -0.5038 -5.3070 -3.2250 0.0000 0.0000 1.8355 -0.5038 -2.6289 0.0000 0.0000 0.0000 0.0000 -1.5171 0.0000 0.0000 0.0000 0.0000 0.0000 -1.2925 -9.0718 0.0000 0.0000 0.0000 0.5047 -1.2925 The reduced input/state matrix Br is -1.5343 1.5343 -0.3614 0.3614 -1.1096 1.1096 -4.5325 -4.5325 -0.7396 -0.7396 The reduced state/output matrix Cr is 1.8971 -0.3055 -2.1124 0.4421 -2.1023 -0.0394 1.1112 -0.3119 0.0000 0.0000 -1.8971 0.3055 2.1124 0.4421 -2.1023 The reduced input/output matrix Dr is 0.0126 -0.0126 0.0005 -0.0005 -0.0126 0.0126