Edit on GitHub

Manipulation of linearizations of polynomial or rational matrices

FunctionDescription
lsbalance!Scaling of a descriptor system based linearization
lsminrealComputation of irreducible descriptor system based linearizations
lsminreal2Computation of irreducible descriptor system based linearizations (potentially more efficient)
lpsminrealComputation of strongly minimal pencil based linearizations
lsequalChecking the equivalence of two descriptor system based linearizations
lpsequalChecking the equivalence of two pencil based linearizations
lsevalEvaluation of the value of the rational matrix corresponding to a descriptor system based linearization
lpsevalEvaluation of the value of the rational matrix corresponding to a pencil based linearization
lps2lsConversion of a pencil based linearization into a descriptor system based linearization
lsbalqualEvaluation of the scaling quality of descriptor system based linearizations
MatrixPencils.lsbalance! — Function
 lsbalance!(A, B, C; withB = true, withC = true, maxred = 16) -> D

Reduce the 1-norm of a system matrix 

         S =  ( A  B )
              ( C  0 )

corresponding to a standard system triple (A,B,C) by balancing. This involves a diagonal similarity transformation inv(D)*A*D to make the rows and columns of 

              diag(inv(D),I) * S * diag(D,I)

as close in norm as possible.

The balancing can be optionally performed on the following particular system matrices: 

    S = A, if withB = false and withC = false ,
    S = ( A  B ), if withC = false,     or    
    S = ( A ), if withB = false .
        ( C )

The keyword argument maxred = s specifies the maximum allowed reduction in the 1-norm of S (in an iteration) if zero rows or columns are present. s must be a positive power of 2, otherwise s is rounded to the nearest power of 2.

Note: This function is a translation of the SLICOT routine TB01ID.f, which represents an extensiom of the LAPACK family *GEBAL.f and also includes the fix proposed in [1].

[1] R.James, J.Langou and B.Lowery. "On matrix balancing and eigenvector computation." ArXiv, 2014, http://arxiv.org/abs/1401.5766

source
 lsbalance!(A, E, B, C; withB = true, withC = true, pow2, maxiter = 100, tol = 1) -> (D1,D2)

Reduce the 1-norm of the matrix 

         S =  ( abs(A)+abs(E)  abs(B) )
              (    abs(C)        0    )

corresponding to a descriptor system triple (A-λE,B,C) by row and column balancing. This involves diagonal similarity transformations D1*(A-λE)*D2 applied iteratively to abs(A)+abs(E) to make the rows and columns of 

              diag(D1,I)  * S * diag(D2,I)

as close in norm as possible.

The balancing can be performed optionally on the following particular system matrices: 

    S = abs(A)+abs(E), if withB = false and withC = false ,
    S = ( abs(A)+abs(E)  abs(B) ), if withC = false,     or    
    S = ( abs(A)+abs(E) ), if withB = false .
        (   abs(C)     )

The keyword argument maxiter = k specifies the maximum number of iterations k allowed in the balancing algorithm.

The keyword argument tol = τ, with τ ≤ 1, specifies the tolerance used in the stopping criterion.

If the keyword argument pow2 = true is specified, then the components of the resulting optimal D1 and D2 are replaced by their nearest integer powers of 2, in which case the scaling of matrices is done exactly in floating point arithmetic. If pow2 = false, the optimal values D1 and D2 are returned.

Note: This function is an extension of the MATLAB function rowcolsums.m of [1].

[1] F.M.Dopico, M.C.Quintana and P. van Dooren, "Diagonal scalings for the eigenstructure of arbitrary pencils", SIMAX, 43:1213-1237, 2022.

source
MatrixPencils.lsminreal — Function
lsminreal(A, B, C; fast = true, atol = 0, rtol, contr = true, obs = true, noseig = true) 
          -> (Ar, Br, Cr, nuc, nuo)

Reduce the linearization (A-λI,B,C,0) of a strictly proper rational matrix to a reduced form (Ar-λI,Br,Cr,0) such that

         -1                 -1
 C*(λI-A)  *B  = Cr*(λI-Ar)  *Br

with the least possible order nr of Ar if contr = true and obs = true. Such a realization is called minimal and satisfies:

 (1) rank[Br Ar-λI] = nr for all λ (controllability);

 (2) rank[Ar-λI; Cr] = nr for all λ (observability).

