#

# analysis.py

#

# Copyright 2010 Enrico Avventi <avventi@Lonewolf>

#

# This program is free software; you can redistribute it and/or modify

# it under the terms of the GNU General Public License version 2 as

# published by the Free Software Foundation.

#

# This program is distributed in the hope that it will be useful,

# but WITHOUT ANY WARRANTY; without even the implied warranty of

# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

# GNU General Public License for more details.

#

# You should have received a copy of the GNU General Public License

# along with this program; if not, write to the Free Software

# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,

# MA 02110-1301, USA.

import numpy as np

from . import _wrapper

from .exceptions import raise_if_slycot_error, SlycotParameterError

def ab01nd(n, m, A, B, jobz='N', tol=0, ldwork=None):

""" Ac,Bc,ncont,indcon,nblk,Z,tau = ab01nd_i(n,m,A,B,[jobz,tol,ldwork])

To find a controllable realization for the linear time-invariant

multi-input system

::

dX/dt = A * X + B * U,

where A and B are N-by-N and N-by-M matrices, respectively,

which are reduced by this routine to orthogonal canonical form

using (and optionally accumulating) orthogonal similarity

transformations. Specifically, the pair ``(A, B)`` is reduced to

the pair ``(Ac, Bc), Ac = Z' * A * Z, Bc = Z' * B``, given by

::

[ Acont * ] [ Bcont ]

Ac = [ ], Bc = [ ],

[ 0 Auncont ] [ 0 ]

and

[ A11 A12 . . . A1,p-1 A1p ] [ B1 ]

[ A21 A22 . . . A2,p-1 A2p ] [ 0 ]

[ 0 A32 . . . A3,p-1 A3p ] [ 0 ]

Acont = [ . . . . . . . ], Bc = [ . ],

[ . . . . . . ] [ . ]

[ . . . . . ] [ . ]

[ 0 0 . . . Ap,p-1 App ] [ 0 ]

where the blocks ``B1, A21, ..., Ap,p-1`` have full row ranks and

`p` is the controllability index of the pair. The size of the

block `Auncont` is equal to the dimension of the uncontrollable

subspace of the pair ``(A, B)``.

Parameters

----------

n : int

The order of the original state-space representation, i.e.

the order of the matrix A. ``n > 0``.

m : int

The number of system inputs, or of columns of B. ``m > 0``.

A : (n, n) array_like

The original state dynamics matrix A.

B : (n, m) array_like

The input matrix B.

jobz : {'N', 'F', 'I'}, optional

Indicates whether the user wishes to accumulate in a matrix Z

the orthogonal similarity transformations for reducing the system,

as follows:

:= 'N': Do not form Z and do not store the orthogonal transformations;

(default)

:= 'F': Do not form Z, but store the orthogonal transformations in

the factored form;

:= 'I': Z is initialized to the unit matrix and the orthogonal

transformation matrix Z is returned.

tol : float, optional

The tolerance to be used in rank determination when transforming

``(A, B)``. If ``tol <= 0`` a default value is used.

ldwork : int, optional

The length of the cache array. ``ldwork >= max(n, 3*m)``.

For optimum performance it should be larger.

default: ``ldwork = max(n, 3*m)``

Returns

-------

Ac : (n, n) ndarray

The leading ncont-by-ncont part contains the upper block

Hessenberg state dynamics matrix Acont in Ac, given by Z'*A*Z,

of a controllable realization for the original system. The

elements below the first block-subdiagonal are set to zero.

Bc : (n, m) ndarray

The leading ncont-by-m part of this array contains the transformed

input matrix Bcont in Bc, given by ``Z'*B``, with all elements but the

first block set to zero.

ncont : int

The order of the controllable state-space representation.

indcon : int

The controllability index of the controllable part of the system

representation.

nblk : (n, ) int ndarray

The leading indcon elements of this array contain the the orders of

the diagonal blocks of Acont.

Z : (n, n) ndarray

- If jobz = 'I', then the leading N-by-N part of this array contains

the matrix of accumulated orthogonal similarity transformations

which reduces the given system to orthogonal canonical form.

- If jobz = 'F', the elements below the diagonal, with the array tau,

represent the orthogonal transformation matrix as a product of

elementary reflectors. The transformation matrix can then be

obtained by calling the LAPACK Library routine DORGQR.

- If jobz = 'N', the array Z is `None`.

tau : (n, ) ndarray

The elements of tau contain the scalar factors of the

elementary reflectors used in the reduction of B and A.

"""

hidden = ' (hidden by the wrapper)'

arg_list = ['jobz', 'n', 'm', 'A', 'LDA'+hidden, 'B', 'LDB'+hidden,

'ncont', 'indcon', 'nblk', 'Z', 'LDZ'+hidden, 'tau', 'tol',

'IWORK'+hidden, 'DWORK'+hidden, 'ldwork', 'info'+hidden]

if ldwork is None:

