AB09IY

Cholesky factors of the frequency-weighted controllability and observability Grammians

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

Purpose

  To compute for given state-space representations
  (A,B,C,0), (AV,BV,CV,DV), and (AW,BW,CW,DW) of the
  transfer-function matrices G, V and W, respectively,
  the Cholesky factors of the frequency-weighted
  controllability and observability Grammians corresponding
  to a frequency-weighted model reduction problem.
  G, V and W must be stable transfer-function matrices with
  the state matrices A, AV, and AW in real Schur form.
  It is assumed that the state space realizations (AV,BV,CV,DV)
  and (AW,BW,CW,DW) are minimal. In case of possible pole-zero
  cancellations in forming V*G and/or G*W, the parameters for the
  choice of frequency-weighted Grammians ALPHAO and/or ALPHAC,
  respectively, must be different from 1.

Specification
      SUBROUTINE AB09IY( DICO, JOBC, JOBO, WEIGHT, N, M, P, NV, PV,
     $                   NW, MW, ALPHAC, ALPHAO, A, LDA, B, LDB, C, LDC,
     $                   AV, LDAV, BV, LDBV, CV, LDCV, DV, LDDV,
     $                   AW, LDAW, BW, LDBW, CW, LDCW, DW, LDDW,
     $                   SCALEC, SCALEO, S, LDS, R, LDR,
     $                   DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER        DICO, JOBC, JOBO, WEIGHT
      INTEGER          INFO, LDA, LDAV, LDAW, LDB, LDBV, LDBW,
     $                 LDC, LDCV, LDCW, LDDV, LDDW, LDR, LDS, LDWORK,
     $                 M, MW, N, NV, NW, P, PV
      DOUBLE PRECISION ALPHAC, ALPHAO, SCALEC, SCALEO
C     .. Array Arguments ..
      DOUBLE PRECISION A(LDA,*), AV(LDAV,*), AW(LDAW,*),
     $                 B(LDB,*), BV(LDBV,*), BW(LDBW,*),
     $                 C(LDC,*), CV(LDCV,*), CW(LDCW,*),
     $                           DV(LDDV,*), DW(LDDW,*),
     $                 DWORK(*), R(LDR,*),   S(LDS,*)

Arguments

Mode Parameters

  DICO    CHARACTER*1
          Specifies the type of the systems as follows:
          = 'C':  G, V and W are continuous-time systems;
          = 'D':  G, V and W are discrete-time systems.

  JOBC    CHARACTER*1
          Specifies the choice of frequency-weighted controllability
          Grammian as follows:
          = 'S': choice corresponding to a combination method [4]
                 of the approaches of Enns [1] and Lin-Chiu [2,3];
          = 'E': choice corresponding to the stability enhanced
                 modified combination method of [4].

  JOBO    CHARACTER*1
          Specifies the choice of frequency-weighted observability
          Grammian as follows:
          = 'S': choice corresponding to a combination method [4]
                 of the approaches of Enns [1] and Lin-Chiu [2,3];
          = 'E': choice corresponding to the stability enhanced
                 modified combination method of [4].

  WEIGHT  CHARACTER*1
          Specifies the type of frequency weighting, as follows:
          = 'N':  no weightings are used (V = I, W = I);
          = 'L':  only left weighting V is used (W = I);
          = 'R':  only right weighting W is used (V = I);
          = 'B':  both left and right weightings V and W are used.

Input/Output Parameters
  N       (input) INTEGER
          The order of the state-space representation of G, i.e.,
          the order of the matrix A.  N >= 0.

  M       (input) INTEGER
          The number of columns of the matrix B and
          the number of rows of the matrices CW and DW.  M >= 0.
          M represents the dimension of the input vector of the
          system with the transfer-function matrix G and
          also the dimension of the output vector of the system
          with the transfer-function matrix W.

  P       (input) INTEGER
          The number of rows of the matrix C and the
          number of columns of the matrices BV and DV.  P >= 0.
          P represents the dimension of the output vector of the
          system with the transfer-function matrix G and
          also the dimension of the input vector of the system
          with the transfer-function matrix V.

  NV      (input) INTEGER
          The order of the matrix AV. Also the number of rows of
          the matrix BV and the number of columns of the matrix CV.
          NV represents the dimension of the state vector of the
          system with the transfer-function matrix V.  NV >= 0.

  PV      (input) INTEGER
          The number of rows of the matrices CV and DV.  PV >= 0.
          PV represents the dimension of the output vector of the
          system with the transfer-function matrix V.

  NW      (input) INTEGER
          The order of the matrix AW. Also the number of rows of
          the matrix BW and the number of columns of the matrix CW.
          NW represents the dimension of the state vector of the
          system with the transfer-function matrix W.  NW >= 0.

  MW      (input) INTEGER
          The number of columns of the matrices BW and DW.  MW >= 0.
          MW represents the dimension of the input vector of the
          system with the transfer-function matrix W.

  ALPHAC  (input) DOUBLE PRECISION
          Combination method parameter for defining the
          frequency-weighted controllability Grammian (see METHOD);
          ABS(ALPHAC) <= 1.

  ALPHAO  (input) DOUBLE PRECISION
          Combination method parameter for defining the
          frequency-weighted observability Grammian (see METHOD);
          ABS(ALPHAO) <= 1.

  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
          The leading N-by-N part of this array must
          contain the state matrix A (of the system with the
          transfer-function matrix G) in a real Schur form.

  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
          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
          state/output matrix C.

  LDC     INTEGER
          The leading dimension of array C.  LDC >= MAX(1,P).

  AV      (input) DOUBLE PRECISION array, dimension (LDAV,NV)
          If WEIGHT = 'L' or 'B', the leading NV-by-NV part of this
          array must contain the state matrix AV (of the system with
          the transfer-function matrix V) in a real Schur form.
          AV is not referenced if WEIGHT = 'R' or 'N'.

  LDAV    INTEGER
          The leading dimension of array AV.
          LDAV >= MAX(1,NV), if WEIGHT = 'L' or 'B';
          LDAV >= 1,         if WEIGHT = 'R' or 'N'.

  BV      (input) DOUBLE PRECISION array, dimension (LDBV,P)
          If WEIGHT = 'L' or 'B', the leading NV-by-P part of this
          array must contain the input matrix BV of the system with
          the transfer-function matrix V.
          BV is not referenced if WEIGHT = 'R' or 'N'.

  LDBV    INTEGER
          The leading dimension of array BV.
          LDBV >= MAX(1,NV), if WEIGHT = 'L' or 'B';
          LDBV >= 1,         if WEIGHT = 'R' or 'N'.

  CV      (input) DOUBLE PRECISION array, dimension (LDCV,NV)
          If WEIGHT = 'L' or 'B', the leading PV-by-NV part of this
          array must contain the output matrix CV of the system with
          the transfer-function matrix V.
          CV is not referenced if WEIGHT = 'R' or 'N'.

  LDCV    INTEGER
          The leading dimension of array CV.
          LDCV >= MAX(1,PV), if WEIGHT = 'L' or 'B';
          LDCV >= 1,         if WEIGHT = 'R' or 'N'.

  DV      (input) DOUBLE PRECISION array, dimension (LDDV,P)
          If WEIGHT = 'L' or 'B', the leading PV-by-P part of this
          array must contain the feedthrough matrix DV of the system
          with the transfer-function matrix V.
          DV is not referenced if WEIGHT = 'R' or 'N'.

  LDDV    INTEGER
          The leading dimension of array DV.
          LDDV >= MAX(1,PV), if WEIGHT = 'L' or 'B';
          LDDV >= 1,         if WEIGHT = 'R' or 'N'.

  AW      (input) DOUBLE PRECISION array, dimension (LDAW,NW)
          If WEIGHT = 'R' or 'B', the leading NW-by-NW part of this
          array must contain the state matrix AW (of the system with
          the transfer-function matrix W) in a real Schur form.
          AW is not referenced if WEIGHT = 'L' or 'N'.

  LDAW    INTEGER
          The leading dimension of array AW.
          LDAW >= MAX(1,NW), if WEIGHT = 'R' or 'B';
          LDAW >= 1,         if WEIGHT = 'L' or 'N'.

  BW      (input) DOUBLE PRECISION array, dimension (LDBW,MW)
          If WEIGHT = 'R' or 'B', the leading NW-by-MW part of this
          array must contain the input matrix BW of the system with
          the transfer-function matrix W.
          BW is not referenced if WEIGHT = 'L' or 'N'.

  LDBW    INTEGER
          The leading dimension of array BW.
          LDBW >= MAX(1,NW), if WEIGHT = 'R' or 'B';
          LDBW >= 1,         if WEIGHT = 'L' or 'N'.

  CW      (input) DOUBLE PRECISION array, dimension (LDCW,NW)
          If WEIGHT = 'R' or 'B', the leading M-by-NW part of this
          array must contain the output matrix CW of the system with
          the transfer-function matrix W.
          CW is not referenced if WEIGHT = 'L' or 'N'.

  LDCW    INTEGER
          The leading dimension of array CW.
          LDCW >= MAX(1,M), if WEIGHT = 'R' or 'B';
          LDCW >= 1,        if WEIGHT = 'L' or 'N'.

  DW      (input) DOUBLE PRECISION array, dimension (LDDW,MW)
          If WEIGHT = 'R' or 'B', the leading M-by-MW part of this
          array must contain the feedthrough matrix DW of the system
          with the transfer-function matrix W.
          DW is not referenced if WEIGHT = 'L' or 'N'.

  LDDW    INTEGER
          The leading dimension of array DW.
          LDDW >= MAX(1,M), if WEIGHT = 'R' or 'B';
          LDDW >= 1,        if WEIGHT = 'L' or 'N'.

  SCALEC  (output) DOUBLE PRECISION
          Scaling factor for the controllability Grammian in (1)
          or (3). See METHOD.

  SCALEO  (output) DOUBLE PRECISION
          Scaling factor for the observability Grammian in (2)
          or (4). See METHOD.

  S       (output) DOUBLE PRECISION array, dimension (LDS,N)
          The leading N-by-N upper triangular part of this array
          contains the Cholesky factor S of the frequency-weighted
          cotrollability Grammian P = S*S'. See METHOD.

  LDS     INTEGER
          The leading dimension of array S.  LDS >= MAX(1,N).

  R       (output) DOUBLE PRECISION array, dimension (LDR,N)
          The leading N-by-N upper triangular part of this array
          contains the Cholesky factor R of the frequency-weighted
          observability Grammian Q = R'*R. See METHOD.

  LDR     INTEGER
          The leading dimension of array R.  LDR >= MAX(1,N).

Workspace
  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, LLEFT, LRIGHT ),
          where
          LLEFT  = (N+NV)*(N+NV+MAX(N+NV,PV)+5)
                           if WEIGHT = 'L' or 'B' and PV > 0;
          LLEFT  = N*(P+5) if WEIGHT = 'R' or 'N' or  PV = 0;
          LRIGHT = (N+NW)*(N+NW+MAX(N+NW,MW)+5)
                           if WEIGHT = 'R' or 'B' and MW > 0;
          LRIGHT = N*(M+5) if WEIGHT = 'L' or 'N' or  MW = 0.
          For optimum performance LDWORK should be larger.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  if the state matrices A and/or AV are not stable or
                not in a real Schur form;
          = 2:  if the state matrices A and/or AW are not stable or
                not in a real Schur form;
          = 3:  eigenvalues computation failure.