The achieved dimensional reductions to fulfill conditions (1) and (2) are returned in nuc and nuo, respectively.

Some reduction steps can be skipped by appropriately selecting the keyword arguments contr and obs.

If contr = false, then the controllability condition (1) is not enforced.

If obs = false, then observability condition (2) is not enforced.

To enforce conditions (1)-(2), orthogonal similarity transformations are performed on the matrices of the original linearization (A-λI,B,C,0) to obtain a minimal linearization using structured pencil reduction algorithms, as the fast versions of the reduction techniques of the full row rank pencil [B A-λI] and full column rank pencil [A-λI;C] proposed in [1].

The underlying pencil manipulation algorithms employ rank determinations based on either the use of rank revealing QR-decomposition with column pivoting, if fast = true, or the SVD-decomposition. The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments atol and rtol, specify, respectively, the absolute and relative tolerances for the nonzero elements of matrices A, B, C.

[1] P. Van Dooreen, The generalized eigenstructure problem in linear system theory, IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

source
lsminreal(A, E, B, C, D; fast = true, atol1 = 0, atol2, rtol, contr = true, obs = true, noseig = true) 
          -> (Ar, Er, Br, Cr, Dr, nuc, nuo, nse)

Reduce the linearization (A-λE,B,C,D) of a rational matrix to a reduced form (Ar-λEr,Br,Cr,Dr) such that

         -1                    -1
 C*(λE-A)  *B + D = Cr*(λEr-Ar)  *Br + Dr

with the least possible order nr of Ar-λEr if contr = true, obs = true and nseig = true. Such a realization is called minimal and satisfies:

 (1) rank[Br Ar-λEr] = nr for all finite λ (finite controllability)

 (2) rank[Br Er] = nr (infinite controllability)

 (3) rank[Ar-λEr; Cr] = nr for all finite λ (finite observability)

 (4) rank[Er; Cr] = nr (infinite observability)

 (5) Ar-λEr has no simple infinite eigenvalues

A realization satisfying only conditions (1)-(4) is called irreducible.

The achieved dimensional reductions to fulfill conditions (1) and (2), conditions (3) and (4), and respectively, condition (5) are returned in nuc, nuo, nse.

Some reduction steps can be skipped by appropriately selecting the keyword arguments contr, obs and nseig.

If contr = false, then the controllability conditions (1) and (2) are not enforced.

If obs = false, then observability condition (3) and (4) are not enforced.

If nseig = false, then condition (5) on the lack of simple infinite eigenvalues is not enforced.

To enforce conditions (1)-(4), orthogonal similarity transformations are performed on the matrices of the original linearization (A-λE,B,C,D) to obtain an irreducible linearization using structured pencil reduction algorithms, as the fast versions of the reduction techniques of the full row rank pencil [B A-λE] and full column rank pencil [A-λE;C] proposed in [1]. To enforce condition (5), residualization formulas (see, e.g., [2, page 329]) are employed which involves matrix inversions.

The underlying pencil manipulation algorithms employ rank determinations based on either the use of rank revealing QR-decomposition with column pivoting, if fast = true, or the SVD-decomposition. The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of matrices A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E.

[1] P. Van Dooreen, The generalized eigenstructure problem in linear system theory, IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[2] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.

source
MatrixPencils.lsminreal2 — Function
lsminreal2(A, E, B, C, D; fast = true, atol1 = 0, atol2 = 0, rtol, finite = true, infinite = true, contr = true, obs = true, noseig = true) 
           -> (Ar, Er, Br, Cr, Dr, nuc, nuo, nse)

Reduce the linearization (A-λE,B,C,D) of a rational matrix to a reduced form (Ar-λEr,Br,Cr,Dr) such that

         -1                    -1
 C*(λE-A)  *B + D = Cr*(λEr-Ar)  *Br + Dr

with the least possible order nr of Ar-λEr if finite = true, infinite = true, contr = true, obs = true and nseig = true. Such a realization is called minimal and satisfies:

 (1) rank[Br Ar-λEr] = nr for all finite λ (finite controllability);

 (2) rank[Br Er] = nr (infinite controllability);

 (3) rank[Ar-λEr; Cr] = nr for all finite λ (finite observability);

 (4) rank[Er; Cr] = nr (infinite observability);

 (5) Ar-λEr has no simple infinite eigenvalues.

A realization satisfying only conditions (1)-(4) is called irreducible.