ldwork = max(n, 3*m)

Ac, Bc, ncont, indcon, nblk, Z, tau, info = _wrapper.ab01nd(

jobz, n, m, A, B, tol=tol, ldwork=ldwork)

raise_if_slycot_error(info, arg_list)

if jobz == "N":

Z = None

return Ac, Bc, ncont, indcon, nblk, Z, tau

def ab04md(type_t, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None):

""" At,Bt,Ct,Dt = ab04md(type_t, n, m, p, A, B, C, D, [alpha, beta,ldwork])

Parameters

----------

type_t : {'D','C'}

Indicates the type of the original system and the

transformation to be performed as follows:

= 'D': discrete-time -> continuous-time;

= 'C': continuous-time -> discrete-time.

n : int

The order of the matrix A, the number of rows of matrix B and

the number of columns of matrix C. It represents the dimension of

the state vector. n > 0.

m : int

The number of columns of matrix B. It represents the dimension of

the input vector. m > 0.

p : int

The number of rows of matrix C. It represents the dimension of

the output vector. p > 0.

A : (n, n) array_like

The leading n-by-n part of this array must contain the system state

matrix A.

B : (n, m) array_like

The leading n-by-m part of this array must contain the system input

matrix B.

C : (p, n) array_like

The leading p-by-n part of this array must contain the system output

matrix C.

D : (p, m) array_like

The leading p-by-m part of this array must contain the system direct

transmission matrix D.

alpha : float, optional

Parameter specifying the bilinear transformation.

Recommended values for stable systems: alpha = 1, alpha != 0,

Default is 1.0.

beta : float, optional

Parameter specifying the bilinear transformation.

Recommended values for stable systems: beta = 1, beta != 0,

Default is 1.0.

ldwork : int, optional

The length of the cache array.

ldwork >= max(1, n), default is max(1, n)

Returns

-------

At : (n, n) ndarray

The state matrix At of the transformed system.

Bt : (n, m) ndarray

The input matrix Bt of the transformed system.

Ct : (p, n) ndarray

The output matrix Ct of the transformed system.

Dt : (p, m) ndarray

The transmission matrix Dt of the transformed system.

Raises

------

SlycotArithmeticError

:info == 1:

If the matrix (ALPHA*I + A) is exactly singular

:info == 2:

If the matrix (BETA*I - A) is exactly singular.

"""

hidden = ' (hidden by the wrapper)'

arg_list = ['type_t', 'n', 'm', 'p', 'alpha', 'beta',

'A', 'LDA'+hidden, 'B', 'LDB'+hidden, 'C', 'LDC'+hidden, 'D', 'LDD'+hidden,

'IWORK'+hidden, 'DWORK'+hidden, 'ldwork', 'info'+hidden]

if ldwork is None:

ldwork = max(1, n)

out = _wrapper.ab04md(type_t, n, m, p, alpha, beta, A, B, C, D, ldwork=ldwork)

info=out[-1]

raise_if_slycot_error(info, arg_list)

return out[:-1]

def ab05md(n1,m1,p1,n2,p2,A1,B1,C1,D1,A2,B2,C2,D2,uplo='U'):

""" n,a,b,c,d = ab05md(n1,m1,p1,n2,p2,a1,b1,c1,d1,a2,b2,c2,d2,[uplo])

To obtain the state-space model (A,B,C,D) for the cascaded

inter-connection of two systems, each given in state-space form.

Parameters

----------

n1 : int

The number of state variables in the first system, i.e. the order

of the matrix A1. n1 > 0.

m1 : int

The number of input variables for the first system. m1 > 0.

p1 : int

The number of output variables from the first system and the number

of input variables for the second system. p1 > 0.

n2 : int

The number of state variables in the second system, i.e. the order

of the matrix A2. n2 > 0.

p2 : int

The number of output variables from the second system. p2 > 0.

A1 : (n1, n1) array_like

The leading n1-by-n1 part of this array must contain the state

transition matrix A1 for the first system.

B1 : (n1, m1) array_like

The leading n1-by-m1 part of this array must contain the input/state

matrix B1 for the first system.

C1 : (p1, n1) array_like

The leading p1-by-n1 part of this array must contain the state/output

matrix C1 for the first system.

D1 : (p1, m1) array_like

The leading p1-by-m1 part of this array must contain the input/output

matrix D1 for the first system.

A2 : (n2, n2) array_like

The leading n2-by-n2 part of this array must contain the state

transition matrix A2 for the second system.

B2 : (n2, p1) array_like

The leading n2-by-p1 part of this array must contain the input/state

matrix B2 for the second system.

C2 : (p2, n2) array_like

The leading p2-by-n2 part of this array must contain the state/output

matrix C2 for the second system.

D2 : (p2, p1) array_like

The leading p2-by-p1 part of this array must contain the input/output

matrix D2 for the second system.

uplo : {'U', 'L'}, optional

Indicates whether the user wishes to obtain the matrix A in

the upper or lower block diagonal form, as follows:

= 'U': Obtain A in the upper block diagonal form;

= 'L': Obtain A in the lower block diagonal form.

Default is `U`.

Returns

-------

n : int

The number of state variables (n1 + n2) in the resulting system,

i.e. the order of the matrix A, the number of rows of B and

the number of columns of C.

A : (n1+n2, n1+n2) ndarray

The leading N-by-N part of this array contains the state transition

matrix A for the cascaded system.

B : (n1+n2, m1) ndarray

The leading n-by-m1 part of this array contains the input/state

matrix B for the cascaded system.

C : (p2, n1+n2) ndarray

The leading p2-by-n part of this array contains the state/output

matrix C for the cascaded system.

D : (p2, m1) ndarray

The leading p2-by-m1 part of this array contains the input/output

matrix D for the cascaded system.

Notes

-----

The implemented methods rely on accuracy enhancing square-root or

balancing-free square-root techniques.

The algorithms require less than 30N^3 floating point operations.

"""

