## TB04BD

### Transfer matrix of a given state-space representation (A,B,C,D), using the pole-zeros method

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

Purpose

```  To compute the transfer function matrix G of a state-space
representation (A,B,C,D) of a linear time-invariant multivariable
system, using the pole-zeros method. Each element of the transfer
function matrix is returned in a cancelled, minimal form, with
numerator and denominator polynomials stored either in increasing
or decreasing order of the powers of the indeterminate.

```
Specification
```      SUBROUTINE TB04BD( JOBD, ORDER, EQUIL, N, M, P, MD, A, LDA, B,
\$                   LDB, C, LDC, D, LDD, IGN, LDIGN, IGD, LDIGD,
\$                   GN, GD, TOL, IWORK, DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
CHARACTER          EQUIL, JOBD, ORDER
DOUBLE PRECISION   TOL
INTEGER            INFO, LDA, LDB, LDC, LDD, LDIGD, LDIGN, LDWORK,
\$                   M, MD, N, P
C     .. Array Arguments ..
DOUBLE PRECISION   A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
\$                   DWORK(*), GD(*), GN(*)
INTEGER            IGD(LDIGD,*), IGN(LDIGN,*), IWORK(*)

```
Arguments

Mode Parameters

```  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 to be a zero matrix.

ORDER   CHARACTER*1
Specifies the order in which the polynomial coefficients
are stored, as follows:
= 'I':  Increasing order of powers of the indeterminate;
= 'D':  Decreasing order of powers of the indeterminate.

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.

```
Input/Output Parameters
```  N       (input) INTEGER
The order of the system (A,B,C,D).  N >= 0.

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

P       (input) INTEGER
The number of the system outputs.  P >= 0.

MD      (input) INTEGER
The maximum degree of the polynomials in G, plus 1. An
upper bound for MD is N+1.  MD >= 1.

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 EQUIL = 'S', the leading N-by-N part of this
array contains the balanced matrix inv(S)*A*S, as returned
by SLICOT Library routine TB01ID.
If EQUIL = 'N', this array is unchanged on exit.

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 input matrix B.
On exit, the contents of B are destroyed: all elements but
those in the first row are set to zero.

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 output matrix C.
On exit, if EQUIL = 'S', the leading P-by-N part of this
array contains the balanced matrix C*S, as returned by
SLICOT Library routine TB01ID.
If EQUIL = 'N', this array is unchanged on exit.

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 matrix D.
If JOBD = 'Z', the array D is not referenced.

LDD     INTEGER
The leading dimension of array D.
LDD >= MAX(1,P), if JOBD = 'D';
LDD >= 1,        if JOBD = 'Z'.

IGN     (output) INTEGER array, dimension (LDIGN,M)
The leading P-by-M part of this array contains the degrees
of the numerator polynomials in the transfer function
matrix G. Specifically, the (i,j) element of IGN contains
the degree of the numerator polynomial of the transfer
function G(i,j) from the j-th input to the i-th output.

LDIGN   INTEGER
The leading dimension of array IGN.  LDIGN >= max(1,P).

IGD     (output) INTEGER array, dimension (LDIGD,M)
The leading P-by-M part of this array contains the degrees
of the denominator polynomials in the transfer function
matrix G. Specifically, the (i,j) element of IGD contains
the degree of the denominator polynomial of the transfer
function G(i,j).

LDIGD   INTEGER
The leading dimension of array IGD.  LDIGD >= max(1,P).

GN      (output) DOUBLE PRECISION array, dimension (P*M*MD)
This array contains the coefficients of the numerator
polynomials, Num(i,j), of the transfer function matrix G.
The polynomials are stored in a column-wise order, i.e.,
Num(1,1), Num(2,1), ..., Num(P,1), Num(1,2), Num(2,2),
..., Num(P,2), ..., Num(1,M), Num(2,M), ..., Num(P,M);
MD memory locations are reserved for each polynomial,
hence, the (i,j) polynomial is stored starting from the
location ((j-1)*P+i-1)*MD+1. The coefficients appear in
increasing or decreasing order of the powers of the
indeterminate, according to ORDER.

GD      (output) DOUBLE PRECISION array, dimension (P*M*MD)
This array contains the coefficients of the denominator
polynomials, Den(i,j), of the transfer function matrix G.
The polynomials are stored in the same way as the
numerator polynomials.

```
Tolerances
```  TOL     DOUBLE PRECISION
The tolerance to be used in determining the
controllability of a single-input system (A,b) or (A',c'),
where b and c' are columns in B and C' (C transposed). If
the user sets TOL > 0, then the given value of TOL is used
as an absolute tolerance; elements with absolute value
less than TOL are considered neglijible. If the user sets
TOL <= 0, then an implicitly computed, default tolerance,
defined by TOLDEF = N*EPS*MAX( NORM(A), NORM(bc) ) is used
instead, where EPS is the machine precision (see LAPACK
Library routine DLAMCH), and bc denotes the currently used
column in B or C' (see METHOD).

```
Workspace
```  IWORK   INTEGER array, dimension (N)

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, N*(N+P) +
MAX( N + MAX( N,P ), N*(2*N+5)))
If N >= P, N >= 1, the formula above can be written as
LDWORK >= N*(3*N + P + 5).
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:  the QR algorithm failed to converge when trying to
compute the zeros of a transfer function;
= 2:  the QR algorithm failed to converge when trying to
compute the poles of a transfer function.
The errors INFO = 1 or 2 are unlikely to appear.

```
Method
```  The routine implements the pole-zero method proposed in .
This method is based on an algorithm for computing the transfer
function of a single-input single-output (SISO) system.
Let (A,b,c,d) be a SISO system. Its transfer function is computed
as follows:

1) Find a controllable realization (Ac,bc,cc) of (A,b,c).
2) Find an observable realization (Ao,bo,co) of (Ac,bc,cc).
3) Compute the r eigenvalues of Ao (the poles of (Ao,bo,co)).
4) Compute the zeros of (Ao,bo,co,d).
5) Compute the gain of (Ao,bo,co,d).

This algorithm can be implemented using only orthogonal
transformations . However, for better efficiency, the
implementation in TB04BD uses one elementary transformation
in Step 4 and r elementary transformations in Step 5 (to reduce
an upper Hessenberg matrix to upper triangular form). These
special elementary transformations are numerically stable
in practice.

In the multi-input multi-output (MIMO) case, the algorithm
computes each element (i,j) of the transfer function matrix G,
for i = 1 : P, and for j = 1 : M. For efficiency reasons, Step 1
is performed once for each value of j (each column of B). The
matrices Ac and Ao result in Hessenberg form.

```
References
```   Varga, A. and Sima, V.
Numerically Stable Algorithm for Transfer Function Matrix
Evaluation.
Int. J. Control, vol. 33, nr. 6, pp. 1123-1133, 1981.

```
Numerical Aspects
```  The algorithm is numerically stable in practice and requires about
20*N**3 floating point operations at most, but usually much less.

```
```  For maximum efficiency of index calculations, GN and GD are
implemented as one-dimensional arrays.

```
Example

