IB01PX

Estimating system matrices B and D using Kronecker products

Purpose

  To build and solve the least squares problem  T*X = Kv,  and
  estimate the matrices B and D of a linear time-invariant (LTI)
  state space model, using the solution  X,  and the singular
  value decomposition information and other intermediate results,
  provided by other routines.

  The matrix  T  is computed as a sum of Kronecker products,

     T = T + kron(Uf(:,(i-1)*m+1:i*m),N_i),  for i = 1 : s,

  (with  T  initialized by zero), where  Uf  is the triangular
  factor of the QR factorization of the future input part (see
  SLICOT Library routine IB01ND),  N_i  is given by the i-th block
  row of the matrix

         [ Q_11  Q_12  ...  Q_1,s-2  Q_1,s-1  Q_1s ]   [ I_L  0  ]
         [ Q_12  Q_13  ...  Q_1,s-1    Q_1s    0   ]   [         ]
     N = [ Q_13  Q_14  ...    Q_1s      0      0   ] * [         ],
         [  :     :            :        :      :   ]   [         ]
         [ Q_1s   0    ...     0        0      0   ]   [  0  GaL ]

  and where

            [   -L_1|1    ]          [ M_i-1 - L_1|i ]
     Q_11 = [             ],  Q_1i = [               ],  i = 2:s,
            [ I_L - L_2|1 ]          [     -L_2|i    ]

  are  (n+L)-by-L  matrices, and  GaL  is built from the first  n
  relevant singular vectors,  GaL = Un(1:L(s-1),1:n),  computed
  by IB01ND.

  The vector  Kv  is vec(K), with the matrix  K  defined by

     K = [ K_1  K_2  K_3  ...  K_s ],

  where  K_i = K(:,(i-1)*m+1:i*m),  i = 1:s,  is  (n+L)-by-m.
  The given matrices are  Uf,  GaL,  and

         [ L_1|1  ...  L_1|s ]
     L = [                   ],   (n+L)-by-L*s,
         [ L_2|1  ...  L_2|s ]

     M = [ M_1  ...  M_s-1 ],      n-by-L*(s-1),  and
     K,                            (n+L)-by-m*s.

  Matrix  M  is the pseudoinverse of the matrix  GaL,  computed by
  SLICOT Library routine IB01PD.

      SUBROUTINE IB01PX( JOB, NOBR, N, M, L, UF, LDUF, UN, LDUN, UL,
     $                   LDUL, PGAL, LDPGAL, K, LDK, R, LDR, X, B, LDB,
     $                   D, LDD, TOL, IWORK, DWORK, LDWORK, IWARN,
     $                   INFO )
C     .. Scalar Arguments ..
      DOUBLE PRECISION   TOL
      INTEGER            INFO, IWARN, L, LDB, LDD, LDK, LDPGAL, LDR,
     $                   LDUF, LDUL, LDUN, LDWORK, M, N, NOBR
      CHARACTER          JOB
C     .. Array Arguments ..
      DOUBLE PRECISION   B(LDB, *), D(LDD, *), DWORK(*), K(LDK, *),
     $                   PGAL(LDPGAL, *), R(LDR, *), UF(LDUF, *),
     $                   UL(LDUL, *), UN(LDUN, *), X(*)
      INTEGER            IWORK( * )