hidden = ' (hidden by the wrapper)'

arg_list = ['uplo', 'OVER'+hidden, 'n1', 'm1', 'p1', 'n2', 'p2', 'A1',

'LDA1'+hidden, 'B1', 'LDB1'+hidden, 'C1', 'LDC1'+hidden, 'D1',

'LDD1'+hidden, 'A2', 'LDA2'+hidden, 'B2', 'LDB2'+hidden, 'C2',

'LDC2'+hidden, 'D2', 'LDD2'+hidden, 'n', 'A', 'LDA'+hidden, 'B',

'LDB'+hidden, 'C', 'LDC'+hidden, 'D', 'LDD'+hidden, 'DWORK'+hidden,

'ldwork', 'info'+hidden ]

out = _wrapper.ab05md(n1,m1,p1,n2,p2,A1,B1,C1,D1,A2,B2,C2,D2,uplo=uplo)

raise_if_slycot_error(out[-1], arg_list)

return out[:-1]

def ab05nd(n1,m1,p1,n2,A1,B1,C1,D1,A2,B2,C2,D2,alpha=1.0,ldwork=None):

""" n,A,B,C,D = ab05nd(n1,m1,p1,n2,A1,B1,C1,D1,A2,B2,C2,D2,[alpha,ldwork])

To obtain the state-space model (A,B,C,D) for the feedback inter-connection

of two systems, each given in state-space form.

Parameters

----------

n1 : int

The number of state variables in the first system, i.e. the order

of the matrix A1. n1 > 0.

m1 : int

The number of input variables for the first system and the number

of output variables from the second system. m1 > 0.

p1 : int

The number of output variables from the first system and the number

of input variables for the second system. p1 > 0.

n2 : int

The number of state variables in the second system, i.e. the order

of the matrix A2. n2 > 0.

A1 : (n1, n1) array_like

The leading n1-by-n1 part of this array must contain the state

transition matrix A1 for the first system.

B1 : (n1, m1) array_like

The leading n1-by-m1 part of this array must contain the input/state

matrix B1 for the first system.

C1 : (p1, n1) array_like

The leading p1-by-n1 part of this array must contain the state/output

matrix C1 for the first system.

D1 : (p1, m1) array_like

The leading p1-by-m1 part of this array must contain the input/output

matrix D1 for the first system.

A2 : (n2, n2) array_like

The leading n2-by-n2 part of this array must contain the state

transition matrix A2 for the second system.

B2 : (n2, p1) array_like

The leading n2-by-p1 part of this array must contain the input/state

matrix B2 for the second system.

C2 : (m1, n2) array_like

The leading m1-by-n2 part of this array must contain the state/output

matrix C2 for the second system.

D2 : (m1, p1) array_like

The leading m1-by-p1 part of this array must contain the input/output

matrix D2 for the second system.

alpha : float, optional

A coefficient multiplying the transfer-function matrix (or the

output equation) of the second system. i.e alpha = +1 corresponds

to positive feedback, and alpha = -1 corresponds to negative

feedback.

Default is `1.0`.

ldwork : int, optional

The length of the cache array. ldwork >= max(p1*p1,m1*m1,n1*p1).

Default is max(p1*p1,m1*m1,n1*p1).

Returns

-------

n : int

The number of state variables (n1 + n2) in the connected system, i.e.

the order of the matrix A, the number of rows of B and the number of

columns of C.

A : (n1+n2, n1+n2) ndarray

The leading n-by-n part of this array contains the state transition

matrix A for the connected system.

B : (n1+n2, m1) ndarray

The leading n-by-m1 part of this array contains the input/state

matrix B for the connected system.

C : (p1, n1, n2) ndarray

The leading p1-by-n part of this array contains the state/output

matrix C for the connected system.

D : (p1, m1) ndarray

The leading p1-by-m1 part of this array contains the input/output

matrix D for the connected system.

Raises

------

SlycotArithmeticError

:1 <= info <= p1:

the system is not completely controllable. That is, the matrix

``(I + ALPHA*D1*D2)`` is exactly singular (the element

``U(i,i)```` of the upper triangular factor of ``LU```

factorization is exactly zero), possibly due to

rounding errors.

"""