The achieved dimensional reductions to fulfill conditions (1) and (2), conditions (3) and (4), and respectively, condition (5) are returned in nuc, nuo, nse.

Some reduction steps can be skipped by appropriately selecting the keyword arguments contr, obs, finite, infinite and nseig.

If contr = false, then the controllability conditions (1) and (2) are not enforced. If contr = true and finite = true, then the finite controllability condition (1) is enforced. If contr = true and finite = false, then the finite controllability condition (1) is not enforced. If contr = true and infinite = true, then the infinite controllability condition (2) is enforced. If contr = true and infinite = false, then the infinite controllability condition (2) is not enforced.

If obs = false, then observability condition (3) and (4) are not enforced. If obs = true and finite = true, then the finite observability condition (3) is enforced. If obs = true and finite = false, then the finite observability condition (3) is not enforced. If obs = true and infinite = true, then the infinite observability condition (4) is enforced. If obs = true and infinite = false, then the infinite observability condition (4) is not enforced.

If nseig = false, then condition (5) on the lack of simple infinite eigenvalues is not enforced.

To enforce conditions (1)-(4), the Procedure GIR in [1, page 328] is employed, which performs orthogonal similarity transformations on the matrices of the original linearization (A-λE,B,C,D) to obtain an irreducible linearization using structured pencil reduction algorithms. To enforce condition (5), residualization formulas (see, e.g., [1, page 329]) are employed which involves matrix inversions.

The underlying pencil manipulation algorithms employ rank determinations based on either the use of rank revealing QR-decomposition with column pivoting, if fast = true, or the SVD-decomposition. The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of matrices A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E.

[1] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.

source
MatrixPencils.lpsminreal — Function
lpsminreal(A, E, B, F, C, G, D, H; fast = true, atol1 = 0, atol2, rtol, contr = true, obs = true)  
           -> (Ar, Er, Br, Fr, Cr, Gr, Dr, Hr, V, W, nuc, nuo)

Reduce the linearization (A-λE,B-λF,C-λG,D-λH) of a rational matrix to a reduced form (Ar-λEr,Br-λFr,Cr-λGr,Dr-λHr) such that, for appropriate invertible upper triangular matrices V and W, 

                  -1                                     -1
 V'*((C-λG)*(λE-A)  *(B-λF) + D-λH)*W = (Cr-λGr)*(λEr-Ar)  *(Br-λFr) + Dr-λHr

with the least possible order nr of Ar-λEr if contr = true and obs = true. Such a realization is called strongly minimal and satisfies:

 (1) rank[Br-λFr Ar-λEr] = nr for all finite and infinite λ (strong controllability)

 (2) rank[Ar-λEr; Cr-λGr] = nr for all finite and infinite λ (strong observability)

The achieved dimensional reductions to fulfill conditions (1) and (2) are returned in nuc and nuo, respectively.

If contr = true, then the strong controllability condition (1) is enforced and W is an invertible upper triangular matrix or W = I if nuc = 0. If contr = false, then the strong controllability condition (1) is not enforced and W = I.

If obs = true, then the strong observability condition (2) is enforced and V is an invertible upper triangular matrix or V = I if nuo = 0. If obs = false, then the strong observability condition (2) is not enforced and V = I.

To enforce conditions (1) and (2), orthogonal similarity transformations are performed on the matrices of the original linearization (A-λE,B-λF,C-λG,D-λH) to obtain a strongly minimal linearization using structured pencil reduction algorithms [1]. The resulting realization (Ar-λEr,Br-λFr,Cr-λGr,Dr-λHr) fulfills the strong controllability and strong observability conditions established in [2].

The underlying pencil manipulation algorithms employ rank determinations based on either the use of rank revealing QR-decomposition with column pivoting, if fast = true, or the SVD-decomposition, if fast = false. The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of matrices A, B, C, D, the absolute tolerance for the nonzero elements of E, F, G, H and the relative tolerance for the nonzero elements of A, B, C, D and E, F, G, H.

[1] F.M. Dopico, M.C. Quintana and P. Van Dooren, Linear system matrices of rational transfer functions, to appear in "Realization and Model Reduction of Dynamical Systems", A Festschrift to honor the 70th birthday of Thanos Antoulas", Springer-Verlag. arXiv:1903.05016

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’, Int. J. Control, Vol.31(5) (1980) 1007–1009.