Method
  Let Pi = Si*Si' and Qo = Ro'*Ro be the Cholesky factored
  controllability and observability Grammians satisfying
  in the continuous-time case

         Ai*Pi + Pi*Ai' +  scalec^2*Bi*Bi' = 0,       (1)

         Ao'*Qo + Qo*Ao +  scaleo^2*Co'*Co = 0,       (2)

  and in the discrete-time case

         Ai*Pi*Ai' - Pi +  scalec^2*Bi*Bi' = 0,       (3)

         Ao'*Qo*Ao - Qo +  scaleo^2*Co'*Co = 0,       (4)

  where

        Ai = ( A  B*Cw ) ,   Bi = ( B*Dw ) ,
             ( 0   Aw  )          (  Bw  )

        Ao = (  A   0  ) ,   Co = ( Dv*C  Cv ) .
             ( Bv*C Av )

  Consider the partitioned Grammians

        Pi = ( P11  P12 )   and    Qo = ( Q11  Q12 ) ,
             ( P12' P22 )               ( Q12' Q22 )

  where P11 and Q11 are the leading N-by-N parts of Pi and Qo,
  respectively, and let P0 and Q0 be non-negative definite matrices
  defined in the combination method [4]
                                     -1
         P0 = P11 - ALPHAC**2*P12*P22 *P21 ,
                                     -1
         Q0 = Q11 - ALPHAO**2*Q12*Q22 *Q21.

  The frequency-weighted controllability and observability
  Grammians, P and Q, respectively, are defined as follows:
  P = P0 if JOBC = 'S' (standard combination method [4]);
  P = P1 >= P0 if JOBC = 'E', where P1 is the controllability
  Grammian defined to enforce stability for a modified combination
  method of [4];
  Q = Q0 if JOBO = 'S' (standard combination method [4]);
  Q = Q1 >= Q0 if JOBO = 'E', where Q1 is the observability
  Grammian defined to enforce stability for a modified combination
  method of [4].

  If JOBC = JOBO = 'S' and ALPHAC = ALPHAO = 0, the choice of
  Grammians corresponds to the method of Enns [1], while if
  ALPHAC = ALPHAO = 1, the choice of Grammians corresponds to the
  method of Lin and Chiu [2,3].

  The routine computes directly the Cholesky factors S and R
  such that P = S*S' and Q = R'*R according to formulas
  developed in [4]. No matrix inversions are involved.

References
  [1] Enns, D.
      Model reduction with balanced realizations: An error bound
      and a frequency weighted generalization.
      Proc. CDC, Las Vegas, pp. 127-132, 1984.

  [2] Lin, C.-A. and Chiu, T.-Y.
      Model reduction via frequency-weighted balanced realization.
      Control Theory and Advanced Technology, vol. 8,
      pp. 341-351, 1992.

  [3] Sreeram, V., Anderson, B.D.O and Madievski, A.G.
      New results on frequency weighted balanced reduction
      technique.
      Proc. ACC, Seattle, Washington, pp. 4004-4009, 1995.

  [4] Varga, A. and Anderson, B.D.O.
      Square-root balancing-free methods for the frequency-weighted
      balancing related model reduction.
      (report in preparation)

Further Comments
  None
Example

Program Text

  None
Program Data
  None
Program Results
  None

Return to Supporting Routines index