hidden = ' (hidden by the wrapper)'

arg_list = ['over'+hidden, 'n1', 'm1', 'p1', 'n2', 'alpha', 'A1', 'LDA1'+hidden,

'B1', 'LDB1'+hidden, 'C1', 'LDC1'+hidden, 'D1', 'LDD1'+hidden, 'A2',

'LDA2'+hidden, 'B2', 'LDB2'+hidden, 'C2', 'LDC2'+hidden, 'D2',

'LDD2'+hidden, 'n', 'A', 'LDA'+hidden, 'B', 'LDB'+hidden, 'C',

'LDC'+hidden, 'D', 'LDD'+hidden, 'IWORK'+hidden, 'DWORK'+hidden,

'ldwork', 'info'+hidden]

if ldwork is None:

ldwork = max(p1*p1,m1*m1,n1*p1)

out = _wrapper.ab05nd(n1,m1,p1,n2,alpha,A1,B1,C1,D1,A2,B2,C2,D2,ldwork=ldwork)

raise_if_slycot_error(out[-1], arg_list, ab05nd, locals())

return out[:-1]

def ab07nd(n,m,A,B,C,D,ldwork=None):

""" Ai,Bi,Ci,Di,rcond = ab07nd(n,m,A,B,C,D,[ldwork])

To compute the inverse (Ai,Bi,Ci,Di) of a given system (A,B,C,D).

Parameters

----------

n : int

The order of the state matrix A. n >= 0.

m : int

The number of system inputs and outputs. m >= 0.

A : (n, n) array_like

The leading n-by-n part of this array must contain the state matrix

A of the original system.

B : (n, m) array_like

The leading n-by-m part of this array must contain the input matrix

B of the original system.

C : (m, n) array_like

The leading m-by-n part of this array must contain the output matrix

C of the original system.

D : (m, m) array_like

The leading m-by-m part of this array must contain the feedthrough

matrix D of the original system.

ldwork : int, optional

The length of the cache array. The default value is max(1,4*m),

for better performance should be larger.

Returns

-------

Ai : (n, n) ndarray

The leading n-by-n part of this array contains the state matrix Ai

of the inverse system.

Bi : (n, m) ndarray

The leading n-by-m part of this array contains the input matrix Bi

of the inverse system.

Ci : (m, n) ndarray

The leading m-by-n part of this array contains the output matrix Ci

of the inverse system.

Di : (m, m) ndarray

The leading m-by-m part of this array contains the feedthrough

matrix Di of the inverse system.

rcond : float

The estimated reciprocal condition number of the feedthrough matrix

D of the original system.

Warns

-----

SlycotResultWarning

:1 <= info <= m:

the matrix `D` is exactly singular; the ({info},{info})

diagonal element is zero, `RCOND` was set to zero;

:info == m+1:

the matrix `D` is numerically singular, i.e., `RCOND`

is less than the relative machine precision, `EPS`

(see LAPACK Library routine DLAMCH). The

calculations have been completed, but the results

could be very inaccurate.

"""

hidden = ' (hidden by the wrapper)'

arg_list = ['n', 'm', 'A', 'LDA' + hidden, 'B', 'LDB' + hidden,

'C', 'LDC' + hidden, 'D', 'LDD' + hidden, 'rcond',

'IWORK' + hidden, 'DWORK' + hidden, 'ldwork', 'INFO' + hidden]

if ldwork is None:

ldwork = max(1, 4*m)

out = _wrapper.ab07nd(n, m, A, B, C, D, ldwork=ldwork)

raise_if_slycot_error(out[-1], arg_list, ab07nd.__doc__, locals())

return out[:-1]

def ab08nd(n,m,p,A,B,C,D,equil='N',tol=0,ldwork=None):

