Purpose
To compute a reduced order controller (Acr,Bcr,Ccr,Dcr) for an
original state-space controller representation (Ac,Bc,Cc,Dc) by
using the frequency-weighted square-root or balancing-free
square-root Balance & Truncate (B&T) or Singular Perturbation
Approximation (SPA) model reduction methods. The algorithm tries
to minimize the norm of the frequency-weighted error
||V*(K-Kr)*W||
where K and Kr are the transfer-function matrices of the original
and reduced order controllers, respectively. V and W are special
frequency-weighting transfer-function matrices constructed
to enforce closed-loop stability and/or closed-loop performance.
If G is the transfer-function matrix of the open-loop system, then
the following weightings V and W can be used:
-1
(a) V = (I-G*K) *G, W = I - to enforce closed-loop stability;
-1
(b) V = I, W = (I-G*K) *G - to enforce closed-loop stability;
-1 -1
(c) V = (I-G*K) *G, W = (I-G*K) - to enforce closed-loop
stability and performance.
G has the state space representation (A,B,C,D).
If K is unstable, only the ALPHA-stable part of K is reduced.
Specification
SUBROUTINE SB16AD( DICO, JOBC, JOBO, JOBMR, WEIGHT, EQUIL, ORDSEL,
$ N, M, P, NC, NCR, ALPHA, A, LDA, B, LDB,
$ C, LDC, D, LDD, AC, LDAC, BC, LDBC, CC, LDCC,
$ DC, LDDC, NCS, HSVC, TOL1, TOL2, IWORK, DWORK,
$ LDWORK, IWARN, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, EQUIL, JOBC, JOBO, JOBMR, ORDSEL, WEIGHT
INTEGER INFO, IWARN, LDA, LDAC, LDB, LDBC, LDC, LDCC,
$ LDD, LDDC, LDWORK, M, N, NC, NCR, NCS, P
DOUBLE PRECISION ALPHA, TOL1, TOL2
C .. Array Arguments ..
INTEGER IWORK(*)
DOUBLE PRECISION A(LDA,*), AC(LDAC,*), B(LDB,*), BC(LDBC,*),
$ C(LDC,*), CC(LDCC,*), D(LDD,*), DC(LDDC,*),
$ DWORK(*), HSVC(*)
Arguments
Mode Parameters
DICO CHARACTER*1
Specifies the type of the original controller as follows:
= 'C': continuous-time controller;
= 'D': discrete-time controller.
JOBC CHARACTER*1
Specifies the choice of frequency-weighted controllability
Grammian as follows:
= 'S': choice corresponding to standard Enns' method [1];
= 'E': choice corresponding to the stability enhanced
modified Enns' method of [2].
JOBO CHARACTER*1
Specifies the choice of frequency-weighted observability
Grammian as follows:
= 'S': choice corresponding to standard Enns' method [1];
= 'E': choice corresponding to the stability enhanced
modified combination method of [2].
JOBMR CHARACTER*1
Specifies the model reduction approach to be used
as follows:
= 'B': use the square-root B&T method;
= 'F': use the balancing-free square-root B&T method;
= 'S': use the square-root SPA method;
= 'P': use the balancing-free square-root SPA method.
WEIGHT CHARACTER*1
Specifies the type of frequency-weighting, as follows:
= 'N': no weightings are used (V = I, W = I);
= 'O': stability enforcing left (output) weighting
-1
V = (I-G*K) *G is used (W = I);
= 'I': stability enforcing right (input) weighting
-1
W = (I-G*K) *G is used (V = I);
= 'P': stability and performance enforcing weightings
-1 -1
V = (I-G*K) *G , W = (I-G*K) are used.
EQUIL CHARACTER*1
Specifies whether the user wishes to preliminarily
equilibrate the triplets (A,B,C) and (Ac,Bc,Cc) 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 NCR is fixed;
= 'A': the resulting order NCR is automatically
determined on basis of the given tolerance TOL1.
Input/Output Parameters
N (input) INTEGER
The order of the open-loop system 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.
NC (input) INTEGER
The order of the controller state-space representation,
i.e., the order of the matrix AC. NC >= 0.
NCR (input/output) INTEGER
On entry with ORDSEL = 'F', NCR is the desired order of
the resulting reduced order controller. 0 <= NCR <= NC.
On exit, if INFO = 0, NCR is the order of the resulting
reduced order controller. For a controller with NCU
ALPHA-unstable eigenvalues and NCS ALPHA-stable
eigenvalues (NCU+NCS = NC), NCR is set as follows:
if ORDSEL = 'F', NCR is equal to
NCU+MIN(MAX(0,NCR-NCU),NCMIN), where NCR is the desired
order on entry, NCMIN is the number of frequency-weighted
Hankel singular values greater than NCS*EPS*S1, EPS is the
machine precision (see LAPACK Library Routine DLAMCH) and
S1 is the largest Hankel singular value (computed in
HSVC(1)); NCR can be further reduced to ensure
HSVC(NCR-NCU) > HSVC(NCR+1-NCU);
if ORDSEL = 'A', NCR is the sum of NCU and the number of
Hankel singular values greater than MAX(TOL1,NCS*EPS*S1).
ALPHA (input) DOUBLE PRECISION
Specifies the ALPHA-stability boundary for the eigenvalues
of the state dynamics matrix AC. For a continuous-time
controller (DICO = 'C'), ALPHA <= 0 is the boundary value
for the real parts of eigenvalues; for a discrete-time
controller (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 of the open-loop
system.
On exit, if INFO = 0 and EQUIL = 'S', the leading N-by-N
part of this array contains the scaled state dynamics
matrix of the open-loop system.
If EQUIL = 'N', this array is unchanged on exit.
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/state matrix B of the open-loop system.
On exit, if INFO = 0 and EQUIL = 'S', the leading N-by-M
part of this array contains the scaled input/state matrix
of the open-loop system.
If EQUIL = 'N', this array is unchanged on exit.
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 state/output matrix C of the open-loop system.
On exit, if INFO = 0 and EQUIL = 'S', the leading P-by-N
part of this array contains the scaled state/output matrix
of the open-loop system.
If EQUIL = 'N', this array is unchanged on exit.
LDC INTEGER
The leading dimension of array C. LDC >= MAX(1,P).
D (input) DOUBLE PRECISION array, dimension (LDD,M)
The leading P-by-M part of this array must contain the
input/output matrix D of the open-loop system.
LDD INTEGER
The leading dimension of array D. LDD >= MAX(1,P).
AC (input/output) DOUBLE PRECISION array, dimension (LDAC,NC)
On entry, the leading NC-by-NC part of this array must
contain the state dynamics matrix Ac of the original
controller.
On exit, if INFO = 0, the leading NCR-by-NCR part of this
array contains the state dynamics matrix Acr of the
reduced controller. The resulting Ac has a
block-diagonal form with two blocks.
For a system with NCU ALPHA-unstable eigenvalues and
NCS ALPHA-stable eigenvalues (NCU+NCS = NC), the leading
NCU-by-NCU block contains the unreduced part of Ac
corresponding to the ALPHA-unstable eigenvalues.
The trailing (NCR+NCS-NC)-by-(NCR+NCS-NC) block contains
the reduced part of Ac corresponding to ALPHA-stable
eigenvalues.
LDAC INTEGER
The leading dimension of array AC. LDAC >= MAX(1,NC).
BC (input/output) DOUBLE PRECISION array, dimension (LDBC,P)
On entry, the leading NC-by-P part of this array must
contain the input/state matrix Bc of the original
controller.
On exit, if INFO = 0, the leading NCR-by-P part of this
array contains the input/state matrix Bcr of the reduced
controller.
LDBC INTEGER
The leading dimension of array BC. LDBC >= MAX(1,NC).
CC (input/output) DOUBLE PRECISION array, dimension (LDCC,NC)
On entry, the leading M-by-NC part of this array must
contain the state/output matrix Cc of the original
controller.
On exit, if INFO = 0, the leading M-by-NCR part of this
array contains the state/output matrix Ccr of the reduced
controller.
LDCC INTEGER
The leading dimension of array CC. LDCC >= MAX(1,M).
DC (input/output) DOUBLE PRECISION array, dimension (LDDC,P)
On entry, the leading M-by-P part of this array must
contain the input/output matrix Dc of the original
controller.
On exit, if INFO = 0, the leading M-by-P part of this
array contains the input/output matrix Dcr of the reduced
controller.
LDDC INTEGER
The leading dimension of array DC. LDDC >= MAX(1,M).
NCS (output) INTEGER
The dimension of the ALPHA-stable part of the controller.
HSVC (output) DOUBLE PRECISION array, dimension (NC)
If INFO = 0, the leading NCS elements of this array
contain the frequency-weighted Hankel singular values,
ordered decreasingly, of the ALPHA-stable part of the
controller.
Tolerances
TOL1 DOUBLE PRECISION
If ORDSEL = 'A', TOL1 contains the tolerance for
determining the order of the reduced controller.
For model reduction, the recommended value is
TOL1 = c*S1, where c is a constant in the
interval [0.00001,0.001], and S1 is the largest
frequency-weighted Hankel singular value of the
ALPHA-stable part of the original controller
(computed in HSVC(1)).
If TOL1 <= 0 on entry, the used default value is
TOL1 = NCS*EPS*S1, where NCS is the number of
ALPHA-stable eigenvalues of Ac and 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 ALPHA-stable part of the given
controller. The recommended value is TOL2 = NCS*EPS*S1.
This value is used by default if TOL2 <= 0 on entry.
If TOL2 > 0 and ORDSEL = 'A', then TOL2 <= TOL1.
Workspace
IWORK INTEGER array, dimension (MAX(1,LIWRK1,LIWRK2))
LIWRK1 = 0, if JOBMR = 'B';
LIWRK1 = NC, if JOBMR = 'F';
LIWRK1 = 2*NC, if JOBMR = 'S' or 'P';
LIWRK2 = 0, if WEIGHT = 'N';
LIWRK2 = 2*(M+P), if WEIGHT = 'O', 'I', or 'P'.
On exit, if INFO = 0, IWORK(1) contains NCMIN, the order
of the computed minimal realization of the stable part of
the controller.
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 >= 2*NC*NC + MAX( 1, LFREQ, LSQRED ),
where
LFREQ = (N+NC)*(N+NC+2*M+2*P)+
MAX((N+NC)*(N+NC+MAX(N+NC,M,P)+7), (M+P)*(M+P+4))
if WEIGHT = 'I' or 'O' or 'P';
LFREQ = NC*(MAX(M,P)+5) if WEIGHT = 'N' and EQUIL = 'N';
LFREQ = MAX(N,NC*(MAX(M,P)+5)) if WEIGHT = 'N' and
EQUIL = 'S';
LSQRED = MAX( 1, 2*NC*NC+5*NC );
For optimum performance LDWORK should be larger.
Warning Indicator
IWARN INTEGER
= 0: no warning;
= 1: with ORDSEL = 'F', the selected order NCR 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
controller; in this case, the resulting NCR is set
equal to NSMIN;
= 2: with ORDSEL = 'F', the selected order NCR
corresponds to repeated singular values for the
ALPHA-stable part of the controller, which are
neither all included nor all excluded from the
reduced model; in this case, the resulting NCR is
automatically decreased to exclude all repeated
singular values;
= 3: with ORDSEL = 'F', the selected order NCR is less
than the order of the ALPHA-unstable part of the
given controller. In this case NCR 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 closed-loop system is not well-posed;
its feedthrough matrix is (numerically) singular;
= 2: the computation of the real Schur form of the
closed-loop state matrix failed;
= 3: the closed-loop state matrix is not stable;
= 4: the solution of a symmetric eigenproblem failed;
= 5: the computation of the ordered real Schur form of Ac
failed;
= 6: the separation of the ALPHA-stable/unstable
diagonal blocks failed because of very close
eigenvalues;
= 7: the computation of Hankel singular values failed.
Method
Let K be the transfer-function matrix of the original linear
controller
d[xc(t)] = Ac*xc(t) + Bc*y(t)
u(t) = Cc*xc(t) + Dc*y(t), (1)
where d[xc(t)] is dxc(t)/dt for a continuous-time system and
xc(t+1) for a discrete-time system. The subroutine SB16AD
determines the matrices of a reduced order controller
d[z(t)] = Acr*z(t) + Bcr*y(t)
u(t) = Ccr*z(t) + Dcr*y(t), (2)
such that the corresponding transfer-function matrix Kr minimizes
the norm of the frequency-weighted error
V*(K-Kr)*W, (3)
where V and W are special stable transfer-function matrices
chosen to enforce stability and/or performance of the closed-loop
system [3] (see description of the parameter WEIGHT).
The following procedure is used to reduce K in conjunction
with the frequency-weighted balancing approach of [2]
(see also [3]):
1) Decompose additively K, of order NC, as
K = K1 + K2,
such that K1 has only ALPHA-stable poles and K2, of order NCU,
has only ALPHA-unstable poles.
2) Compute for K1 a B&T or SPA frequency-weighted approximation
K1r of order NCR-NCU using the frequency-weighted balancing
approach of [1] in conjunction with accuracy enhancing
techniques specified by the parameter JOBMR.
3) Assemble the reduced model Kr as
Kr = K1r + K2.
For the reduction of the ALPHA-stable part, several accuracy
enhancing techniques can be employed (see [2] for details).
If JOBMR = 'B', the square-root B&T method of [1] is used.
If JOBMR = 'F', the balancing-free square-root version of the
B&T method [1] is used.
If JOBMR = 'S', the square-root version of the SPA method [2,3]
is used.
If JOBMR = 'P', the balancing-free square-root version of the
SPA method [2,3] is used.
For each of these methods, two left and right truncation matrices
are determined using the Cholesky factors of an input
frequency-weighted controllability Grammian P and an output
frequency-weighted observability Grammian Q.
P and Q are determined as the leading NC-by-NC diagonal blocks
of the controllability Grammian of K*W and of the
observability Grammian of V*K. Special techniques developed in [2]
are used to compute the Cholesky factors of P and Q directly
(see also SLICOT Library routine SB16AY).
The frequency-weighted Hankel singular values HSVC(1), ....,
HSVC(NC) are computed as the square roots of the eigenvalues
of the product P*Q.
References
[1] Enns, D.
Model reduction with balanced realizations: An error bound
and a frequency weighted generalization.
Proc. 23-th CDC, Las Vegas, pp. 127-132, 1984.
[2] Varga, A. and Anderson, B.D.O.
Square-root balancing-free methods for frequency-weighted
balancing related model reduction.
(report in preparation)
[3] Anderson, B.D.O and Liu, Y.
Controller reduction: concepts and approaches.
IEEE Trans. Autom. Control, Vol. 34, pp. 802-812, 1989.
Numerical Aspects
The implemented methods rely on accuracy enhancing square-root techniques.Further Comments
NoneExample
Program Text
* SB16AD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX, MMAX, PMAX, NCMAX
PARAMETER ( NMAX = 20, MMAX = 20, PMAX = 20,
$ NCMAX = 20 )
INTEGER MPMAX, NNCMAX
PARAMETER ( MPMAX = MMAX + PMAX, NNCMAX = NMAX + NCMAX )
INTEGER LDA, LDB, LDC, LDD, LDAC, LDBC, LDCC, LDDC
PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX,
$ LDD = PMAX, LDAC = NCMAX, LDBC = NCMAX,
$ LDCC = PMAX, LDDC = PMAX )
INTEGER LIWORK
PARAMETER ( LIWORK = 2*MAX( NCMAX, MPMAX ) )
INTEGER LDWORK
PARAMETER ( LDWORK = 2*NCMAX*NCMAX +
$ NNCMAX*( NNCMAX + 2*MPMAX ) +
$ MAX( NNCMAX*( NNCMAX +
$ MAX( NNCMAX, MMAX, PMAX ) + 7 ),
$ MPMAX*( MPMAX + 4 ) ) )
* .. Local Scalars ..
DOUBLE PRECISION ALPHA, TOL1, TOL2
INTEGER I, INFO, IWARN, J, M, N, NCR, NCS, NC, P
CHARACTER*1 DICO, EQUIL, JOBC, JOBO, JOBMR, ORDSEL, WEIGHT
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ D(LDD,MMAX), DWORK(LDWORK), HSVC(NMAX),
$ AC(LDAC,NCMAX), BC(LDBC,PMAX), CC(LDCC,NMAX),
$ DC(LDDC,PMAX)
INTEGER IWORK(LIWORK)
* .. External Subroutines ..
EXTERNAL SB16AD
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. 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, NC, NCR, ALPHA, TOL1, TOL2, DICO,
$ JOBC, JOBO, JOBMR, WEIGHT, EQUIL, ORDSEL
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), J = 1,M ), I = 1, N )
IF( 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 )
IF( NC.LT.0 .OR. NC.GT.NCMAX ) THEN
WRITE ( NOUT, FMT = 99986 ) NC
ELSE
IF( NC.GT.0 ) THEN
READ ( NIN, FMT = * )
$ ( ( AC(I,J), J = 1,NC ), I = 1,NC )
READ ( NIN, FMT = * )
$ ( ( BC(I,J), J = 1,P ), I = 1, NC )
READ ( NIN, FMT = * )
$ ( ( CC(I,J), J = 1,NC ), I = 1,M )
END IF
READ ( NIN, FMT = * )
$ ( ( DC(I,J), J = 1,P ), I = 1,M )
END IF
* Find a reduced ssr for (AC,BC,CC,DC).
CALL SB16AD( DICO, JOBC, JOBO, JOBMR, WEIGHT, EQUIL,
$ ORDSEL, N, M, P, NC, NCR, ALPHA, A, LDA,
$ B, LDB, C, LDC, D, LDD, AC, LDAC, BC, LDBC,
$ CC, LDCC, DC, LDDC, NCS, HSVC, TOL1, TOL2,
$ IWORK, DWORK, LDWORK, IWARN, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
IF( IWARN.NE.0) WRITE ( NOUT, FMT = 99984 ) IWARN
WRITE ( NOUT, FMT = 99997 ) NCR
WRITE ( NOUT, FMT = 99987 )
WRITE ( NOUT, FMT = 99995 ) ( HSVC(J), J = 1, NCS )
IF( NCR.GT.0 ) WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, NCR
WRITE ( NOUT, FMT = 99995 ) ( AC(I,J), J = 1,NCR )
20 CONTINUE
IF( NCR.GT.0 ) WRITE ( NOUT, FMT = 99993 )
DO 40 I = 1, NCR
WRITE ( NOUT, FMT = 99995 ) ( BC(I,J), J = 1,P )
40 CONTINUE
IF( NCR.GT.0 ) WRITE ( NOUT, FMT = 99992 )
DO 60 I = 1, M
WRITE ( NOUT, FMT = 99995 ) ( CC(I,J), J = 1,NCR )
60 CONTINUE
WRITE ( NOUT, FMT = 99991 )
DO 70 I = 1, M
WRITE ( NOUT, FMT = 99995 ) ( DC(I,J), J = 1,P )
70 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' SB16AD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB16AD = ',I2)
99997 FORMAT (/' The order of reduced controller = ',I2)
99996 FORMAT (/' The reduced controller state dynamics matrix Ac is ')
99995 FORMAT (20(1X,F8.4))
99993 FORMAT (/' The reduced controller input/state matrix Bc is ')
99992 FORMAT (/' The reduced controller state/output matrix Cc is ')
99991 FORMAT (/' The reduced controller input/output matrix Dc 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 weighted ALPHA-stable',
$ ' part are')
99986 FORMAT (/' NC is out of range.',/' NC = ',I5)
99984 FORMAT (' IWARN on exit from SB16AD = ',I2)
END
Program Data
SB16AD EXAMPLE PROGRAM DATA (Continuous system)
3 1 1 3 2 0.0 0.1E0 0.0 C S S F I N F
-1. 0. 4.
0. 2. 0.
0. 0. -3.
1.
1.
1.
1. 1. 1.
0.
-26.4000 6.4023 4.3868
32.0000 0 0
0 8.0000 0
-16
0
0
9.2994 1.1624 0.1090
0
Program Results
SB16AD EXAMPLE PROGRAM RESULTS The order of reduced controller = 2 The Hankel singular values of weighted ALPHA-stable part are 3.8253 0.2005 The reduced controller state dynamics matrix Ac is 9.1900 0.0000 0.0000 -34.5297 The reduced controller input/state matrix Bc is -11.9593 86.3137 The reduced controller state/output matrix Cc is 2.8955 -1.3566 The reduced controller input/output matrix Dc is 0.0000