Mode Parameters

  JOB     CHARACTER*1
          Specifies which of the matrices B and D should be
          computed, as follows:
          = 'B':  compute the matrix B, but not the matrix D;
          = 'D':  compute both matrices B and D.

  NOBR    (input) INTEGER
          The number of block rows,  s,  in the input and output
          Hankel matrices processed by other routines.  NOBR > 1.

  N       (input) INTEGER
          The order of the system.  NOBR > N > 0.

  M       (input) INTEGER
          The number of system inputs.  M >= 0.

  L       (input) INTEGER
          The number of system outputs.  L > 0.

  UF      (input/output) DOUBLE PRECISION array, dimension
          ( LDUF,M*NOBR )
          On entry, the leading  M*NOBR-by-M*NOBR  upper triangular
          part of this array must contain the upper triangular
          factor of the QR factorization of the future input part,
          as computed by SLICOT Library routine IB01ND.
          The strict lower triangle need not be set to zero.
          On exit, the leading  M*NOBR-by-M*NOBR  upper triangular
          part of this array is unchanged, and the strict lower
          triangle is set to zero.

  LDUF    INTEGER
          The leading dimension of the array  UF.
          LDUF >= MAX( 1, M*NOBR ).

  UN      (input) DOUBLE PRECISION array, dimension ( LDUN,N )
          The leading  L*(NOBR-1)-by-N  part of this array must
          contain the matrix  GaL,  i.e., the leading part of the
          first  N  columns of the matrix  Un  of relevant singular
          vectors.

  LDUN    INTEGER
          The leading dimension of the array  UN.
          LDUN >= L*(NOBR-1).

  UL      (input/output) DOUBLE PRECISION array, dimension
          ( LDUL,L*NOBR )
          On entry, the leading  (N+L)-by-L*NOBR  part of this array
          must contain the given matrix  L.
          On exit, if  M > 0,  the leading  (N+L)-by-L*NOBR  part of
          this array is overwritten by the matrix
          [ Q_11  Q_12  ...  Q_1,s-2  Q_1,s-1  Q_1s ].

  LDUL    INTEGER
          The leading dimension of the array  UL.  LDUL >= N+L.

  PGAL    (input) DOUBLE PRECISION array, dimension
          ( LDPGAL,L*(NOBR-1) )
          The leading  N-by-L*(NOBR-1)  part of this array must
          contain the pseudoinverse of the matrix  GaL,  computed by
          SLICOT Library routine IB01PD.

  LDPGAL  INTEGER
          The leading dimension of the array  PGAL.  LDPGAL >= N.

  K       (input) DOUBLE PRECISION array, dimension ( LDK,M*NOBR )
          The leading  (N+L)-by-M*NOBR  part of this array must
          contain the given matrix  K.

  LDK     INTEGER
          The leading dimension of the array  K.  LDK >= N+L.

  R       (output) DOUBLE PRECISION array, dimension ( LDR,M*(N+L) )
          The leading  (N+L)*M*NOBR-by-M*(N+L)  part of this array
          contains details of the complete orthogonal factorization
          of the coefficient matrix  T  of the least squares problem
          which is solved for getting the system matrices B and D.

  LDR     INTEGER
          The leading dimension of the array  R.
          LDR >= MAX( 1, (N+L)*M*NOBR ).

  X       (output) DOUBLE PRECISION array, dimension
          ( (N+L)*M*NOBR )
          The leading  M*(N+L)  elements of this array contain the
          least squares solution of the system  T*X = Kv.
          The remaining elements are used as workspace (to store the
          corresponding part of the vector Kv = vec(K)).

  B       (output) DOUBLE PRECISION array, dimension ( LDB,M )
          The leading N-by-M part of this array contains the system
          input matrix.

  LDB     INTEGER
          The leading dimension of the array B.  LDB >= N.

  D       (output) DOUBLE PRECISION array, dimension ( LDD,M )
          If  JOB = 'D',  the leading L-by-M part of this array
          contains the system input-output matrix.
          If  JOB = 'B',  this array is not referenced.

  LDD     INTEGER
          The leading dimension of the array D.
          LDD >= L, if  JOB = 'D';
          LDD >= 1, if  JOB = 'B'.

  TOL     DOUBLE PRECISION
          The tolerance to be used for estimating the rank of
          matrices. If the user sets  TOL > 0,  then the given value
          of  TOL  is used as a lower bound for the reciprocal
          condition number;  an m-by-n matrix whose estimated
          condition number is less than  1/TOL  is considered to
          be of full rank.  If the user sets  TOL <= 0,  then an
          implicitly computed, default tolerance, defined by
          TOLDEF = m*n*EPS,  is used instead, where  EPS  is the
          relative machine precision (see LAPACK Library routine
          DLAMCH).

  IWORK   INTEGER array, dimension ( M*(N+L) )

  DWORK   DOUBLE PRECISION array, dimension ( LDWORK )
          On exit, if  INFO = 0,  DWORK(1) returns the optimal value
          of LDWORK,  and, if  M > 0,  DWORK(2)  contains the
          reciprocal condition number of the triangular factor of
          the matrix  T.
          On exit, if  INFO = -26,  DWORK(1)  returns the minimum
          value of LDWORK.

  LDWORK  INTEGER
          The length of the array DWORK.
          LDWORK >= MAX( (N+L)*(N+L), 4*M*(N+L)+1 ).
          For good performance,  LDWORK  should be larger.

  IWARN   INTEGER
          = 0:  no warning;
          = 4:  the least squares problem to be solved has a
                rank-deficient coefficient matrix.

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value.

  The matrix  T  is computed, evaluating the sum of Kronecker
  products, and then the linear system  T*X = Kv  is solved in a
  least squares sense. The matrices  B  and  D  are then directly
  obtained from the least squares solution.

  [1] Verhaegen M., and Dewilde, P.
      Subspace Model Identification. Part 1: The output-error
      state-space model identification class of algorithms.
      Int. J. Control, 56, pp. 1187-1210, 1992.

  [2] Van Overschee, P., and De Moor, B.
      N4SID: Two Subspace Algorithms for the Identification
      of Combined Deterministic-Stochastic Systems.
      Automatica, Vol.30, No.1, pp. 75-93, 1994.

  [3] Van Overschee, P.
      Subspace Identification : Theory - Implementation -
      Applications.
      Ph. D. Thesis, Department of Electrical Engineering,
      Katholieke Universiteit Leuven, Belgium, Feb. 1995.

  The implemented method is numerically stable.

  None

Program Text

  None

  None

  None