""" nu,rank,dinfz,nkror,nkrol,infz,kronr,kronl,Af,Bf = ab08nd(n,m,p,A,B,C,D,[equil,tol,ldwork])

To construct for a linear multivariable system described by a state-space

model (A,B,C,D) a regular pencil (Af - lambda*Bf ) which has the invariant

zeros of the system as generalized eigenvalues.

The routine also computes the orders of the infinite zeros and the

right and left Kronecker indices of the system (A,B,C,D).

Parameters

----------

n : int

The number of state variables. n >= 0.

m : int

The number of system inputs. m >= 0.

p : int

The number of system outputs. p >= 0.

A : (n, n) array_like

The leading n-by-n part of this array must contain the state

dynamics matrix A of the system.

B : (n, m) array_like

The leading n-by-m part of this array must contain the input/state

matrix B of the system.

C : (p, n) array_like

The leading p-by-n part of this array must contain the state/output

matrix C of the system.

D : (p, m) array_like

The leading p-by-m part of this array must contain the direct

transmission matrix D of the system.

equil : {'S', 'N'}, optional

Specifies whether the user wishes to balance the compound matrix

as follows:

= 'S': Perform balancing (scaling);

= 'N': Do not perform balancing.

Default is `N`.

tol : float, optional

A tolerance used in rank decisions to determine the effective rank,

which is defined as the order of the largest leading (or trailing)

triangular submatrix in the QR (or RQ) factorization with column

(or row) pivoting whose estimated condition number is less than 1/tol.

Default is `0.0`.

ldwork : int, optional

The length of the cache array. The default value is n + 3*max(m,p),

for better performance should be larger.

Default is None.

Returns

-------

nu : int

The number of (finite) invariant zeros.

rank : int

The normal rank of the transfer function matrix.

dinfz : int

The maximum degree of infinite elementary divisors.

nkror : int

The number of right Kronecker indices.

nkrol : int

The number of left Kronecker indices.

infz : (n, ) ndarray

The leading dinfz elements of infz contain information on the

infinite elementary divisors as follows: the system has infz(i)

infinite elementary divisors of degree i, where i = 1,2,...,dinfz.

kronr :(max(n,m)+1, ) ndarray

the leading nkror elements of this array contain the right kronecker

(column) indices.

kronl : (max(n,p)+1, ) ndarray

the leading nkrol elements of this array contain the left kronecker

(row) indices.

Af : (max(1,n+m), n+min(p,m)) ndarray

the leading nu-by-nu part of this array contains the coefficient

matrix Af of the reduced pencil. the remainder of the leading

(n+m)-by-(n+min(p,m)) part is used as internal workspace.

Bf : (max(1,n+p), n+m) ndarray

The leading nu-by-nu part of this array contains the coefficient

matrix Bf of the reduced pencil. the remainder of the leading

(n+p)-by-(n+m) part is used as internal workspace.

"""

hidden = ' (hidden by the wrapper)'

arg_list = ['equil', 'n', 'm', 'p', 'A', 'LDA'+hidden, 'B', 'LDB'+hidden,

'C', 'LDC'+hidden, 'D', 'LDD'+hidden, 'nu', 'rank', 'dinfz', 'nkror',

'nkrol', 'infz', 'kronr', 'kronl', 'Af', 'LDAF'+hidden, 'Bf',

'LDBF'+hidden, 'tol', 'IWORK'+hidden, 'DWORK'+hidden, 'ldwork',

'INFO'+hidden]

if ldwork is None:

ldwork = n+3*max(m,p) #only an upper bound

out = _wrapper.ab08nd(n,m,p,A,B,C,D,equil=equil,tol=tol,ldwork=ldwork)

raise_if_slycot_error(out[-1], arg_list)

return out[:-1]

def ab08nz(n, m, p, A, B, C, D, equil='N', tol=0., lzwork=None):

""" nu,rank,dinfz,nkror,nkrol,infz,kronr,kronl,Af,Bf = ab08nz(n,m,p,A,B,C,D,[equil,tol,lzwork])

To construct for a linear multivariable system described by a state-space

model (A,B,C,D) a regular pencil (Af - lambda*Bf ) which has the invariant

zeros of the system as generalized eigenvalues.

The routine also computes the orders of the infinite zeros and the

right and left Kronecker indices of the system (A,B,C,D).

Parameters

----------

n : int

The number of state variables. n >= 0.

m : int

The number of system inputs. m >= 0.

p : int

The number of system outputs. p >= 0.

A : (n, n) array_like

The leading n-by-n part of this array must contain the state

dynamics matrix A of the system.

B : (n, m) array_like

The leading n-by-m part of this array must contain the input/state

matrix B of the system.

C : (p, n) array_like

The leading p-by-n part of this array must contain the state/output

matrix C of the system.

D : (p, m) array_like

The leading p-by-m part of this array must contain the direct

transmission matrix D of the system.

equil : {'S', 'N'}, optional

Specifies whether the user wishes to balance the compound matrix

as follows:

= 'S': Perform balancing (scaling);

= 'N': Do not perform balancing.

Default is `N`.

tol : float, optional

A tolerance used in rank decisions to determine the effective rank,

which is defined as the order of the largest leading (or trailing)

triangular submatrix in the QR (or RQ) factorization with column

(or row) pivoting whose estimated condition number is less than 1/tol.

If tol is set to less than SQRT((N+P)*(N+M))*EPS

then the tolerance is taken as SQRT((N+P)*(N+M))*EPS,

where EPS is the machine precision (see LAPACK Library

Routine DLAMCH).

Default is 0.0.

lzwork : int, optional

The length of the internal cache array ZWORK. The default value is

calculated to

MAX( 1,

MIN(P,M) + MAX(3*M-1,N),

MIN(P,N) + MAX(3*P-1,N+P,N+M),

MIN(M,N) + MAX(3*M-1,N+M) )

For optimum performance lzwork should be larger.

If lzwork = -1, then a workspace query is assumed;

the routine only calculates the optimal size of the

ZWORK array, and returns this value in lzwork_opt

Default is None.

Returns

-------

nu : int

The number of (finite) invariant zeros.

rank : int

The normal rank of the transfer function matrix.

dinfz : int

The maximum degree of infinite elementary divisors.

nkror : int

The number of right Kronecker indices.

nkrol : int

The number of left Kronecker indices.

infz : (n, ) ndarray

The leading dinfz elements of infz contain information on the

infinite elementary divisors as follows: the system has infz(i)

infinite elementary divisors of degree i, where i = 1,2,...,dinfz.

kronr : (max(n,m)+1, ) ndarray

the leading nkror elements of this array contain the right kronecker

(column) indices.

kronl : (max(n,p)+1, ) ndarray

the leading nkrol elements of this array contain the left kronecker

(row) indices.

Af : (max(1,n+m), n+min(p,m)) ndarray

the leading nu-by-nu part of this array contains the coefficient

matrix Af of the reduced pencil. the remainder of the leading

(n+m)-by-(n+min(p,m)) part is used as internal workspace.

Bf : (max(1,n+p), n+m) ndarray

The leading nu-by-nu part of this array contains the coefficient

matrix Bf of the reduced pencil. the remainder of the leading

(n+p)-by-(n+m) part is used as internal workspace.

lzwork_opt : int

The optimal value of lzwork.

"""

