Purpose
To compute a reduced order model (Ar,Br,Cr,Dr) for an original state-space representation (A,B,C,D) by using either the square-root or the balancing-free square-root Singular Perturbation Approximation (SPA) model reduction method in conjunction with stable coprime factorization techniques.Specification
SUBROUTINE AB09GD( DICO, JOBCF, FACT, JOBMR, EQUIL, ORDSEL, N, M,
$ P, NR, ALPHA, A, LDA, B, LDB, C, LDC, D, LDD,
$ NQ, HSV, TOL1, TOL2, TOL3, IWORK, DWORK,
$ LDWORK, IWARN, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, EQUIL, FACT, JOBCF, JOBMR, ORDSEL
INTEGER INFO, IWARN, LDA, LDB, LDC, LDD, LDWORK, M, N,
$ NQ, NR, P
DOUBLE PRECISION ALPHA, TOL1, TOL2, TOL3
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.
JOBCF CHARACTER*1
Specifies whether left or right coprime factorization is
to be used as follows:
= 'L': use left coprime factorization;
= 'R': use right coprime factorization.
FACT CHARACTER*1
Specifies the type of coprime factorization to be computed
as follows:
= 'S': compute a coprime factorization with prescribed
stability degree ALPHA;
= 'I': compute a coprime factorization with inner
denominator.
JOBMR CHARACTER*1
Specifies the model reduction approach to be used
as follows:
= 'B': use the square-root Balance & Truncate method;
= 'N': use the balancing-free square-root
Balance & Truncate method.
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(NR,NQ,NMIN), where NR
is the desired order on entry, NQ is the order of the
computed coprime factorization of the given system, and
NMIN is the order of a minimal realization of the extended
system (see METHOD); NMIN is determined as the number of
Hankel singular values greater than NQ*EPS*HNORM(Ge),
where EPS is the machine precision (see LAPACK Library
Routine DLAMCH) and HNORM(Ge) is the Hankel norm of the
extended system (computed in HSV(1));
if ORDSEL = 'A', NR is equal to the number of Hankel
singular values greater than MAX(TOL1,NQ*EPS*HNORM(Ge)).
ALPHA (input) DOUBLE PRECISION
If FACT = 'S', the desired stability degree for the
factors of the coprime factorization (see SLICOT Library
routines SB08ED/SB08FD).
ALPHA < 0 for a continuous-time system (DICO = 'C'), and
0 <= ALPHA < 1 for a discrete-time system (DICO = 'D').
If FACT = 'I', ALPHA is not used.
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, if INFO = 0, the leading NR-by-NR part of this
array contains the state dynamics matrix Ar of the reduced
order system.
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).
NQ (output) INTEGER
The order of the computed extended system Ge (see METHOD).
HSV (output) DOUBLE PRECISION array, dimension (N)
If INFO = 0, it contains the NQ Hankel singular values of
the extended system Ge ordered decreasingly (see METHOD).
Tolerances
TOL1 DOUBLE PRECISION
If ORDSEL = 'A', TOL1 contains the tolerance for
determining the order of reduced extended system.
For model reduction, the recommended value is
TOL1 = c*HNORM(Ge), where c is a constant in the
interval [0.00001,0.001], and HNORM(Ge) is the
Hankel-norm of the extended system (computed in HSV(1)).
The value TOL1 = NQ*EPS*HNORM(Ge) is used by default if
TOL1 <= 0 on entry, where EPS is the machine precision
(see LAPACK Library Routine DLAMCH).
If ORDSEL = 'F', the value of TOL1 is ignored.
TOL2 DOUBLE PRECISION
The tolerance for determining the order of a minimal
realization of the extended system Ge (see METHOD).
The recommended value is TOL2 = NQ*EPS*HNORM(Ge).
This value is used by default if TOL2 <= 0 on entry.
If TOL2 > 0, then TOL2 <= TOL1.
TOL3 DOUBLE PRECISION
The absolute tolerance level below which the elements of
B or C are considered zero (used for controllability or
observability tests by the coprime factorization method).
If the user sets TOL3 <= 0, then an implicitly computed,
default tolerance TOLDEF is used:
TOLDEF = N*EPS*NORM(C'), if JOBCF = 'L', or
TOLDEF = N*EPS*NORM(B), if JOBCF = 'R',
where EPS is the machine precision, and NORM(.) denotes
the 1-norm of a matrix.
Workspace
IWORK INTEGER array, dimension (MAX(1,2*N,PM))
where PM = P, if JOBCF = 'L',
PM = M, if JOBCF = 'R'.
On exit with INFO = 0, IWORK(1) contains the order of the
minimal realization of the system.
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,LW1) if JOBCF = 'L' and FACT = 'S',
LDWORK >= MAX(1,LW2) if JOBCF = 'L' and FACT = 'I',
LDWORK >= MAX(1,LW3) if JOBCF = 'R' and FACT = 'S',
LDWORK >= MAX(1,LW4) if JOBCF = 'R' and FACT = 'I', where
LW1 = N*(2*MAX(M,P) + P) + MAX(M,P)*(MAX(M,P) + P) +
MAX( N*P+MAX(N*(N+5), 5*P, 4*M), LWR ),
LW2 = N*(2*MAX(M,P) + P) + MAX(M,P)*(MAX(M,P) + P) +
MAX( N*P+MAX(N*(N+5), P*(P+2), 4*P, 4*M), LWR ),
LW3 = (N+M)*(M+P) + MAX( 5*M, 4*P, LWR ),
LW4 = (N+M)*(M+P) + MAX( M*(M+2), 4*M, 4*P, LWR ), and
LWR = 2*N*N + N*(MAX(N,M+P)+5) + N*(N+1)/2.
For optimum performance LDWORK should be larger.
Warning Indicator
IWARN INTEGER
= 0: no warning;
= 10*K+I:
I = 1: with ORDSEL = 'F', the selected order NR is
greater than the order of the computed coprime
factorization 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;
K > 0: K violations of the numerical stability
condition occured when computing the coprime
factorization using pole assignment (see SLICOT
Library routines SB08CD/SB08ED, SB08DD/SB08FD).
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 a real Schur form failed;
= 2: a failure was detected during the ordering of the
real Schur form of A, or in the iterative process
for reordering the eigenvalues of Z'*(A + H*C)*Z
(or Z'*(A + B*F)*Z) along the diagonal; see SLICOT
Library routines SB08CD/SB08ED (or SB08DD/SB08FD);
= 3: the matrix A has an observable or controllable
eigenvalue on the imaginary axis if DICO = 'C' or
on the unit circle if DICO = 'D';
= 4: the computation of Hankel singular values failed.
Method
Let be the 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, and let G be the corresponding
transfer-function matrix. The subroutine AB09GD determines
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)
with the transfer-function matrix Gr, by using the
singular perturbation approximation (SPA) method in conjunction
with a left coprime factorization (LCF) or a right coprime
factorization (RCF) technique:
1. Compute the appropriate stable coprime factorization of G:
-1 -1
G = R *Q (LCF) or G = Q*R (RCF).
2. Perform the model reduction algorithm on the extended system
( Q )
Ge = ( Q R ) (LCF) or Ge = ( R ) (RCF)
to obtain a reduced extended system with reduced factors
( Qr )
Ger = ( Qr Rr ) (LCF) or Ger = ( Rr ) (RCF).
3. Recover the reduced system from the reduced factors as
-1 -1
Gr = Rr *Qr (LCF) or Gr = Qr*Rr (RCF).
The approximation error for the extended system satisfies
HSV(NR) <= INFNORM(Ge-Ger) <= 2*[HSV(NR+1) + ... + HSV(NQ)],
where INFNORM(G) is the infinity-norm of G.
If JOBMR = 'B', the balancing-based square-root SPA method of [1]
is used for model reduction.
If JOBMR = 'N', the balancing-free square-root SPA method of [2]
is used for model reduction.
By setting TOL1 = TOL2, the routine can be used to compute
Balance & Truncate approximations.
If FACT = 'S', the stable coprime factorization with prescribed
stability degree ALPHA is computed by using the algorithm of [3].
If FACT = 'I', the stable coprime factorization with inner
denominator is computed by using the algorithm of [4].
References
[1] Liu Y. and Anderson B.D.O.
Singular Perturbation Approximation of Balanced Systems.
Int. J. Control, Vol. 50, pp. 1379-1405, 1989.
[2] Varga A.
Balancing-free square-root algorithm for computing singular
perturbation approximations.
Proc. 30-th IEEE CDC, Brighton, Dec. 11-13, 1991, Vol. 2,
pp. 1062-1065.
[3] Varga A.
Coprime factors model reduction method based on square-root
balancing-free techniques.
System Analysis, Modelling and Simulation, Vol. 11,
pp. 303-311, 1993.
[4] Varga A.
A Schur method for computing coprime factorizations with
inner denominators and applications in model reduction.
Proc. ACC'93, San Francisco, CA, pp. 2130-2131, 1993.
Numerical Aspects
The implemented methods rely on accuracy enhancing square-root or
balancing-free square-root techniques.
3
The algorithms require less than 30N floating point operations.
Further Comments
NoneExample
Program Text
* AB09GD 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( 2*NMAX, MMAX, PMAX ) )
* The formula below uses that NMAX = MMAX = PMAX.
INTEGER LDWORK
PARAMETER ( LDWORK = 10*NMAX*NMAX + 5*NMAX )
* .. Local Scalars ..
DOUBLE PRECISION ALPHA, TOL1, TOL2, TOL3
INTEGER I, INFO, IWARN, J, M, N, NQ, NR, P
CHARACTER*1 DICO, EQUIL, FACT, JOBCF, JOBMR, 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 AB09GD
* .. 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, NR, ALPHA, TOL1, TOL2, TOL3,
$ DICO, JOBCF, FACT, JOBMR, 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,D).
CALL AB09GD( DICO, JOBCF, FACT, JOBMR, EQUIL, ORDSEL,
$ N, M, P, NR, ALPHA, A, LDA, B, LDB, C, LDC,
$ D, LDD, NQ, HSV, TOL1, TOL2, TOL3, 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, NQ )
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 80 I = 1, P
WRITE ( NOUT, FMT = 99995 ) ( D(I,J), J = 1, M )
80 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' AB09GD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from AB09GD = ',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 of coprime factors are')
END
Program Data
AB08GD EXAMPLE PROGRAM DATA (Continuous system) 7 2 3 0 -1.e-1 .1 1.E-10 1.E-10 C L I B S A -0.04165 0.0000 4.9200 0.4920 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.0545 0.0000 0.0000 0.0000 0.0000 0.0000 -0.49200 0.004165 0.0000 4.9200 0.0000 0.0000 0.0000 0.0000 0.5210 -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
AB09GD EXAMPLE PROGRAM RESULTS The order of reduced model = 5 The Hankel singular values of coprime factors are 13.6047 9.4106 1.7684 0.7456 0.6891 0.0241 0.0230 The reduced state dynamics matrix Ar is 0.0521 -0.1491 0.0032 -0.0242 0.0181 0.2341 0.2615 0.0009 -0.0171 -0.0362 0.1170 0.0076 -0.5471 2.0904 -0.8098 0.1675 -0.2122 -2.2113 -2.4097 -0.4139 0.0390 -0.5061 0.8787 -0.3166 -3.2955 The reduced input/state matrix Br is 1.0449 0.5863 -0.0490 -1.9210 -1.3930 0.0540 -1.7206 -0.8039 0.6358 -2.0542 The reduced state/output matrix Cr is 0.4331 0.1125 -1.2534 1.6965 -0.6773 1.1171 0.3963 0.6102 0.1213 -0.0841 0.0736 -1.8815 0.0134 0.8457 2.0413 The reduced input/output matrix Dr is 0.0480 0.0003 -0.0017 0.0001 0.0005 0.0460