source
MatrixPencils.lsequal — Function
 lsequal(A1, E1, B1, C1, D1, A2, E2, B2, C2, D2; 
         fastrank = true, atol1 = 0, atol2 = 0, rtol = min(atol1,atol2)>0 ? 0 : n*ϵ) -> flag::Bool

Check if two descriptor system linearizations (A1-λE1,B1,C1,D1) and (A2-λE2,B2,C2,D2) satisfy the equivalence condition 

            -1                      -1
 C1*(λE1-A1)  *B1 + D1 = C2*(λE2-A2)  *B2 + D2

The ckeck is performed by computing the normal rank n of the structured matrix pencil M - λN

          | A1-λE1   0    |   B1  | 
          |   0    A2-λE2 |   B2  | 
 M - λN = |---------------|-------|
          |   C1     -C2  | D1-D2 |

and verifying that n = n1+n2, where n1 and n2 are the orders of the square matrices A1 and A2, respectively.

If fastrank = true, the rank is evaluated by counting how many singular values of M - γ N have magnitude greater than max(max(atol1,atol2), rtol*σ₁), where σ₁ is the largest singular value of M - γ N and γ is a randomly generated value [1]. If fastrank = false, the rank is evaluated as nr + ni + nf + nl, where nr and nl are the sums of right and left Kronecker indices, respectively, while ni and nf are the number of infinite and finite eigenvalues, respectively. The sums nr+ni and nf+nl, are determined from an appropriate Kronecker-like form (KLF) exhibiting the spliting of the right and left structures of the pencil M - λN. For efficiency purpose, the reduction to the relevant KLF is only partially performed using rank decisions based on rank revealing SVD-decompositions.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of M, the absolute tolerance for the nonzero elements of N and the relative tolerance for the nonzero elements of M and N. The default relative tolerance is k*ϵ, where k is the size of the smallest dimension of M, and ϵ is the machine epsilon of the element type of M.

[1] A. Varga, On checking null rank conditions of rational matrices, arXiv:1812.11396, 2018.

source
MatrixPencils.lpsequal — Function
 lpsequal(A1, E1, B1, F1, C1, G1, D1, H1, 
          A2, E2, B2, F2, C2, G2, D2, H2; fastrank = true, atol1 = 0, atol2 = 0, rtol = min(atol1,atol2)>0 ? 0 : n*ϵ) 
          -> flag::Bool

Check if two pencil based linearizations (A1-λE1,B1-λF1,C1-λG1,D1-λH1) and (A2-λE2,B2-λF2,C2-λG2,D2-λH2) satisfy the equivalence condition 

                  -1                                      -1
 (C1-λG1)*(λE1-A1)  *(B1-λF1) + D1-λH1 = (C2-λG2)*(λE2-A2)  *(B2-λF2) + D2-λH2

The ckeck is performed by computing the normal rank n of the structured matrix pencil M - λN

          | A1-λE1    0    |    B1-λF1      | 
          |   0     A2-λE2 |    B2-λF2      | 
 M - λN = |----------------|----------------|
          | C1-λG1 -C2+λG2 | D1-D2-λ(H1-H2) |

and verifying that n = n1+n2, where n1 and n2 are the orders of the square matrices A1 and A2, respectively.

If fastrank = true, the rank is evaluated by counting how many singular values of M - γ N have magnitude greater than max(max(atol1,atol2), rtol*σ₁), where σ₁ is the largest singular value of M - γ N and γ is a randomly generated value [1]. If fastrank = false, the rank is evaluated as nr + ni + nf + nl, where nr and nl are the sums of right and left Kronecker indices, respectively, while ni and nf are the number of infinite and finite eigenvalues, respectively. The sums nr+ni and nf+nl, are determined from an appropriate Kronecker-like form (KLF) exhibiting the spliting of the right and left structures of the pencil M - λN. For efficiency purpose, the reduction to the relevant KLF is only partially performed using rank decisions based on rank revealing SVD-decompositions.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of M, the absolute tolerance for the nonzero elements of N, and the relative tolerance for the nonzero elements of M and N. The default relative tolerance is k*ϵ, where k is the size of the smallest dimension of M, and ϵ is the machine epsilon of the element type of M.

[1] A. Varga, On checking null rank conditions of rational matrices, arXiv:1812.11396, 2018.

source
MatrixPencils.lseval — Function
lseval(A, E, B, C, D, val; atol1, atol2, rtol, fast = true) -> Gval