hidden = ' (hidden by the wrapper)'

arg_list = ['equil', 'n', 'm', 'p',

'a', 'lda' + hidden, 'b', 'ldb' + hidden,

'c', 'ldc' + hidden, 'd', 'ldd' + hidden,

'nu', 'rank', 'dinfz', 'nkror', 'nkrol', 'infz', 'kronr',

'kronl', 'af', 'ldaf' + hidden, 'bf', 'ldbf' + hidden,

'tol', 'iwork' + hidden, 'dwork' + hidden, 'zwork',

'lzwork', 'info']

if lzwork is None:

lzwork = max(min(p, m) + max(3*m-1, n),

min(p, n) + max(3*p-1, n+p, n+m),

min(m, n) + max(3*m-1, n+m))

nu, rank, dinfz, nkror, nkrol, infz, kronr, kronl, Af, Bf, zwork, info \

= _wrapper.ab08nz(n, m, p, A, B, C, D,

equil=equil, tol=tol, lzwork=lzwork)

raise_if_slycot_error(info, arg_list)

return (nu, rank, dinfz, nkror, nkrol, infz, kronr, kronl, Af, Bf,

int(zwork[0].real))

def ab09ad(dico,job,equil,n,m,p,A,B,C,nr=None,tol=0,ldwork=None):

""" nr,Ar,Br,Cr,hsv = ab09ad(dico,job,equil,n,m,p,A,B,C,[nr,tol,ldwork])

Compute reduced order State-Space-Model (Ar, Br, Cr) for a stable system

(A, B, C) by using either the square-root or the balancing-free square-

root Balance & truncate (B & T) model reduction method.

Parameters

----------

dico : {'D', 'C'}

Indicate whether the system is discrete `D` or continuous `C`

job : {'B', 'N'}

Balance `B` or not `N`

equil : {'S', 'N'}

Scale `S` or not `N`

n : int

The number of state variables. n >= 0.

m : int

The number of system inputs. m >= 0.

p : int

The number of system outputs. p >= 0.

A : (n, n) array_like

The leading n-by-n part of this array must contain the state

dynamics matrix A of the system.

B : (n, m) array_like

The leading n-by-m part of this array must contain the input/state

matrix B of the system.

C : (p, n) array_like

The leading p-by-n part of this array must contain the

state/output matrix C of the system.

nr : int, optional

`nr` is the desired order of the resulting reduced order

system. ``0 <= nr <= n``. Automatically determined by `tol` if

``nr is None`` and returned. See return object `nr`.

Default is None.

tol : float, optional

If ``nr is None``, `tol`contains the tolerance for determining the

order of the reduced system. For model reduction, th recommended

value is ``tol = c * HNORM(A, B, C)``, where `c` is a constan in the

interval ``[0.00001, 0.001]`` and ``HNORM(A, B, C)`` is the

Hankel-Norm of the given sysstem (computed in ``HSV(1)``). For

computing a minimal realization, the recommended value is

``tol = n * eps * HNORM(A, B, C)``, where `eps` is the machine

precision (see LAPACK Library Routine `DLAMCH`). This value is

used by default if ``tol <= 0`` on entry. If `nr` is specified,

the value of `tol` is ignored. Default is `0.0`.

ldwork : int, optional

The length of the cache array. The default value is

``n*(2*n+max(n,m,p)+5) + n*(n+1)/2 ~= 3.5*n**2 + 5*n``,

a larger value should lead to better performance.

Default is None.

Returns

-------

nr : int

`nr` is the order of the resulting reduced order model.

`nr` is set as follows:

If on input ``nr is not None``, `nr` is equal to ``MIN(nr,NMIN)``,

where `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 on input ``nr is None``, `nr` is equal to the number of Hankel

singular values greater than ``MAX(tol,n*eps*HNORM(A,B,C))``.

Ar : (nr, nr) ndarray

This array contains the state dynamics matrix `Ar` of the reduced

order system.

Br : (nr, m) ndarray

This array contains the input/state matrix `Br` of the reduced

order system.

Cr : (p, nr) ndarray

This array contains the state/output matrix `Cr` of the reduced

order system.

hsv : (n, ) ndarray

If ``INFO = 0``, it contains the Hankel singular values of

the original system ordered decreasingly. ``HSV(1)`` is the

Hankel norm of the system.

Raises

------

SlycotArithmeticError

:info == 1:

The reduction of A to the real Schur form failed

:info == 2 and dico == 'C':

The state matrix A is not stable

:info == 2 and dico == 'D':

The state matrix A is not convergent

:info == 3:

The computation of Hankel singular values failed

Warns

-----

SlycotResultWarning

:iwarn == 1:

The selected order {nr} is greater

than the order of a minimal realization of the

given system. `nr` was set automatically to {Nr}

corresponding to the order of a minimal realization

of the system

"""