Program Text

```*     TB04BD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
INTEGER          NIN, NOUT
PARAMETER        ( NIN = 5, NOUT = 6 )
INTEGER          NMAX, MMAX, PMAX, MDMAX
PARAMETER        ( NMAX = 20, MMAX = 20, PMAX = 20,
\$                   MDMAX = NMAX + 1 )
INTEGER          PMNMAX
PARAMETER        ( PMNMAX = PMAX*MMAX*MDMAX )
INTEGER          LDA, LDB, LDC, LDD, LDIGD, LDIGN
PARAMETER        ( LDA = NMAX, LDB = NMAX, LDC = PMAX,
\$                   LDD = PMAX, LDIGD = PMAX, LDIGN = PMAX )
INTEGER          LIWORK
PARAMETER        ( LIWORK = NMAX )
INTEGER          LDWORK
PARAMETER        ( LDWORK = NMAX*( NMAX + PMAX ) +
\$                            MAX( NMAX + MAX( NMAX, PMAX ),
\$                                 NMAX*( 2*NMAX + 5 ) ) )
*     .. Local Scalars ..
DOUBLE PRECISION TOL
INTEGER          I, IJ, INFO, J, K, M, MD, N, P
CHARACTER*1      JOBD, ORDER, EQUIL
CHARACTER*132    ULINE
*     .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
\$                 D(LDD,MMAX), DWORK(LDWORK), GD(PMNMAX),
\$                 GN(PMNMAX)
INTEGER          IGD(LDIGD,MMAX), IGN(LDIGN,MMAX), IWORK(LIWORK)
*     .. External Functions ..
LOGICAL          LSAME
EXTERNAL         LSAME
*     .. External Subroutines ..
EXTERNAL         TB04BD
*     .. 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, TOL, JOBD, ORDER, EQUIL
MD = N + 1
ULINE(1:20) = ' '
DO 20 I = 21, 132
ULINE(I:I) = '-'
20 CONTINUE
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99991 ) 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 = 99990 ) M
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), I = 1,N ), J = 1,M )
IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
WRITE ( NOUT, FMT = 99989 ) 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 the transfer matrix T(s) of (A,B,C,D).
CALL TB04BD( JOBD, ORDER, EQUIL, N, M, P, MD, A, LDA, B,
\$                      LDB, C, LDC, D, LDD, IGN, LDIGN, IGD, LDIGD,
\$                      GN, GD, TOL, IWORK, DWORK, LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
IF ( LSAME( ORDER, 'I' ) ) THEN
WRITE ( NOUT, FMT = 99997 )
ELSE
WRITE ( NOUT, FMT = 99996 )
END IF
WRITE ( NOUT, FMT = 99995 )
DO 60 J = 1, M
DO 40 I = 1, P
IJ = ( (J-1)*P + I-1 )*MD + 1
WRITE ( NOUT, FMT = 99994 ) I, J,
\$                    ( GN(K), K = IJ,IJ+IGN(I,J) )
WRITE ( NOUT, FMT = 99993 )
\$                          ULINE(1:7*(IGD(I,J)+1)+21)
WRITE ( NOUT, FMT = 99992 )
\$                        ( GD(K), K = IJ,IJ+IGD(I,J) )
40                CONTINUE
60             CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TB04BD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TB04BD = ',I2)
99997 FORMAT (/' The polynomial coefficients appear in increasing',
\$         ' order'/' of the powers of the indeterminate')
99996 FORMAT (/' The polynomial coefficients appear in decreasing',
\$         ' order'/' of the powers of the indeterminate')
99995 FORMAT (/' The coefficients of polynomials in the transfer matri',
\$       'x T(s) are ')
99994 FORMAT (/' element (',I2,',',I2,') is ',20(1X,F6.2))
99993 FORMAT (1X,A)
99992 FORMAT (20X,20(1X,F6.2))
99991 FORMAT (/' N is out of range.',/' N = ',I5)
99990 FORMAT (/' M is out of range.',/' M = ',I5)
99989 FORMAT (/' P is out of range.',/' P = ',I5)
END
```
Program Data
``` TB04BD EXAMPLE PROGRAM DATA
3     2     2  0.0         D     I     N
-1.0   0.0   0.0
0.0  -2.0   0.0
0.0   0.0  -3.0
0.0   1.0  -1.0
1.0   1.0   0.0
0.0   1.0   1.0
1.0   1.0   1.0
1.0   0.0
0.0   1.0
```
Program Results
``` TB04BD EXAMPLE PROGRAM RESULTS

The polynomial coefficients appear in increasing order
of the powers of the indeterminate

The coefficients of polynomials in the transfer matrix T(s) are

element ( 1, 1) is    7.00   5.00   1.00
----------------------
6.00   5.00   1.00

element ( 2, 1) is    1.00
----------------------
6.00   5.00   1.00

element ( 1, 2) is    1.00
---------------
2.00   1.00

element ( 2, 2) is    5.00   5.00   1.00
----------------------
2.00   3.00   1.00
```