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 for the ALPHA-stable part of the system.Specification
SUBROUTINE AB09ND( DICO, JOB, EQUIL, ORDSEL, N, M, P, NR, ALPHA,
$ A, LDA, B, LDB, C, LDC, D, LDD, NS, HSV, TOL1,
$ TOL2, IWORK, DWORK, LDWORK, IWARN, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, EQUIL, JOB, ORDSEL
INTEGER INFO, IWARN, LDA, LDB, LDC, LDD, LDWORK,
$ M, N, NR, NS, P
DOUBLE PRECISION ALPHA, 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.
JOB CHARACTER*1
Specifies the model reduction approach to be used
as follows:
= 'B': use the square-root SPA method;
= 'N': use the balancing-free square-root SPA 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. For a system with NU ALPHA-unstable
eigenvalues and NS ALPHA-stable eigenvalues (NU+NS = N),
NR is set as follows: if ORDSEL = 'F', NR is equal to
NU+MIN(MAX(0,NR-NU),NMIN), where NR is the desired order
on entry, and NMIN is the order of a minimal realization
of the ALPHA-stable part of the given system; NMIN is
determined as the number of Hankel singular values greater
than NS*EPS*HNORM(As,Bs,Cs), where EPS is the machine
precision (see LAPACK Library Routine DLAMCH) and
HNORM(As,Bs,Cs) is the Hankel norm of the ALPHA-stable
part of the given system (computed in HSV(1));
if ORDSEL = 'A', NR is the sum of NU and the number of
Hankel singular values greater than
MAX(TOL1,NS*EPS*HNORM(As,Bs,Cs)).
ALPHA (input) DOUBLE PRECISION
Specifies the ALPHA-stability boundary for the eigenvalues
of the state dynamics matrix A. For a continuous-time
system (DICO = 'C'), ALPHA <= 0 is the boundary value for
the real parts of eigenvalues, while for a discrete-time
system (DICO = 'D'), 0 <= ALPHA <= 1 represents the
boundary value for the moduli of eigenvalues.
The ALPHA-stability domain does not include the boundary.
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.
The resulting A has a block-diagonal form with two blocks.
For a system with NU ALPHA-unstable eigenvalues and
NS ALPHA-stable eigenvalues (NU+NS = N), the leading
NU-by-NU block contains the unreduced part of A
corresponding to ALPHA-unstable eigenvalues in an
upper real Schur form.
The trailing (NR+NS-N)-by-(NR+NS-N) block contains
the reduced part of A corresponding to ALPHA-stable
eigenvalues.
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).
NS (output) INTEGER
The dimension of the ALPHA-stable subsystem.
HSV (output) DOUBLE PRECISION array, dimension (N)
If INFO = 0, the leading NS elements of HSV contain the
Hankel singular values of the ALPHA-stable part of the
original system ordered decreasingly.
HSV(1) is the Hankel norm of the ALPHA-stable subsystem.
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(As,Bs,Cs), where c is a constant in the
interval [0.00001,0.001], and HNORM(As,Bs,Cs) is the
Hankel-norm of the ALPHA-stable part of the given system
(computed in HSV(1)).
If TOL1 <= 0 on entry, the used default value is
TOL1 = NS*EPS*HNORM(As,Bs,Cs), where NS is the number of
ALPHA-stable eigenvalues of A and EPS is the machine
precision (see LAPACK Library Routine DLAMCH).
This value is appropriate to compute a minimal realization
of the ALPHA-stable part.
If ORDSEL = 'F', the value of TOL1 is ignored.
TOL2 DOUBLE PRECISION
The tolerance for determining the order of a minimal
realization of the ALPHA-stable part of the given system.
The recommended value is TOL2 = NS*EPS*HNORM(As,Bs,Cs).
This value is used by default if TOL2 <= 0 on entry.
If TOL2 > 0, then TOL2 <= TOL1.
Workspace
IWORK INTEGER array, dimension (MAX(1,2*N))
On exit, if INFO = 0, IWORK(1) contains the order of the
minimal realization of the ALPHA-stable part 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,N*(2*N+MAX(N,M,P)+5) + N*(N+1)/2).
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 NSMIN, the sum of the order of the
ALPHA-unstable part and the order of a minimal
realization of the ALPHA-stable part of the given
system. In this case, the resulting NR is set equal
to NSMIN.
= 2: with ORDSEL = 'F', the selected order NR is less
than the order of the ALPHA-unstable part of the
given system. In this case NR is set equal to the
order of the ALPHA-unstable part.
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: the computation of the ordered real Schur form of A
failed;
= 2: the separation of the ALPHA-stable/unstable diagonal
blocks failed because of very close eigenvalues;
= 3: the computation of Hankel singular values failed.
Method
Let be the following 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 AB09ND 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+NS-N) <= INFNORM(G-Gr) <= 2*[HSV(NR+NS-N+1)+...+HSV(NS)],
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 following procedure is used to reduce a given G:
1) Decompose additively G as
G = G1 + G2
such that G1 = (As,Bs,Cs,D) has only ALPHA-stable poles and
G2 = (Au,Bu,Cu,0) has only ALPHA-unstable poles.
2) Determine G1r, a reduced order approximation of the
ALPHA-stable part G1.
3) Assemble the reduced model Gr as
Gr = G1r + G2.
To reduce the ALPHA-stable part G1, if JOB = 'B', the square-root
balancing-based SPA method of [1] is used, and for an ALPHA-stable
system, the resulting reduced model is balanced.
If JOB = 'N', the balancing-free square-root SPA method of [2]
is used to reduce the ALPHA-stable part G1.
By setting TOL1 = TOL2, the routine can be used to compute
Balance & Truncate approximations as well.
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.
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
* AB09ND 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 = 2*NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = NMAX*( 2*NMAX +
$ MAX( NMAX, MMAX, PMAX ) + 5 ) +
$ ( NMAX*( NMAX + 1 ) )/2 )
* .. Local Scalars ..
DOUBLE PRECISION ALPHA, TOL1, TOL2
INTEGER I, INFO, IWARN, J, M, N, NR, NS, P
CHARACTER*1 DICO, EQUIL, JOB, 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 AB09ND
* .. 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,
$ DICO, JOB, 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 AB09ND( DICO, JOB, EQUIL, ORDSEL, N, M, P, NR,
$ ALPHA, A, LDA, B, LDB, C, LDC, D, LDD,
$ NS, 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,NS )
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 (' AB09ND EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from AB09ND = ',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 ALPHA-stable part are')
END
Program Data
AB09ND EXAMPLE PROGRAM DATA (Continuous system) 7 2 3 0 -.6D0 1.D-1 1.E-14 C N 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
AB09ND EXAMPLE PROGRAM RESULTS The order of reduced model = 5 The Hankel singular values of ALPHA-stable part are 1.9178 0.8621 0.7666 0.0336 0.0246 The reduced state dynamics matrix Ar is -0.5181 -1.1084 0.0000 0.0000 0.0000 8.8157 -0.5181 0.0000 0.0000 0.0000 0.0000 0.0000 0.5847 0.0000 1.9230 0.0000 0.0000 0.0000 -1.6606 0.0000 0.0000 0.0000 -4.3823 0.0000 -3.2922 The reduced input/state matrix Br is -1.2837 1.2837 -0.7522 0.7522 -0.6379 -0.6379 2.0656 -2.0656 -3.9315 -3.9315 The reduced state/output matrix Cr is -0.1380 -0.6445 -0.6416 -0.6293 0.2526 0.6246 0.0196 0.0000 0.4107 0.0000 0.1380 0.6445 -0.6416 0.6293 0.2526 The reduced input/output matrix Dr is 0.0582 -0.0090 0.0015 -0.0015 -0.0090 0.0582