hidden = ' (hidden by the wrapper)'

arg_list = ['dico', 'job', 'equil', 'ordsel', 'n', 'm', 'p', 'nr',

'A', 'lda' + hidden, 'B', 'ldb' + hidden, 'C', 'ldc' + hidden,

'hsv', 'tol', 'iwork' + hidden, 'dwork ' + hidden, 'ldwork',

'iwarn', 'info']

if ldwork is None:

ldwork = max(1, n*(2*n+max(n, max(m, p))+5)+n*(n+1)/2)

if nr is None:

ordsel = 'A'

nr = 0 # order will be computed by the routine

else:

ordsel = 'F'

out = _wrapper.ab09ad(dico, job, equil, ordsel,

n, m, p, nr, A, B, C, tol, ldwork)

Nr, A, B, C, hsv = out[:-2]

raise_if_slycot_error(out[-2:], arg_list, ab09ad.__doc__, locals())

return Nr, A[:Nr, :Nr], B[:Nr, :], C[:, :Nr], hsv

def ab09ax(dico,job,n,m,p,A,B,C,nr=None,tol=0.0,ldwork=None):

"""``nr,Ar,Br,Cr,hsv,T,Ti = ab09ad(dico,job,equil,n,m,p,nr,A,B,C,[nr,tol,ldwork])``

To compute a reduced order model ``(Ar,Br,Cr)`` for a stable original

state-space representation ``(A,B,C)`` by using either the square-root

or the balancing-free square-root Balance & Truncate model

reduction method. The state dynamics matrix `A` of the original

system is an upper quasi-triangular matrix in *real Schur canonical

form.* The matrices of the reduced order system are computed using

the truncation formulas:

``Ar = TI * A * T , Br = TI * B , Cr = C * T`` .

Parameters

----------

dico : {'D', 'C'}

Indicate whether the system is discrete `D` or continuous `C`

job : {'B', 'N'}

Balance `B` or not `N`

n : int

The number of state variables. n >= 0.

m : int

The number of system inputs. m >= 0.

p : int

The number of system outputs. p >= 0.

A : (n, n) array_like

The leading n-by-n part of this array must contain the state

dynamics matrix A of the system *in real Schur form.*

B : (n, m) array_like

The leading n-by-m part of this array must contain the input/state

matrix B of the system.

C : (p, n) array_like

The leading p-by-n part of this array must contain the

state/output matrix C of the system.

nr : int, optional

`nr` is the desired order of the resulting reduced order

system. ``0 <= nr <= n``. Automatically determined by `tol` if

``nr is None`` and returned. See return object `nr`.

Default is None.

tol : float, optional

If ``nr is None``, `tol`contains the tolerance for determining the

order of the reduced system. For model reduction, the recommended

value is ``tol = 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 sysstem (computed in ``HSV(1)``). For

computing a minimal realization, the recommended value is

``tol = n * eps * HNORM(A, B, C)``, where `eps` is the machine

precision (see LAPACK Library Routine `DLAMCH`). This value is

used by default if ``tol <= 0`` on entry. If `nr` is specified,

the value of `tol` is ignored. Default is `0.0`.

ldwork : int, optional

The length of the cache array. The default value is

``n*(2*n+max(n,m,p)+5) + n*(n+1)/2 ~= 3.5*n**2 + 5*n``,

a larger value should lead to better performance.

Default is None.

Returns

-------

nr : int

`nr` is the order of the resulting reduced order model.

`nr` is set as follows:

If on input ``nr is not None``, `nr` is equal to ``MIN(nr,NMIN)``,

where `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 on input ``nr is None``, `nr` is equal to the number of Hankel

singular values greater than ``MAX(tol,n*eps*HNORM(A,B,C))``.

Ar : (nr, nr) ndarray

This array contains the state dynamics matrix `Ar` of the reduced

order system.

Br : (nr, m) ndarray

Tthis array contains the input/state matrix `Br` of the reduced

order system.

Cr : (p, nr) ndarray

This array contains the state/output matrix `Cr` of the reduced

order system.

hsv : (n, ) ndarray

If ``INFO = 0``, it contains the Hankel singular values of

the original system ordered decreasingly. ``HSV(1)`` is the

Hankel norm of the system.

T : (n, nr) ndarray

This array contains the right truncation matrix `T` of the reduced

order system.

Ti : (nr, n) ndarray

This array contains the left truncation matrix `Ti` of the reduced

order system.

Raises

------

SlycotArithmeticError

:info == 1 and dico == 'C':

The state matrix A is not stable

:info == 1 and dico == 'D':

The state matrix A is not convergent

:info == 2:

The computation of Hankel singular values failed

Warns

-----

SlycotResultWarning

:iwarn == 1:

The selected order {nr} is greater

than the order of a minimal realization of the

given system. `nr` was set automatically to {Nr}

corresponding to the order of a minimal realization

of the system

"""

