**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,Dc) using a coprime factorization based controller reduction approach. For reduction, either the square-root or the balancing-free square-root versions of the Balance & Truncate (B&T) or Singular Perturbation Approximation (SPA) model reduction methods are used in conjunction with stable coprime factorization techniques.

SUBROUTINE SB16BD( DICO, JOBD, JOBMR, JOBCF, EQUIL, ORDSEL, $ N, M, P, NCR, A, LDA, B, LDB, C, LDC, D, LDD, $ F, LDF, G, LDG, DC, LDDC, HSV, TOL1, TOL2, $ IWORK, DWORK, LDWORK, IWARN, INFO ) C .. Scalar Arguments .. CHARACTER DICO, EQUIL, JOBCF, JOBD, JOBMR, ORDSEL INTEGER INFO, IWARN, LDA, LDB, LDC, LDD, LDDC, $ LDF, LDG, LDWORK, M, N, NCR, P DOUBLE PRECISION TOL1, TOL2 C .. Array Arguments .. INTEGER IWORK(*) DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*), $ DC(LDDC,*), DWORK(*), F(LDF,*), G(LDG,*), HSV(*)

**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: = '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; = 'S': use the square-root SPA method; = 'P': use the balancing-free square-root SPA method. 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. EQUIL CHARACTER*1 Specifies whether the user wishes to perform a preliminary equilibration before performing order reduction as follows: = 'S': perform equilibration (scaling); = 'N': do not perform equilibration. 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 TOL1.

N (input) INTEGER The order of the open-loop 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,NMIN), where NCR is the desired order on entry, and NMIN is the order of a minimal realization of an extended system Ge (see METHOD); NMIN is determined as the number of Hankel singular values greater than N*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', NCR is equal to the number of Hankel singular values greater than MAX(TOL1,N*EPS*HNORM(Ge)). 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) DOUBLE PRECISION array, dimension (LDB,M) The leading N-by-M part of this array must contain the original input/state matrix B. LDB INTEGER The leading dimension of array B. LDB >= MAX(1,N). C (input) DOUBLE PRECISION array, dimension (LDC,N) The leading P-by-N part of this array must contain the original state/output matrix C. LDC INTEGER The leading dimension of array C. LDC >= MAX(1,P). D (input) DOUBLE PRECISION array, dimension (LDD,M) 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 state/output 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). DC (output) DOUBLE PRECISION array, dimension (LDDC,P) If INFO = 0, the leading M-by-P part of this array contains the input/output matrix Dc of the reduced controller. LDDC INTEGER The leading dimension of array DC. LDDC >= MAX(1,M). HSV (output) DOUBLE PRECISION array, dimension (N) If INFO = 0, it contains the N Hankel singular values of the extended system ordered decreasingly (see METHOD).

TOL1 DOUBLE PRECISION If ORDSEL = 'A', TOL1 contains the tolerance for determining the order of the 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 = N*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 coprime factorization controller (see METHOD). The recommended value is TOL2 = N*EPS*HNORM(Ge) (see METHOD). This value is used by default if TOL2 <= 0 on entry. If TOL2 > 0 and ORDSEL = 'A', then TOL2 <= TOL1.