Evaluate Gval, the value of the rational matrix G(λ) = C*inv(λE-A)*B+D for λ = val. The computed Gval has infinite entries if val is a pole (finite or infinite) of G(λ). If val is finite and val*E-A is singular or if val = Inf and E is singular, then the entries of Gval are evaluated separately for minimal realizations of each input-output channel.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of matrices A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E.

The computation of minimal realizations of individual input-output channels relies on pencil manipulation algorithms, which employ rank determinations based on either the use of rank revealing QR-decomposition with column pivoting, if fast = true, or the SVD-decomposition. The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

source
MatrixPencils.lpseval — Function
lpseval(A, E, B, F, C, G, D, H, val; atol1, atol2, rtol, fast = true) -> Gval

Evaluate Gval, the value of the rational matrix G(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH for λ = val. The computed Gval has infinite entries if val is a pole (finite or infinite) of G(λ). If val is finite and val*E-A is singular or val is infinite and E is singular, then the entries of Gval are evaluated separately for each input-output chanel employing descriptor system based minimal realizations.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of matrices A, B, C, D, the absolute tolerance for the nonzero elements of E, F, G, H, and the relative tolerance for the nonzero elements of A, B, C, D, E, F, G, and H.

The computation of minimal realizations of individual input-output chanels relies on pencil manipulation algorithms, which employ rank determinations based on either the use of rank revealing QR-decomposition with column pivoting, if fast = true, or the SVD-decomposition. The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

source
MatrixPencils.lps2ls — Function
 lps2ls(A, E, B, F, C, G, D, H; compacted = false, atol1 = 0, atol2 = 0, atol3 = 0, rtol = min(atol1,atol2,atol3)>0 ? 0 : n*ϵ) -> (Ad,Ed,Bd,Cd,Dd)

Construct an input-output equivalent descriptor system linearizations (Ad-λdE,Bd,Cd,Dd) to a pencil based linearization (A-λE,B-λF,C-λG,D-λH) satisfying 

            -1                        -1
 Cd*(λEd-Ad)  *Bd + Dd = (C-λG)*(λE-A)  *(B-λF) + D-λH .

If compacted = true, a compacted linearization is determined by exploiting possible rank defficiencies of the matrices F, G, and H.

The keyword arguments atol1, atol2, atol3, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of F, the absolute tolerance for the nonzero elements of G, the absolute tolerance for the nonzero elements of H and the relative tolerance for the nonzero elements of F, G and H. The default relative tolerance is n*ϵ, where n is the size of A, and ϵ is the machine epsilon of the element type of A.

source
MatrixPencils.lsbalqual — Function
qs = lsbalqual(A, B, C; SysMat = false)

Compute the 1-norm based scaling quality of a standard system triple (A,B,C).

If SysMat = false, the resulting qs is computed as 

    qs = max(qS(A),qS(B),qS(C)) ,

where qS(⋅) is the scaling quality measure defined in Definition 5.5 of [1] for nonnegative matrices. This definition has been extended to also cover matrices with zero rows or columns.

If SysMat = true, the resulting qs is computed as 

    qs = qS(S) ,

where S is the system matrix defined as 

         S =  ( A  B )
              ( C  0 )

A large value of qs indicates a possible poorly scaled state-space model.

[1] F.M.Dopico, M.C.Quintana and P. van Dooren, "Diagonal scalings for the eigenstructure of arbitrary pencils", SIMAX, 43:1213-1237, 2022.

source
qs = lsbalqual(A, E, B, C; SysMat = false)

Compute the 1-norm based scaling quality of a descriptor system triple (A-λE,B,C).

If SysMat = false, the resulting qs is computed as 

    qs = max(qS(A),qS(E),qS(B),qS(C)) ,

where qS(⋅) is the scaling quality measure defined in Definition 5.5 of [1] for nonnegative matrices. This definition has been extended to also cover matrices with zero rows or columns.

If SysMat = true, the resulting qs is computed as 

    qs = qS(S) ,

where S is the system matrix defined as 

         S =  ( abs(A)+abs(E)  abs(B) )
              (    abs(C)        0    )

A large value of qs indicates a possible poorly scaled state-space model.

[1] F.M.Dopico, M.C.Quintana and P. van Dooren, "Diagonal scalings for the eigenstructure of arbitrary pencils", SIMAX, 43:1213-1237, 2022.

source