hidden = ' (hidden by the wrapper)'

arg_list = ['dico', 'job', 'ordsel', 'n', 'm', 'p', 'nr',

'A', 'lda' + hidden, 'B', 'ldb' + hidden, 'C', 'ldc' + hidden,

'hsv', 'T', 'ldt' + hidden, 'Ti', 'ldti' + hidden, 'tol',

'iwork' + hidden, 'dwork' + hidden, 'ldwork', 'iwarn', 'info']

if ldwork is None:

ldwork = max(1, n*(2*n + max(n, max(m, p))+5)+n*(n+1)/2)

if nr is None:

ordsel = 'A'

nr = 0 # order will be computed by the routine

else:

ordsel = 'F'

out = _wrapper.ab09ax(dico, job, ordsel, n, m, p, nr, A, B, C, tol, ldwork)

Nr, A, B, C, hsv, T, Ti = out[:-2]

raise_if_slycot_error(out[-2:], arg_list, ab09ax.__doc__, locals())

return Nr, A[:Nr, :Nr], B[:Nr, :], C[:, :Nr], hsv, T[:, :Nr], Ti[:Nr, :]

def ab09bd(dico,job,equil,n,m,p,A,B,C,D,nr=None,tol1=0,tol2=0,ldwork=None):

""" nr,Ar,Br,Cr,Dr,hsv = ab09bd(dico,job,equil,n,m,p,A,B,C,D,[nr,tol1,tol2,ldwork])

To compute a reduced order model (Ar,Br,Cr,Dr) for a stable

original state-space representation (A,B,C,D) by using either the

square-root or the balancing-free square-root Singular

Perturbation Approximation (SPA) model reduction method.

Must supply either nr or tolerance values.

Parameters

----------

dico : {'C', 'D'}

Specifies the type of the original system as follows:

= 'C': continuous-time system;

= 'D': discrete-time system.

job : {'B', 'N'}

Specifies the model reduction approach to be used

as follows:

= 'B': use the square-root SPA method;

= 'N': use the balancing-free square-root SPA method.

equil : {'S', 'N'}

Specifies whether the user wishes to preliminarily

equilibrate the triplet (A,B,C) as follows:

= 'S': perform equilibration (scaling);

= 'N': do not perform equilibration.

n : int

The order of the original state-space representation, i.e.

the order of the matrix A. n >= 0.

m : int

The number of system inputs. m >= 0.

p : int

The number of system outputs. p >= 0.

A : (n, n) array_like

On entry, the leading n-by-n part of this array must

contain the state dynamics matrix A.

B : (n, m) array_like

On entry, the leading n-by-m part of this array must

contain the original input/state matrix B.

C : (p, n) array_like

On entry, the leading p-by-n part of this array must

contain the original state/output matrix C.

D : (p, m) array_like

On entry, the leading p-by-m part of this array must

contain the original input/output matrix D.

nr : int, optional

nr is the desired order of

the resulting reduced order system. 0 <= nr <= n.

Default is None.

tol1 : float, optional

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.

Default is `0.0`.

tol2 : float, optional

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.

Default is `0.0`.

ldwork : int, optional

The length of the cache array. The default value is n + 3*max(m,p),

for better performance should be larger.

Default is None.

Returns

-------

nr : int

nr is the order of the resulting reduced order model.

nr is set as follows:

if ordsel = 'F', nr is equal to min(nr,nmin), where nr