Purpose
To compute, for a given open-loop model (A,B,C,D), and for given state feedback gain F and full observer gain G, such that A+B*F and A+G*C are stable, a reduced order controller model (Ac,Bc,Cc) using a coprime factorization based controller reduction approach. For reduction of coprime factors, a stability enforcing frequency-weighted model reduction is performed using either the square-root or the balancing-free square-root versions of the Balance & Truncate (B&T) model reduction method.Specification
SUBROUTINE SB16CD( DICO, JOBD, JOBMR, JOBCF, ORDSEL, N, M, P, NCR,
$ A, LDA, B, LDB, C, LDC, D, LDD, F, LDF, G, LDG,
$ HSV, TOL, IWORK, DWORK, LDWORK, IWARN, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, JOBCF, JOBD, JOBMR, ORDSEL
INTEGER INFO, IWARN, LDA, LDB, LDC, LDD,
$ LDF, LDG, LDWORK, M, N, NCR, P
DOUBLE PRECISION TOL
C .. Array Arguments ..
INTEGER IWORK(*)
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
$ DWORK(*), F(LDF,*), G(LDG,*), HSV(*)
Arguments
Mode Parameters
DICO CHARACTER*1
Specifies the type of the open-loop system as follows:
= 'C': continuous-time system;
= 'D': discrete-time system.
JOBD CHARACTER*1
Specifies whether or not a non-zero matrix D appears
in the given state space model, as follows:
= 'D': D is present;
= 'Z': D is assumed a zero matrix.
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.
JOBCF CHARACTER*1
Specifies whether left or right coprime factorization
of the controller is to be used as follows:
= 'L': use left coprime factorization;
= 'R': use right coprime factorization.
ORDSEL CHARACTER*1
Specifies the order selection method as follows:
= 'F': the resulting controller order NCR is fixed;
= 'A': the resulting controller order NCR is
automatically determined on basis of the given
tolerance TOL.
Input/Output Parameters
N (input) INTEGER
The order of the original state-space representation, i.e.
the order of the matrix A. N >= 0.
N also represents the order of the original state-feedback
controller.
M (input) INTEGER
The number of system inputs. M >= 0.
P (input) INTEGER
The number of system outputs. P >= 0.
NCR (input/output) INTEGER
On entry with ORDSEL = 'F', NCR is the desired order of
the resulting reduced order controller. 0 <= NCR <= N.
On exit, if INFO = 0, NCR is the order of the resulting
reduced order controller. NCR is set as follows:
if ORDSEL = 'F', NCR is equal to MIN(NCR,NCRMIN), where
NCR is the desired order on entry, and NCRMIN is the
number of Hankel-singular values greater than N*EPS*S1,
where EPS is the machine precision (see LAPACK Library
Routine DLAMCH) and S1 is the largest Hankel singular
value (computed in HSV(1)); NCR can be further reduced
to ensure HSV(NCR) > HSV(NCR+1);
if ORDSEL = 'A', NCR is equal to the number of Hankel
singular values greater than MAX(TOL,N*EPS*S1).
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 NCR-by-NCR part of this
array contains the state dynamics matrix Ac of the reduced
controller.
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 open-loop system input/state matrix B.
On exit, this array is overwritten with a NCR-by-M
B&T approximation of the 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 open-loop system state/output matrix C.
On exit, this array is overwritten with a P-by-NCR
B&T approximation of the matrix C.
LDC INTEGER
The leading dimension of array C. LDC >= MAX(1,P).
D (input) DOUBLE PRECISION array, dimension (LDD,M)
On entry, if JOBD = 'D', the leading P-by-M part of this
array must contain the system direct input/output
transmission matrix D.
The array D is not referenced if JOBD = 'Z'.
LDD INTEGER
The leading dimension of array D.
LDD >= MAX(1,P), if JOBD = 'D';
LDD >= 1, if JOBD = 'Z'.
F (input/output) DOUBLE PRECISION array, dimension (LDF,N)
On entry, the leading M-by-N part of this array must
contain a stabilizing state feedback matrix.
On exit, if INFO = 0, the leading M-by-NCR part of this
array contains the output/state matrix Cc of the reduced
controller.
LDF INTEGER
The leading dimension of array F. LDF >= MAX(1,M).
G (input/output) DOUBLE PRECISION array, dimension (LDG,P)
On entry, the leading N-by-P part of this array must
contain a stabilizing observer gain matrix.
On exit, if INFO = 0, the leading NCR-by-P part of this
array contains the input/state matrix Bc of the reduced
controller.
LDG INTEGER
The leading dimension of array G. LDG >= MAX(1,N).
HSV (output) DOUBLE PRECISION array, dimension (N)
If INFO = 0, HSV contains the N frequency-weighted
Hankel singular values ordered decreasingly (see METHOD).
Tolerances
TOL DOUBLE PRECISION
If ORDSEL = 'A', TOL contains the tolerance for
determining the order of reduced controller.
The recommended value is TOL = c*S1, where c is a constant
in the interval [0.00001,0.001], and S1 is the largest
Hankel singular value (computed in HSV(1)).
The value TOL = N*EPS*S1 is used by default if
TOL <= 0 on entry, where EPS is the machine precision
(see LAPACK Library Routine DLAMCH).
If ORDSEL = 'F', the value of TOL is ignored.
Workspace
IWORK INTEGER array, dimension (LIWORK)
LIWORK = 0, if JOBMR = 'B';
LIWORK = N, if JOBMR = 'F'.
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*N*N + MAX( 1, 2*N*N + 5*N, N*MAX(M,P),
N*(N + MAX(N,MP) + MIN(N,MP) + 6)),
where MP = M, if JOBCF = 'L';
MP = P, if JOBCF = 'R'.
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 the order of a minimal realization
of the controller;
= 2: with ORDSEL = 'F', the selected order NCR
corresponds to repeated singular values, which are
neither all included nor all excluded from the
reduced controller. In this case, the resulting NCR
is set automatically to the largest value such that
HSV(NCR) > HSV(NCR+1).
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: eigenvalue computation failure;
= 2: the matrix A+G*C is not stable;
= 3: the matrix A+B*F is not stable;
= 4: the Lyapunov equation for computing the
observability Grammian is (nearly) singular;
= 5: the Lyapunov equation for computing the
controllability Grammian is (nearly) singular;
= 6: 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 Go(d) be the open-loop
transfer-function matrix
-1
Go(d) = C*(d*I-A) *B + D .
Let F and G be the state feedback and observer gain matrices,
respectively, chosen such that A+BF and A+GC are stable matrices.
The controller has a transfer-function matrix K(d) given by
-1
K(d) = F*(d*I-A-B*F-G*C-G*D*F) *G .
The closed-loop transfer function matrix is given by
-1
Gcl(d) = Go(d)(I+K(d)Go(d)) .
K(d) can be expressed as a left coprime factorization (LCF)
-1
K(d) = M_left(d) *N_left(d),
or as a right coprime factorization (RCF)
-1
K(d) = N_right(d)*M_right(d) ,
where M_left(d), N_left(d), N_right(d), and M_right(d) are
stable transfer-function matrices.
The subroutine SB16CD determines the matrices of a reduced
controller
d[z(t)] = Ac*z(t) + Bc*y(t)
u(t) = Cc*z(t), (2)
with the transfer-function matrix Kr, using the following
stability enforcing approach proposed in [1]:
(1) If JOBCF = 'L', the frequency-weighted approximation problem
is solved
min||[M_left(d)-M_leftr(d) N_left(d)-N_leftr(d)][-Y(d)]|| ,
[ X(d)]
where
-1
G(d) = Y(d)*X(d)
is a RCF of the open-loop system transfer-function matrix.
The B&T model reduction technique is used in conjunction
with the method proposed in [1].
(2) If JOBCF = 'R', the frequency-weighted approximation problem
is solved
min || [ -U(d) V(d) ] [ N_right(d)-N_rightr(d) ] || ,
[ M_right(d)-M_rightr(d) ]
where
-1
G(d) = V(d) *U(d)
is a LCF of the open-loop system transfer-function matrix.
The B&T model reduction technique is used in conjunction
with the method proposed in [1].
If ORDSEL = 'A', the order of the controller is determined by
computing the number of Hankel singular values greater than
the given tolerance TOL. The Hankel singular values are
the square roots of the eigenvalues of the product of
two frequency-weighted Grammians P and Q, defined as follows.
If JOBCF = 'L', then P is the controllability Grammian of a system
of the form (A+BF,B,*,*), and Q is the observability Grammian of a
system of the form (A+GC,*,F,*). This choice corresponds to an
input frequency-weighted order reduction of left coprime
factors [1].
If JOBCF = 'R', then P is the controllability Grammian of a system
of the form (A+BF,G,*,*), and Q is the observability Grammian of a
system of the form (A+GC,*,C,*). This choice corresponds to an
output frequency-weighted order reduction of right coprime
factors [1].
For the computation of truncation matrices, the B&T approach
is used in conjunction with accuracy enhancing techniques.
If JOBMR = 'B', the square-root B&T method of [2,4] is used.
If JOBMR = 'F', the balancing-free square-root version of the
B&T method [3,4] is used.
References
[1] Liu, Y., Anderson, B.D.O. and Ly, O.L.
Coprime factorization controller reduction with Bezout
identity induced frequency weighting.
Automatica, vol. 26, pp. 233-249, 1990.
[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.
[3] Varga, A.
Efficient minimal realization procedure based on balancing.
Proc. of IMACS/IFAC Symp. MCTS, Lille, France, May 1991,
A. El Moudui, P. Borne, S. G. Tzafestas (Eds.), Vol. 2,
pp. 42-46, 1991.
[4] 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.
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
* SB16CD 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, LDDC, LDF, LDG
PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX,
$ LDD = PMAX, LDDC = MMAX, LDF = MMAX, LDG = NMAX )
INTEGER LDWORK, LIWORK, MPMAX
PARAMETER ( LIWORK = 2*NMAX, MPMAX = MAX( MMAX, PMAX ) )
PARAMETER ( LDWORK = 2*NMAX*NMAX +
$ MAX( 2*NMAX*NMAX + 5*NMAX,
$ NMAX*( NMAX + MAX( NMAX, MPMAX )
$ + MIN( NMAX, MPMAX ) + 6 ) )
$ )
CHARACTER DICO, JOBCF, JOBD, JOBMR, ORDSEL
INTEGER I, INFO, IWARN, J, M, N, NCR, P
DOUBLE PRECISION TOL
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ D(LDD,MMAX), DWORK(LDWORK),
$ F(LDF,NMAX), G(LDG,PMAX), HSV(NMAX)
INTEGER IWORK(LIWORK)
* .. External Subroutines ..
EXTERNAL SB16CD
* .. 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, NCR, TOL,
$ DICO, JOBD, JOBMR, JOBCF, 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 )
READ ( NIN, FMT = * ) ( ( F(I,J), J = 1,N ), I = 1,M )
READ ( NIN, FMT = * ) ( ( G(I,J), J = 1,P ), I = 1,N )
* Find a reduced ssr for (A,B,C,D).
CALL SB16CD( DICO, JOBD, JOBMR, JOBCF, ORDSEL, N, M, P,
$ NCR, A, LDA, B, LDB, C, LDC, D, LDD, F, LDF,
$ G, LDG, HSV, TOL, IWORK, DWORK, LDWORK,
$ IWARN, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 ) NCR
WRITE ( NOUT, FMT = 99987 )
WRITE ( NOUT, FMT = 99995 ) ( HSV(J), J = 1,N )
IF( NCR.GT.0 ) WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, NCR
WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,NCR )
20 CONTINUE
IF( NCR.GT.0 ) WRITE ( NOUT, FMT = 99993 )
DO 40 I = 1, NCR
WRITE ( NOUT, FMT = 99995 ) ( G(I,J), J = 1,P )
40 CONTINUE
IF( NCR.GT.0 ) WRITE ( NOUT, FMT = 99992 )
DO 60 I = 1, M
WRITE ( NOUT, FMT = 99995 ) ( F(I,J), J = 1,NCR )
60 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' SB16CD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB16CD = ',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 ')
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 frequency-weighted Hankel singular values are:')
END
Program Data
SB16CD EXAMPLE PROGRAM DATA (Continuous system)
8 1 1 2 0.1E0 C D F R F
0 1.0000 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 -0.0150 0.7650 0 0 0 0
0 0 -0.7650 -0.0150 0 0 0 0
0 0 0 0 -0.0280 1.4100 0 0
0 0 0 0 -1.4100 -0.0280 0 0
0 0 0 0 0 0 -0.0400 1.850
0 0 0 0 0 0 -1.8500 -0.040
0.0260
-0.2510
0.0330
-0.8860
-4.0170
0.1450
3.6040
0.2800
-.996 -.105 0.261 .009 -.001 -.043 0.002 -0.026
0.0
4.472135954999638e-002 6.610515358414598e-001 4.698598960657579e-003 3.601363251422058e-001 1.032530880771415e-001 -3.754055214487997e-002 -4.268536964759344e-002 3.287284547842979e-002
4.108939884667451e-001
8.684600000000012e-002
3.852317308197148e-004
-3.619366874815911e-003
-8.803722876359955e-003
8.420521094001852e-003
1.234944428038507e-003
4.263205617645322e-003
Program Results
SB16CD EXAMPLE PROGRAM RESULTS The order of reduced controller = 2 The frequency-weighted Hankel singular values are: 3.3073 0.7274 0.1124 0.0784 0.0242 0.0182 0.0101 0.0094 The reduced controller state dynamics matrix Ac is -0.4334 0.4884 -0.1950 -0.1093 The reduced controller input/state matrix Bc is -0.4231 -0.1785 The reduced controller state/output matrix Cc is -0.0326 -0.2307