IWORK INTEGER array, dimension (LIWORK) LIWORK = 0, if ORDSEL = 'F' and NCR = N. Otherwise, LIWORK = MAX(PM,M), if JOBCF = 'L', LIWORK = MAX(PM,P), if JOBCF = 'R', where PM = 0, if JOBMR = 'B', PM = N, if JOBMR = 'F', PM = MAX(1,2*N), if JOBMR = 'S' or 'P'. 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 >= P*N, if ORDSEL = 'F' and NCR = N. Otherwise, LDWORK >= (N+M)*(M+P) + MAX(LWR,4*M), if JOBCF = 'L', LDWORK >= (N+P)*(M+P) + MAX(LWR,4*P), if JOBCF = 'R', where LWR = MAX(1,N*(2*N+MAX(N,M+P)+5)+N*(N+1)/2). For optimum performance LDWORK should be larger.

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.

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = 1: the reduction of A+G*C to a real Schur form failed; = 2: the matrix A+G*C is not stable (if DICO = 'C'), or not convergent (if DICO = 'D'); = 3: the computation of Hankel singular values failed; = 4: the reduction of A+B*F to a real Schur form failed; = 5: the matrix A+B*F is not stable (if DICO = 'C'), or not convergent (if DICO = 'D').

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 so that A+B*F and A+G*C 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 SB16BD determines the matrices of a reduced controller d[z(t)] = Ac*z(t) + Bc*y(t) u(t) = Cc*z(t) + Dc*y(t), (2) with the transfer-function matrix Kr as follows: (1) If JOBCF = 'L', the extended system Ge(d) = [ N_left(d) M_left(d) ] is reduced to Ger(d) = [ N_leftr(d) M_leftr(d) ] by using either the B&T or SPA methods. The reduced order controller Kr(d) is computed as -1 Kr(d) = M_leftr(d) *N_leftr(d) ; (2) If JOBCF = 'R', the extended system Ge(d) = [ N_right(d) ] is reduced to [ M_right(d) ] Ger(d) = [ N_rightr(d) ] by using either the [ M_rightr(d) ] B&T or SPA methods. The reduced order controller Kr(d) is computed as -1 Kr(d) = N_rightr(d)* M_rightr(d) . If ORDSEL = 'A', the order of the controller is determined by computing the number of Hankel singular values greater than the given tolerance TOL1. The Hankel singular values are the square roots of the eigenvalues of the product of the controllability and observability Grammians of the extended system Ge. 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.

[1] 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. [2] 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. [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] 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.

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.

None

**Program Text**

* SB16BD 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, MAXMP, MPMAX PARAMETER ( LIWORK = 2*NMAX, MAXMP = MAX( MMAX, PMAX ), $ MPMAX = MMAX + PMAX ) PARAMETER ( LDWORK = ( NMAX + MAXMP )*MPMAX + $ MAX ( NMAX*( 2*NMAX + $ MAX( NMAX, MPMAX ) + 5 ) $ + ( NMAX*( NMAX + 1 ) )/2, $ 4*MAXMP ) ) CHARACTER DICO, EQUIL, JOBCF, JOBD, JOBMR, ORDSEL INTEGER I, INFO, IWARN, J, M, N, NCR, P DOUBLE PRECISION TOL1, TOL2 * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX), $ D(LDD,MMAX), DC(LDDC,PMAX), DWORK(LDWORK), $ F(LDF,NMAX), G(LDG,PMAX), HSV(NMAX) INTEGER IWORK(LIWORK) * .. External Subroutines .. EXTERNAL SB16BD * .. 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, NCR, TOL1, TOL2, $ DICO, JOBD, JOBMR, JOBCF, 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 ) 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 SB16BD( DICO, JOBD, JOBMR, JOBCF, EQUIL, ORDSEL, N, $ M, P, NCR, A, LDA, B, LDB, C, LDC, D, LDD, $ F, LDF, G, LDG, DC, LDDC, HSV, TOL1, TOL2, $ 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 WRITE ( NOUT, FMT = 99991 ) DO 80 I = 1, P WRITE ( NOUT, FMT = 99995 ) ( DC(I,J), J = 1,M ) 80 CONTINUE END IF END IF END IF END IF STOP * 99999 FORMAT (' SB16BD EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from SB16BD = ',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 extended system are:') END

SB16BD EXAMPLE PROGRAM DATA (Continuous system) 8 1 1 4 0.1E0 0.0 C D F L S 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.4721e-002 6.6105e-001 4.6986e-003 3.6014e-001 1.0325e-001 -3.7541e-002 -4.2685e-002 3.2873e-002 4.1089e-001 8.6846e-002 3.8523e-004 -3.6194e-003 -8.8037e-003 8.4205e-003 1.2349e-003 4.2632e-003

SB16BD EXAMPLE PROGRAM RESULTS The order of reduced controller = 4 The Hankel singular values of extended system are: 4.9078 4.8745 3.8455 3.7811 1.2289 1.1785 0.5176 0.1148 The reduced controller state dynamics matrix Ac is 0.5946 -0.7336 0.1914 -0.3368 0.5960 -0.0184 -0.1088 0.0207 1.2253 0.2043 0.1009 -1.4948 -0.0330 -0.0243 1.3440 0.0035 The reduced controller input/state matrix Bc is 0.0015 -0.0202 0.0159 -0.0544 The reduced controller state/output matrix Cc is 0.3534 0.0274 0.0337 -0.0320 The reduced controller input/output matrix Dc is 0.0000