AB09CD

Optimal Hankel-norm approximation based model reduction for stable systems

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

  To compute a reduced order model (Ar,Br,Cr,Dr) for a stable
  original state-space representation (A,B,C,D) by using the
  optimal Hankel-norm approximation method in conjunction with
  square-root balancing.

Specification
      SUBROUTINE AB09CD( DICO, EQUIL, ORDSEL, N, M, P, NR, A, LDA, B,
     $                   LDB, C, LDC, D, LDD, HSV, TOL1, TOL2, IWORK,
     $                   DWORK, LDWORK, IWARN, INFO )
C     .. Scalar Arguments ..
      CHARACTER         DICO, EQUIL, ORDSEL
      INTEGER           INFO, IWARN, LDA, LDB, LDC, LDD, LDWORK,
     $                  M, N, NR, P
      DOUBLE PRECISION  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.

  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(MAX(0,NR-KR+1),NMIN),
          where KR is the multiplicity of the Hankel singular value
          HSV(NR+1), NR is the desired order on entry, and NMIN is
          the order of a minimal realization of the given system;
          NMIN is determined as the number of Hankel singular values
          greater than N*EPS*HNORM(A,B,C), where EPS is the machine
          precision (see LAPACK Library Routine DLAMCH) and
          HNORM(A,B,C) is the Hankel norm of the system (computed
          in HSV(1));
          if ORDSEL = 'A', NR is equal to the number of Hankel
          singular values greater than MAX(TOL1,N*EPS*HNORM(A,B,C)).

  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 in a real Schur form.

  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).

  HSV     (output) DOUBLE PRECISION array, dimension (N)
          If INFO = 0, it contains the Hankel singular values of
          the original system ordered decreasingly. HSV(1) is the
          Hankel norm of the system.

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(A,B,C), where c is a constant in the
          interval [0.00001,0.001], and HNORM(A,B,C) is the
          Hankel-norm of the given system (computed in HSV(1)).
          For computing a minimal realization, the recommended
          value is TOL1 = N*EPS*HNORM(A,B,C), where EPS is the
          machine precision (see LAPACK Library Routine DLAMCH).
          This value is used by default if TOL1 <= 0 on entry.
          If ORDSEL = 'F', the value of TOL1 is ignored.

  TOL2    DOUBLE PRECISION
          The tolerance for determining the order of a minimal
          realization of the given system. The recommended value is
          TOL2 = N*EPS*HNORM(A,B,C). This value is used by default
          if TOL2 <= 0 on entry.
          If TOL2 > 0, then TOL2 <= TOL1.

Workspace
  IWORK   INTEGER array, dimension (LIWORK)
          LIWORK = MAX(1,M),   if DICO = 'C';
          LIWORK = MAX(1,N,M), if DICO = 'D'.
          On exit, if INFO = 0, IWORK(1) contains NMIN, the order of
          the computed minimal realization.

  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( LDW1, LDW2 ), where
          LDW1 = N*(2*N+MAX(N,M,P)+5) + N*(N+1)/2,
          LDW2 = N*(M+P+2) + 2*M*P + MIN(N,M) +
                 MAX( 3*M+1, MIN(N,M)+P ).
          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 the order of a minimal realization 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.

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 the real Schur form failed;
          = 2:  the state matrix A is not stable (if DICO = 'C')
                or not convergent (if DICO = 'D');
          = 3:  the computation of Hankel singular values failed;
          = 4:  the computation of stable projection failed;
          = 5:  the order of computed stable projection differs
                from the order of Hankel-norm approximation.

Method
  Let be the stable 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 AB09CD 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) <= INFNORM(G-Gr) <= 2*[HSV(NR+1) + ... + HSV(N)],

  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 optimal Hankel-norm approximation method of [1], based on the
  square-root balancing projection formulas of [2], is employed.

References
  [1] Glover, K.
      All optimal Hankel norm approximation of linear
      multivariable systems and their L-infinity error bounds.
      Int. J. Control, Vol. 36, pp. 1145-1193, 1984.

  [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.

Numerical Aspects
  The implemented methods rely on an accuracy enhancing square-root
  technique.
                                      3
  The algorithms require less than 30N  floating point operations.

Further Comments
  None
Example

Program Text

*     AB09CD 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( NMAX, MMAX ) )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = MAX( NMAX*( 2*NMAX +
     $                                 MAX( NMAX, MMAX, PMAX ) + 5 ) +
     $                               ( NMAX*( NMAX + 1 ) )/2,
     $                            NMAX*( MMAX + PMAX + 2 ) +
     $                            2*MMAX*PMAX +
     $                            MIN( NMAX, MMAX ) + MAX( 3*MMAX + 1,
     $                            MIN( NMAX, MMAX ) + PMAX ) ) )
*     .. Local Scalars ..
      DOUBLE PRECISION TOL1, TOL2
      INTEGER          I, INFO, IWARN, J, M, N, NR, P
      CHARACTER*1      DICO, EQUIL, 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         AB09CD
*     .. 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, NR, TOL1, TOL2, DICO, 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).
               CALL AB09CD( DICO, EQUIL, ORDSEL, N, M, P, NR,
     $                      A, LDA, B, LDB, C, LDC, D, LDD, 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,N )
                  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 (' AB09CD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from AB09CD = ',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 are')
      END
Program Data
 AB09CD EXAMPLE PROGRAM DATA (Continuous system)
  7     2     3     0   1.E-1  1.E-14      C     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.0000
Program Results
 AB09CD EXAMPLE PROGRAM RESULTS

 The order of reduced model =  5

 The Hankel singular values are
   2.5139   2.0846   1.9178   0.7666   0.5473   0.0253   0.0246

 The reduced state dynamics matrix Ar is 
  -0.5038  -5.3070  -3.2250   0.0000   0.0000
   1.8355  -0.5038  -2.6289   0.0000   0.0000
   0.0000   0.0000  -1.5171   0.0000   0.0000
   0.0000   0.0000   0.0000  -1.2925  -9.0718
   0.0000   0.0000   0.0000   0.5047  -1.2925

 The reduced input/state matrix Br is 
  -1.5343   1.5343
  -0.3614   0.3614
  -1.1096   1.1096
  -4.5325  -4.5325
  -0.7396  -0.7396

 The reduced state/output matrix Cr is 
   1.8971  -0.3055  -2.1124   0.4421  -2.1023
  -0.0394   1.1112  -0.3119   0.0000   0.0000
  -1.8971   0.3055   2.1124   0.4421  -2.1023

 The reduced input/output matrix Dr is 
   0.0126  -0.0126
   0.0005  -0.0005
  -0.0126   0.0126

Return to index