Advanced operations on transfer function matrices

gsdec Additive spectral decompositions.
grnull Right nullspace basis of a transfer function matrix.
glnull Left nullspace basis of a transfer function matrix.
grange Range space basis of a transfer function matrix.
gcrange Coimage space basis of a transfer function matrix.
grsol Solution of the linear rational matrix equation G(λ)*X(λ) = F(λ).
glsol Solution of the linear rational matrix equation X(λ)*G(λ) = F(λ).
grmcover1 Right minimum dynamic cover of Type 1 based order reduction.
glmcover1 Left minimum dynamic cover of Type 1 based order reduction.
grmcover2 Right minimum dynamic cover of Type 2 based order reduction.
glmcover2 Left minimum dynamic cover of Type 2 based order reduction.
ginv Generalized inverses.

DescriptorSystems.gsdec — Function

gsdec(sys; job = "finite", prescale, smarg, fast = true,  
      atol = 0,  atol1 = atol, atol2 = atol, rtol = nϵ) -> (sys1, sys2)

Compute for the descriptor system sys = (A-λE,B,C,D) with the transfer function matrix G(λ), the additive spectral decomposition G(λ) = G1(λ) + G2(λ) such that G1(λ), the transfer function matrix of the descriptor system sys1 = (A1-λE1,B1,C1,D1), has only poles in a certain domain of interest Cg of the complex plane and G2(λ), the transfer function matrix of the descriptor system sys2 = (A2-λE2,B2,C2,0), has only poles outside of Cg.

If prescale = true, a preliminary balancing of the descriptor system pair (A,E) is performed. The default setting is prescale = MatrixPencils.balqual(sys.A,sys.E) > 10000, where the function pbalqual from the MatrixPencils package evaluates the scaling quality of the linear pencil A-λE.

The keyword argument smarg, if provided, specifies the stability margin for the stable eigenvalues of A-λE, such that, in the continuous-time case, the stable eigenvalues have real parts less than or equal to smarg, and in the discrete-time case, the stable eigenvalues have moduli less than or equal to smarg. If smarg = missing, the used default values are: smarg = -sqrt(ϵ), for a continuous-time system, and smarg = 1-sqrt(ϵ), for a discrete-time system), where ϵ is the machine precision of the working accuracy.

The keyword argument job, in conjunction with smarg, defines the domain of interest Cg, as follows:

for job = "finite", Cg is the whole complex plane without the point at infinity, and sys1 has only finite poles and sys2 has only infinite poles (default); the resulting A2 is nonsingular and upper triangular, while the resulting E2 is nilpotent and upper triangular;

for job = "infinite", Cg is the point at infinity, and sys1 has only infinite poles and sys2 has only finite poles and is the strictly proper part of sys; the resulting A1 is nonsingular and upper triangular, while the resulting E1 is nilpotent and upper triangular;

for job = "stable", Cg is the stability domain of eigenvalues defined by smarg, and sys1 has only stable poles and sys2 has only unstable and infinite poles; the resulting pairs (A1,E1) and (A2,E2) are in generalized Schur form with E1 upper triangular and nonsingular and E2 upper triangular;

for job = "unstable", Cg is the complement of the stability domain of the eigenvalues defined by smarg, and sys1 has only unstable and infinite poles and sys2 has only stable poles; the resulting pairs (A1,E1) and (A2,E2) are in generalized Schur form with E1 upper triangular and E2 upper triangular and nonsingular.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of A, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A and E. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sys. The keyword argument atol can be used to simultaneously set atol1 = atol, atol2 = atol.

The separation of the finite and infinite eigenvalues is performed using rank decisions based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

source

DescriptorSystems.grnull — Function

grnull(sys; polynomial = false, simple = false, inner = false, fast = true, poles = missing, sdeg = missing,  
       atol = 0, atol1 = atol, atol2 = atol, rtol, offset = sqrt(ϵ) ) -> (sysrnull, info)

Determine for the descriptor systems sys = (A-λE,B,C,D) with the p x m transfer function matrix G(λ), the descriptor system sysrnull = (Ar-λEr,Br,Cr,Dr) with the transfer function matrix Nr(λ) such that Nr(λ) is a minimal rational right nullspace basis of G(λ) and satisfies G(λ)*Nr(λ) = 0.

For the call with

grnull(sys, p2; polynomial = false, simple = false, inner = false, fast = true, poles = missing, sdeg = missing,  
       atol = 0, atol1 = atol, atol2 = atol, rtol, offset = sqrt(ϵ) ) -> (sysrnull, info)

sys contains the compound system sys = [sys1; sys2], with G(λ), the transfer function matrix of sys1, and G2(λ), the transfer function matrix of sys2, and has the descriptor realization sys = (A-λE,B,[C;C2],[D;D2]), where sys2 has p2 outputs. The resulting sysrnull contains the compound system [sysrnull1; sys2*sysrnull1] = (Ar-λEr,Br,[Cr;Cr2],[Dr;Dr2]), where sysrnull1 = (Ar-λEr,Br,Cr,Dr) has the transfer function matrix Nr(λ), which is a rational right nullspace basis of G(λ) satisfying G(λ)*Nr(λ) = 0 and sys2*sysrnull1 = (Ar-λEr,Br,Cr2,Dr2) has the transfer function matrix G2(λ)*Nr(λ).

The returned named tuple info has the components info.nrank, info.stdim, info.degs, info.fnorm and info.tcond.

If polynomial = false, the resulting sysrnull has a proper transfer function matrix, while for polynomial = true the resulting sysrnull has a polynomial transfer function matrix. The resulting basis Nr(λ) contains m-r basis vectors, where r = rank G(λ). The rank r is returned in info.nrank. If simple = true, the resulting basis is simple and satisfies the condition that the sum of the number of poles of the m-r basis vectors is equal to the number of poles of Nr(λ) (i.e., its McMillan degree) .

For a non-simple proper basis, the realization (Ar-λEr,Br,Cr,Dr) is controllable and the pencil [Br Ar-λEr] is in a controllable staircase form. The column dimensions of the full row rank diagonal blocks are returned in info.stdim and the corresponding right Kronecker indices are returned in info.degs. For a simple basis, the regular pencil Ar-λEr is block diagonal, with the i-th block of size info.deg[i] (the i-th right Kronecker index) for a proper basis and info.deg[i]+1 for a polynomial basis. The dimensions of the diagonal blocks are returned in this case in info.stdim, while the increasing numbers of poles of the basis vectors are returned in info.degs. For the i-th basis vector vi(λ) (i.e., the i-th column of Nr(λ)) a minimal realization can be explicitly constructed as (Ari-λEri,Bri,Cr,Dr[:,i]), where Ari, Eri and Bri are the i-th diagonal blocks of Ar, Er, and Br, respectively, and Dr[:,i] is the i-th column of Dr. The corresponding realization of G2(λ)*vi(λ) can be constructed as (Ari-λEri,Bri,Cr2,Dr2[:,i]), where Dr2[:,i] is the i-th column of Dr2.

For a proper basis, the poles of Nr(λ) can be freely assigned, by assigning the eigenvalues of the pencil Ar-λEr. The vector poles, specified as a keyword argument, can be used to specify the desired eigenvalues, alternatively to or jointly with enforcing a desired stability degree sdeg of the real parts of the eigenvalues, for a continuous-time system, or the moduli of eigenvalues, for a discrete-time system. If inner = true, the resulting basis Nr(λ) is inner, i.e., Nr(λ)'*Nr(λ) = I, where Nr(s)' = transpose(Nr(-s)) for a continuous-time system with λ = s and Nr(z)' = transpose(Nr(1/z)) for a discrete-time system with λ = z. If the proper basis is simple, each of the resulting individual basis vector is inner. If sys2 has poles on the boundary of the appropriate stability domain Cs, which are not poles of sys1 too, then there exists no inner Nr(λ) such that G2(λ)*Nr(λ) is stable. An offset can be specified via the keyword parameter offset = β to be used to assess the existence of zeros on the stability domain boundary. Accordingly, for a continuous-time system, the boundary of Cs contains the complex numbers with real parts within the interval [-β,β], while for a discrete-time system, the boundary of Cs contains the complex numbers with moduli within the interval [1-β,1+β]. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of 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 default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sys. The keyword argument atol can be used to simultaneously set atol1 = atol and atol2 = atol.

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false. The computation of simple bases involves the solution of several Type 1 minimum dynamic cover problems. This computation involves using non-orthogonal transformations whose worst condition number is returned in info.tcond, in conjunction with using feedback gains, whose norms are returned in info.fnorm. High values of these quantities indicate a potential loss of numerical stability of computations.

Note: The resulting realization of sysrnull is minimal provided the realization of sys is minimal. However, sysrnull1 is a minimal basis only if the realization (A-lambda E,B,C,D) of sys1 is minimal. In this case, info.degs are the degrees of the vectors of a minimal polynomial basis or, if simple = true, of the resulting minimal simple proper basis.

Method: The computation of a minimal proper right nullspace basis is based on [1]; see also [2]. For the computation of a minimal simple proper right nullspace basis the method of [3] is emloyed to compute a simple basis from a minimal proper basis. For the computation of an inner proper right nullspace basis, the inner factor of an inner-outer factorization of Nr(λ) is explicitly constructed using formulas given in [4].

References:

[1] T.G.J. Beelen. New algorithms for computing the Kronecker structure of a pencil with applications to systems and control theory. Ph. D. Thesis, Technical University Eindhoven, 1987.

[2] A. Varga. On computing least order fault detectors using rational nullspace bases. IFAC SAFEPROCESS'03 Symposium, Washington DC, USA, 2003.

[3] A. Varga. On computing nullspace bases – a fault detection perspective. Proc. IFAC 2008 World Congress, Seoul, Korea, pages 6295–6300, 2008.

[4] K. Zhou, J. C. Doyle, and K. Glover. Robust and Optimal Control. Prentice Hall, 1996.

source

DescriptorSystems.glnull — Function

glnull(sys; polynomial = false, simple = false, coinner = false, fast = true, poles = missing, sdeg = missing,  
       atol = 0, atol1 = atol, atol2 = atol, rtol, offset = sqrt(ϵ) ) -> (syslnull, info)

Determine for the descriptor systems sys = (A-λE,B,C,D) with the p x m transfer function matrix G(λ), the descriptor system syslnull = (Al-λEl,Bl,Cl,Dl) with the transfer function matrix Nl(λ) such that Nl(λ) is a minimal rational left nullspace basis of G(λ) and satisfies Nl(λ)*G(λ) = 0.

For the call with

glnull(sys, m2; polynomial = false, simple = false, coinner = false, fast = true, poles = missing, sdeg = missing,  
       atol = 0, atol1 = atol, atol2 = atol, rtol, offset = sqrt(ϵ) ) -> (syslnull, info)

sys contains the compound system sys = [sys1 sys2], with G(λ), the transfer function matrix of sys1, and G2(λ), the transfer function matrix of sys2, and has the descriptor realization sys = (A-λE,[B B2],C,[D D2]), where sys2 has m2 inputs. The resulting syslnull contains the compound system [syslnull1 syslnull1*sys2] = (Al-λEl,[Bl Bl2],Cr,[Dl Dl2]), where syslnull1 = (Al-λEl,Bl,Cl,Dl) has the transfer function matrix Nl(λ), which is a rational left nullspace basis of G(λ) satisfying Nl(λ)*G(λ) = 0 and syslnull1*sys2 = (Al-λEl,Bl2,Cl,Dl2) has the transfer function matrix Nl(λ)*G2(λ).

The returned named tuple info has the components info.nrank, info.stdim, info.degs, info.fnorm and info.tcond.

If polynomial = false, the resulting syslnull has a proper transfer function matrix, while for polynomial = true the resulting syslnull has a polynomial transfer function matrix. The resulting basis Nl(λ) contains p-r basis vectors, where r = rank G(λ). The rank r is returned in info.nrank. If simple = true, the resulting basis is simple and satisfies the condition that the sum of the number of poles of the p-r basis vectors is equal to the number of poles of Nl(λ) (i.e., its McMillan degree) .

For a non-simple proper basis, the realization (Al-λEl,Bl,Cl,Dl) is observable and the pencil [Al-λEl; Cl] is in an observable staircase form. The row dimensions of the full column rank diagonal blocks are returned in info.stdim and the corresponding left Kronecker indices are returned in info.degs. For a simple basis, the regular pencil Al-λEl is block diagonal, with the i-th block of size info.stdim[i]. The increasing numbers of poles of the basis vectors are returned in info.degs. For the i-th basis vector vi(λ) (i.e., the i-th row of Nl(λ)) a minimal realization can be explicitly constructed as (Ali-λEli,Bl,Cli,Dl[i,:]), where Ali, Eli and Cli are the i-th diagonal blocks of Al, El, and Cl, respectively, and Dl[i,:] is the i-th row of Dl. The corresponding realization of vi(λ)*G2(λ) can be constructed as (Ali-λEli,Bl2,Cl2,Dl2[i,:]), where Dl2[i,:] is the i-th row of Dl2.

For a proper basis, the poles of Nl(λ) can be freely assigned, by assigning the eigenvalues of the pencil Al-λEl. The vector poles, specified as a keyword argument, can be used to specify the desired eigenvalues, alternatively to or jointly with enforcing a desired stability degree sdeg of the real parts of the eigenvalues, for a continuous-time system, or the moduli of eigenvalues, for a discrete-time system. If coinner = true, the resulting basis Nl(λ) is coinner, i.e., Nl(λ)*Nl(λ)' = I, where Nl(s)' = transpose(Nl(-s)) for a continuous-time system with λ = s and Nl(z)' = transpose(Nl(1/z)) for a discrete-time system with λ = z. If the proper basis is simple, each of the resulting individual basis vector is inner. If sys2 has poles on the boundary of the appropriate stability domain Cs, which are not poles of sys1 too, then there exists no inner Nl(λ) such that Nl(λ)*G2(λ) is stable. An offset can be specified via the keyword parameter offset = β to be used to assess the existence of zeros on the stability domain boundary. Accordingly, for a continuous-time system, the boundary of Cs contains the complex numbers with real parts within the interval [-β,β], while for a discrete-time system, the boundary of Cs contains the complex numbers with moduli within the interval [1-β,1+β]. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

Note: The resulting realization of syslnull is minimal provided the realization of sys is minimal. However, syslnull1 is a minimal basis only if the realization (A-lambda E,B,C,D) of sys1 is minimal. In this case, info.degs are the degrees of the vectors of a minimal polynomial basis or, if simple = true, of the resulting minimal simple proper basis.

Method: The computation method for the computation of a right nullspace basis is applied to the dual of descriptor system sys. The computation of a minimal proper right nullspace basis is based on [1]; see also [2]. For the computation of a minimal simple proper right nullspace basis the method of [3] is emloyed to compute a simple basis from a minimal proper basis. For the computation of an inner proper right nullspace basis, the inner factor of an inner-outer factorization of Nl(λ) is explicitly constructed using formulas given in [4].

References:

[1] T.G.J. Beelen. New algorithms for computing the Kronecker structure of a pencil with applications to systems and control theory. Ph. D. Thesis, Technical University Eindhoven, 1987.

[2] A. Varga. On computing least order fault detectors using rational nullspace bases. IFAC SAFEPROCESS'03 Symposium, Washington DC, USA, 2003.

[3] A. Varga. On computing nullspace bases – a fault detection perspective. Proc. IFAC 2008 World Congress, Seoul, Korea, pages 6295–6300, 2008.

[4] K. Zhou, J. C. Doyle, and K. Glover. Robust and Optimal Control. Prentice Hall, 1996.

source

DescriptorSystems.grange — Function

grange(sys; zeros = "none", atol = 0, atol1 = atol, atol2 = atol, rtol, 
       fast = true, offset = sqrt(ϵ)) -> (sysr, sysx, info)

Compute for the descriptor system sys = (A-λE,B,C,D) with the transfer function matrix G(λ), the proper descriptor system sysr = (Ar-λEr,Br,Cr,Dr) with a full column rank transfer function matrix R(λ) such that Range(G(λ)) = Range(R(λ)) and the descriptor system sysx = (A-λE,B,Cx,Dx) with the full row rank transfer function matrix X(λ), which satisfies

 G(λ) = R(λ)*X(λ) ,

representing a full rank factorization of G(λ). The number of columns of R(λ) is the normal rank r of G(λ). The columns of R(λ) form a rational basis of the range (or image) space of the rational matrix G(λ). A selected set of zeros of G(λ) are included as zeros of R(λ).

The resulting named triple info contains (nrank, nfuz, niuz), where info.nrank = r, the normal rank of G(λ), info.nfuz is the number of finite zeros of sys on the boundary of the stability domain Cs, and info.niuz is the number of infinite zeros of sys in the continuous-time case and is set to 0 in the discrete-time case.

Depending on the value of the keyword parameter zeros, the following options can be selected for the zeros of G(λ) to be included in R(λ):

 "none"       - include no zeros (default) 
 "all"        - include all zeros of `sys`
 "unstable"   - include all unstable zeros of `sys`
 "s-unstable" - include all strictly unstable zeros of `sys`, both finite and infinite
 "stable"     - include all stable zeros of `sys`
 "finite"     - include all finite zeros of `sys`
 "infinite"   - include all infinite zeros of `sys`

If inner = true, the resulting basis R(λ) is inner, i.e., R(λ)'*R(λ) = I, where R(s)' = transpose(R(-s)) for a continuous-time system with λ = s and R(z)' = transpose(R(1/z)) for a discrete-time system with λ = z. This option can be used only in conjunction with zeros = "none" or zeros = "unstable".

For a continuous-time system sys, the stability domain Cs is defined as the set of complex numbers with real parts at most -β, while for a discrete-time system sys, Cs is the set of complex numbers with moduli at most 1-β (i.e., the interior of a disc of radius 1-β centered in the origin). The boundary offset β to be used to assess the stability of zeros and their number on the boundary of Cs can be specified via the keyword parameter offset = β. Accordingly, for a continuous-time system, the boundary of Cs contains the complex numbers with real parts within the interval [-β,β], while for a discrete-time system, the boundary of Cs contains the complex numbers with moduli within the interval [1-β,1+β]. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

The keyword arguments atol1, atol2 and rtol, specify, respectively, the absolute tolerance for the nonzero elements of A, B, C and D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, E, B, C and D. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sys. The keyword argument atol can be used to simultaneously set atol1 = atol, atol2 = atol.

For the assessment of zeros, the system pencil [A-λE B; C D] is reduced to a special Kronecker-like form (see [2]). In this reduction, the performed rank decisions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

Method: The range computation method is described in [1] and is based on the reduction algorithm of [2], which has been adapted to deal with several zero selection options. The computation of the involved Kronecker-like form is based on the algorithm of [3].

References:

[1] Varga, A. A note on computing the range of rational matrices. arXiv:1707.0048, https://arxiv.org/abs/1707.0048, 2017.

[2] C. Oara. Constructive solutions to spectral and inner–outer factorizations respect to the disk. Automatica, 41, pp. 1855–1866, 2005.

[3] C. Oara and P. Van Dooren. An improved algorithm for the computation of structural invariants of a system pencil and related geometric aspects. Syst. Control Lett., 30:39–48, 1997.

source

DescriptorSystems.gcrange — Function

gcrange(sys; zeros = "none", coinner = false, atol = 0, atol1 = atol, atol2 = atol, rtol, 
        fast = true, offset = sqrt(ϵ)) -> (sysr, sysx, info)

Compute for the descriptor system sys = (A-λE,B,C,D) with the transfer function matrix G(λ), the proper descriptor system sysr = (Ar-λEr,Br,Cr,Dr) with a full row rank transfer function matrix R(λ) such that Coimage(G(λ)) = Coimage(R(λ)) and the descriptor system sysx = (A-λE,B,Cx,Dx) with the full column rank transfer function matrix X(λ), which satisfies

 G(λ) = X(λ)*R(λ) ,

representing a full rank factorization of G(λ). The number of rows of R(λ) is the normal rank r of G(λ). The rows of R(λ) form a rational basis of the coimage space of the rational matrix G(λ). A selected set of zeros of G(λ) are included as zeros of R(λ).

The following options can be selected via the keyword parameter zeros for which zeros of G(λ) to be included in R(λ):

 "none"       - include no zeros (default) 
 "all"        - include all zeros of `sys`
 "unstable"   - include all unstable zeros of `sys`
 "s-unstable" - include all strictly unstable zeros of `sys`, both finite and infinite
 "stable"     - include all stable zeros of `sys`
 "finite"     - include all finite zeros of `sys`
 "infinite"   - include all infinite zeros of `sys`

If coinner = true, the resulting basis R(λ) is coinner, i.e., R(λ)*R(λ)' = I, where R(s)' = transpose(R(-s)) for a continuous-time system with λ = s and R(z)' = transpose(R(1/z)) for a discrete-time system with λ = z. This option can be used only in conjunction with zeros = "none" or zeros = "unstable".

For the assessment of zeros, the dual system pencil transpose([A-λE B; C D]) is reduced to a special Kronecker-like form (see [2]). In this reduction, the performed rank decisions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

Method: The range computation method described in [1], is applied to the dual descriptor system realization corresponding to the transpose of the rational matrix G(λ). The underlying pencil reduction algorithm of [2], has been adapted to deal with several zero selection options. The computation of the involved Kronecker-like form is based on the algorithm of [3].

References:

[1] Varga, A. A note on computing the range of rational matrices. arXiv:1707.0048, https://arxiv.org/abs/1707.0048, 2017.

[2] C. Oara. Constructive solutions to spectral and inner–outer factorizations respect to the disk. Automatica, 41, pp. 1855–1866, 2005.

[3] C. Oara and P. Van Dooren. An improved algorithm for the computation of structural invariants of a system pencil and related geometric aspects. Syst. Control Lett., 30:39–48, 1997.

source

DescriptorSystems.grsol — Function

grsol(sysg, sysf; poles = missing, sdeg = missing, mindeg = false, solgen = false, minreal = true, fast = true, 
      atol = 0, atol1 = atol, atol2 = atol, rtol, ) -> (sysx, info, sysgen)

Determine for the descriptor systems sysg = (Ag-λEg,Bg,Cg,Dg) and sysf = (Af-λEf,Bf,Cf,Df) with the transfer function matrices G(λ) and F(λ), respectively, the descriptor system sysx with the transfer function matrix X(λ) such that X(λ) is the solution of the linear rational equation

G(λ)X(λ) = F(λ) .      (1)

If solgen = true, the descriptor system sysgen is determined representing a generator of all solutions of (1). Its transfer function matrix has the form GEN(λ) = [ X0(λ) XN(λ) ], such that any X(λ) can be generated as

X(λ) = X0(λ) + XN(λ)*Z(λ) ,

where X0(λ) is a particular solution satisfying G(λ)X0(λ) = F(λ), XN(λ) is a proper rational right nullspace basis of G(λ) satisfying G(λ)XN(λ) = 0, and Z(λ) is an arbitrary rational matrix with suitable dimensions. If solgen = false, sysgen is set to nothing.

The call with

grsol(sysgf, mf; poles = missing, sdeg = missing, mindeg = false, solgen = false, minreal = true, fast = true, 
      atol = 0, atol1 = atol, atol2 = atol, rtol, ) -> (sysx, info, sysgen)

uses the compound descriptor system sysgf = (A-λE,[Bg Bf],C,[Dg Df]), where Bf has mf columns, to define the descriptor systems sysg = (A-λE,Bg,C,Dg) and sysf = (A-λE,Bf,C,Df) (i.e., Ag-λEg = Af-λEf = A-λE and Cg = Cf = C).

The generator sysgen has a descriptor system realization sysgen = (A0-λE0,[B0 BN],C0,[D0 DN]), which is usually not minimal (uncontrollable and/or non-dynamic modes present), with

               ( Ar-λEr    *       *    )  
   A0-λE0    = (   0     Af-λEf    *    ) , 
               (   0       0     Ai-λEi ) 

               ( B1 | Br )
   [B0 | BN] = ( B2 | 0  ),  Cg  =   ( Cr   *    *  ) ,
               ( B3 | 0  )

with Er, Ef and Ai invertible and upper triangular, Ei nillpotent and upper triangular, and DN full row rank. The dimensions of the diagonal blocks of A0-λE0 are returned in the named tuple info as the components info.nr, info.nf, and info.ninf, respectively.

A minimal order descriptor system realization of the proper basis XN(λ) is (Ar-λEr,Br,Cr,DN), where Br and DN have mr columns (returned in info.mr), representing the dimension of the right nullspace basis. The normal rank nrank of G(λ) is returned in info.nrank.

If mindeg = false, the solution sysx is determined in the form sysx = (A0+BN*F-λE0,B0,C0+DN*F,D0), where the matrix F = 0, unless a nonzero stabilizing gain is used such that Ar+Br*F-λEr has stable eigenvalues. The vector poles specified as a keyword argument, can be used to specify the desired eigenvalues alternatively to or jointly with enforcing a desired stability degree sdeg of eigenvalues. The dimension nr of Ar is the number of freely assignable poles of the solution X(λ) and is returned in info.nr. The eigenvalues of Af-λEf contain the finite zeros of G(λ), while the zeros of Ai-λEi contain the infinite zeros of G(λ). The norm of the employed gain F is returned in info.fnorm. If G(λ) has infinite zeros, then the solution X(λ) may have infinite poles. The integer vector info.rdeg contains the relative column degrees of X(λ) (i.e., the numbers of integrators/delays needed to make each column of X(λ) proper).

If mindeg = true, a minimum degree solution is determined as X(λ) = X0(λ) + XN(λ)*Z(λ), where Z(λ) is determined using order reduction based on a Type 2 minimum dynamic cover. This computation involves using non-orthogonal transformations whose worst condition number is returned in info.tcond, in conjunction with using feedback and feedforward gains, whose norms are returned in info.fnorm. High values of these quantities indicate a potential loss of numerical stability of computations.

If minreal = true, the computed realization sysx is minimal.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of Ag, Bg, Cg, Dg, Af, Bf, Cf, Df, the absolute tolerance for the nonzero elements of Eg and Ef, and the relative tolerance for the nonzero elements of Ag, Bg, Cg, Dg, Af, Bf, Cf, Df, Eg and Ef. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the maximal order of the systems sysg and sysf. The keyword argument atol can be used to simultaneously set atol1 = atol, atol2 = atol.

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

Method: The method of [1] to solve rational systems is used.

References:

[1] A. Varga, "Computation of least order solutions of linear rational equations", Proc. MTNS'04, Leuven, Belgium, 2004.

source

DescriptorSystems.glsol — Function

glsol(sysg, sysf; poles = missing, sdeg = missing, mindeg = false, solgen = false, minreal = true, fast = true, 
      atol = 0, atol1 = atol, atol2 = atol, rtol, ) -> (sysx, info, sysgen)

X(λ)G(λ) = F(λ) .      (1)

If solgen = true, the descriptor system sysgen is determined representing a generator of all solutions of (1). Its transfer function matrix has the form GEN(λ) = [ X0(λ); XN(λ) ], such that any X(λ) can be generated as

X(λ) = X0(λ) + Z(λ)*XN(λ) ,

where X0(λ) is a particular solution satisfying X0(λ)G(λ) = F(λ), XN(λ) is a proper rational left nullspace basis of G(λ) satisfying XN(λ)G(λ) = 0, and Z(λ) is an arbitrary rational matrix with suitable dimensions. If solgen = false, sysgen is set to nothing.

The call with

glsol(sysgf, pf; poles = missing, sdeg = missing, mindeg = false, solgen = false, minreal = true, fast = true, 
      atol = 0, atol1 = atol, atol2 = atol, rtol, ) -> (sysx, info, sysgen)

uses the compound descriptor system sysgf = (A-λE,B,[Cg; Cf],[Dg; Df]), where Cf has pf rows, to define the descriptor systems sysg = (A-λE,B,Cg,Dg) and sysf = (A-λE,B,Cf,Df) (i.e., Ag-λEg = Af-λEf = A-λE and Bg = Bf = B).

The generator sysgen has a descriptor system realization sysgen = (A0-λE0,B0, [C0; CN],[D0; DN]), which is usually not minimal (unobservable and/or non-dynamic modes present), with

             ( Ai-λEi    *       *    )  
   A0-λE0  = (   0     Af-λEf    *    ) , 
             (   0       0     Al-λEl ) 

             ( *  )
       B0  = ( *  ),   ( C0 )  = ( C1 C2 C3 ) 
             ( Bl )    ( CN )    ( 0  0  Cl )

with El, Ef and Ai invertible and upper triangular, Ei nillpotent and upper triangular, and DN full column rank. The dimensions of the diagonal blocks of A0-λE0 are returned in the named tuple info as the components info.nf, info.ninf and info.nl, respectively.

A minimal order descriptor system realization of the proper basis XN(λ) is (Al-λEl,Bl,Cl,DN), where Cl and DN have pr columns (returned in info.pr), representing the dimension of the left nullspace basis. The normal rank nrank of G(λ) is returned in info.nrank.

If mindeg = false, the solution sysx is determined in the form sysx = (A0+F*CN-λE0,B0+F*DN,C0,D0), where the matrix F = 0, unless a nonzero stabilizing gain is used such that Al+F*Bl-λEl has stable eigenvalues. The vector poles specified as a keyword argument, can be used to specify the desired eigenvalues alternatively to or jointly with enforcing a desired stability degree sdeg of eigenvalues. The dimension nl of Al is the number of freely assignable poles of the solution X(λ) and is returned in info.nl. The eigenvalues of Af-λEf contain the finite zeros of G(λ), while the zeros of Ai-λEi contain the infinite zeros of G(λ). The norm of the employed gain F is returned in info.fnorm. If G(λ) has infinite zeros, then the solution X(λ) may have infinite poles. The integer vector info.rdeg contains the relative row degrees of X(λ) (i.e., the numbers of integrators/delays needed to make each row of X(λ) proper).

If mindeg = true, a minimum degree solution is determined as X(λ) = X0(λ) + Z(λ)XN(λ), where Z(λ) is determined using order reduction based on a Type 2 minimum dynamic cover. This computation involves using non-orthogonal transformations whose worst condition number is returned in info.tcond, in conjunction with using feedback and feedforward gains, whose norms are returned in info.fnorm. High values of these quantities indicate a potential loss of numerical stability of computations.

If minreal = true, the computed realization sysx is minimal.

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

Method: The dual of method of [1] to solve rational systems is used.

References:

[1] A. Varga, "Computation of least order solutions of linear rational equations", Proc. MTNS'04, Leuven, Belgium, 2004.

source

DescriptorSystems.grmcover1 — Function

grmcover1(sys1, sys2; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

Determine for the proper descriptor systems sys1 = (A1-λE1,B1,C1,D1) and sys2 = (A2-λE2,B2,C2,D2) with the transfer function matrices X1(λ) and X2(λ), respectively, using a right minimum dynamic cover of Type 1 based order reduction, the descriptor systems sysx and sysy with the transfer function matrices X(λ) and Y(λ), respectively, such that

X(λ) = X1(λ) + X2(λ)*Y(λ) ,

and sysx has order less than the order of sys1.

The call with

grmcover1(sys, m1; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

uses the compound descriptor system sys = (A-λE,[B1 B2],C,[D1 D2]), where B1 and D1 have m1 columns, to define the proper descriptor systems sys1 = (A-λE,B1,C,D1) and sys2 = (A-λE,B2,C,D2) (i.e., A1-λE1 = A2-λE2 = A-λE and C1 = C2 = C).

The resulting descriptor systems sysx and sysy have controllable realizations of the form sysx = (Ar-λEr,Br,Cr1,D1) and sysy = (Ar-λEr,Br,Cr2,0), where the pencil [Br Ar-λEr] is in a (controllability) staircase form, with νr[i] x νr[i-1] full row rank diagonal blocks, for i = 1, ..., nr, with νr[0] := m1.

The resulting named triple info contains (stdim, tcond, fnorm), where info.stdim = νr is a vector which contains the row dimensions of the blocks of the staircase form [Br Ar-λEr], info.tcond is the maximum of the Frobenius-norm condition numbers of the employed non-orthogonal transformation matrices, and info.fnorm is the Frobenius-norm of the (internally) employed state-feedback to reduce the order. Large values of info.tcond or info.fnorm indicate possible loss of numerical stability of computations.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of A1, B1, C1, D1, A2, B2, C2, D2, the absolute tolerance for the nonzero elements of E1 and E2, and the relative tolerance for the nonzero elements of A1, B1, C1, D1, A2, B2, C2, D2, E1 and E2. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the maximal order of the systems sys1 and sys2. The keyword argument atol can be used to simultaneously set atol1 = atol, atol2 = atol.

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

Note: grmcover1 also works for arbitrary descriptor system sys1, if sys2 is proper. For an improper system sys1, the order reduction is performed only for the proper part of sys1, while the polynomial part of sys1 is included without modification in the resulting realization of sysx. In this case, info.stdim = νr contains the information corresponding to the proper part of sysx.

Method: The method of [1] is used to compute Type 1 minimum dynamic covers for standard systems and the method of [2] for proper descriptor systems. The resulting order (McMillan degree) of sysx is the least achievable one provided the realization of sys2 is maximally observable (i.e., the pair (A2+B2*F-λE2,C2+D2*F) is observable for any F).

References:

[1] A. Varga, Reliable algorithms for computing minimal dynamic covers, Proc. CDC'03, Maui, Hawaii, 2003.

[2] A. Varga. Reliable algorithms for computing minimal dynamic covers for descriptor systems. Proc. MTNS Symposium, Leuven, Belgium, 2004.

source

DescriptorSystems.glmcover1 — Function

glmcover1(sys1, sys2; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

Determine for the proper descriptor systems sys1 = (A1-λE1,B1,C1,D1) and sys2 = (A2-λE2,B2,C2,D2) with the transfer function matrices X1(λ) and X2(λ), respectively, using a left minimum dynamic cover of Type 1 based order reduction, the descriptor systems sysx and sysy with the transfer function matrices X(λ) and Y(λ), respectively, such that

X(λ) = X1(λ) + Y(λ)*X2(λ) ,

and sysx has order less than the order of sys1.

The call with

glmcover1(sys, p1; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

uses the compound descriptor system sys = (A-λE,B,[C1; C2],[D1; D2]), where C1 and D1 have p1 rows, to define the proper descriptor systems sys1 = (A-λE,B,C1,D1) and sys2 = (A-λE,B,C2,D2) (i.e., A1-λE1 = A2-λE2 = A-λE and B1 = B2 = B).

The resulting descriptor systems sysx and sysy have observable realizations of the form sysx = (Ao-λEo,Bo1,Co,D1) and sysy = (Ao-λEo,Bo2,Co,0), where the pencil [Ao-λEo; Co] is in a (observability) staircase form, with νl[i] x νl[i+1] full row rank diagonal blocks, for i = 1, ..., nl, with νl[nl+1] := p1.

The resulting named triple info contains (stdim, tcond, fnorm), where info.stdim = νl is a vector which contains the column dimensions of the blocks of the staircase form [Ao-λEo; Co], info.tcond is the maximum of the Frobenius-norm condition numbers of the employed non-orthogonal transformation matrices, and info.fnorm is the Frobenius-norm of the (internally) employed output-injection gain to reduce the order. Large values of info.tcond or info.fnorm indicate possible loss of numerical stability of computations.

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

Note: glmcover1 also works for arbitrary descriptor system sys1, if sys2 is proper. For an improper system sys1, the order reduction is performed only for the proper part of sys1, while the polynomial part of sys1 is included without modification in the resulting realization of sysx. In this case, info.stdim = νl contains the information corresponding to the proper part of sysx.

Method: The dual of method of [1] is used to compute Type 1 minimum dynamic covers for standard systems and the dual of method of [2] for proper descriptor systems. The resulting McMillan degree of sysx is the least achievable one provided the realization of sys2 is maximally controllable (i.e., the pair (A2+F*C2-λE2,B2+F*D2) is controllable for any F).

References:

[1] A. Varga, Reliable algorithms for computing minimal dynamic covers, Proc. CDC'03, Maui, Hawaii, 2003.

[2] A. Varga. Reliable algorithms for computing minimal dynamic covers for descriptor systems. Proc. MTNS Symposium, Leuven, Belgium, 2004.

source

DescriptorSystems.grmcover2 — Function

grmcover2(sys1, sys2; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

Determine for the proper descriptor systems sys1 = (A1-λE1,B1,C1,D1) and sys2 = (A2-λE2,B2,C2,D2) with the transfer function matrices X1(λ) and X2(λ), respectively, using a right minimum dynamic cover of Type 2 based order reduction, the descriptor systems sysx and sysy with the transfer function matrices X(λ) and Y(λ), respectively, such that

X(λ) = X1(λ) + X2(λ)*Y(λ) ,

and sysx has order less than the order of sys1.

The call with

grmcover2(sys, m1; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

uses the compound descriptor system sys = (A-λE,[B1 B2],C,[D1 D2]), where B1 and D1 haves m1 columns, to define the proper descriptor systems sys1 = (A-λE,B1,C,D1) and sys2 = (A-λE,B2,C,D2) (i.e., A1-λE1 = A2-λE2 =: A-λE and C1 = C2 =: C).

The resulting descriptor systems sysx and sysy have controllable realizations of the form sysx = (Ar-λEr,Br,Cr1,Dr1) and sysy = (Ar-λEr,Br,Cr2,Dr2), where the pencil [Br Ar-λEr] is in a (controllability) staircase form, with νr[i] x νr[i-1] full row rank diagonal blocks, for i = 1, ..., nr, with νr[0] := m1.

The resulting named triple info contains (stdim, tcond, fnorm, gnorm), where info.stdim = νr is a vector which contains the row dimensions of the blocks of the staircase form [Br Ar-λEr], info.tcond is the maximum of the Frobenius-norm condition numbers of the employed non-orthogonal transformation matrices, info.fnorm is the Frobenius-norm of the (internally) employed state-feedback gain to reduce the order, info.gnorm is the Frobenius-norm of the (internally) employed feedforward gain to reduce the order. Large values of info.tcond, info.fnorm or info.gnorm indicate possible loss of numerical stability of computations.

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

Note: grmcover2 also works for arbitrary descriptor system sys1, if sys2 is proper. For an improper system sys1, the order reduction is performed only for the proper part of sys1, while the polynomial part of sys1 is included without modification in the resulting realization of sysx. In this case, info.stdim = νr contains the information corresponding to the proper part of sysx.

Method: The method of [1] is used to compute Type 2 minimum dynamic covers for standard systems and the method of [2] for proper descriptor systems. The resulting McMillan degree of sysx is the least achievable one provided the realization of sys2 is maximally observable (i.e., the pair (A2+B2*F-λE2,C2+D2*F) is observable for any F).

References:

[1] A. Varga, Reliable algorithms for computing minimal dynamic covers, Proc. CDC'03, Maui, Hawaii, 2003.

[2] A. Varga. Reliable algorithms for computing minimal dynamic covers for descriptor systems. Proc. MTNS Symposium, Leuven, Belgium, 2004.

source

DescriptorSystems.glmcover2 — Function

glmcover2(sys1, sys2; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

Determine for the proper descriptor systems sys1 = (A1-λE1,B1,C1,D1) and sys2 = (A2-λE2,B2,C2,D2) with the transfer function matrices X1(λ) and X2(λ), respectively, using a left minimum dynamic cover of Type 2 based order reduction, the descriptor systems sysx and sysy with the transfer function matrices X(λ) and Y(λ), respectively, such that

X(λ) = X1(λ) + Y(λ)*X2(λ) ,

and sysx has order less than the order of sys1.

The call with

glmcover2(sys, p1; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

uses the compound descriptor system sys = (A-λE,B, [C1; C2],[D1; D2]), where C1 and D1 have p1 rows and E is invertible, to define the proper descriptor systems sys1 = (A-λE,B,C1,D1) and sys2 = (A-λE,B,C2,D2) (i.e., A1-λE1 = A2-λE2 = A-λE and B1 = B2 = B).

The resulting descriptor systems sysx and sysy have observable realizations of the form sysx = (Ao-λEo,Bo1,Co,Do1) and sysy = (Ao-λEo,Bo2,Co,Do2), where the pencil [Ao-λEo; Co] is in a (observability) staircase form, with νl[i] x νl[i+1] full row rank diagonal blocks, for i = 1, ..., nl, with νl[nl+1] := p1.

The resulting named triple info contains (stdim, tcond, fnorm, gnorm), where info.stdim = νl is a vector which contains the column dimensions of the blocks of the staircase form [Ao-λEo; Co], info.tcond is the maximum of the Frobenius-norm condition numbers of the employed non-orthogonal transformation matrices, info.fnorm is the Frobenius-norm of the (internally) employed output-injection gain to reduce the order, and info.gnorm is the Frobenius-norm of the (internally) employed output-feedforward gain. Large values of info.tcond,info.fnorm or info.gnorm indicate possible loss of numerical stability of computations.

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

Note: glmcover2 also works for arbitrary descriptor system sys1, if sys2 is proper. For an improper system sys1, the order reduction is performed only for the proper part of sys1, while the polynomial part of sys1 is included without modification in the resulting realization of sysx. In this case, info.stdim = νl contains the information corresponding to the proper part of sysx.

Method: The dual of method of [1] is used to compute Type 2 minimum dynamic covers for standard systems and the dual of method of [2] for proper descriptor systems. The resulting McMillan degree of sysx is the least achievable one provided the realization of sys2 is maximally controllable (i.e., the pair (A2+F*C2-λE2,B2+F*D2) is controllable for any F).

References:

[1] A. Varga, Reliable algorithms for computing minimal dynamic covers, Proc. CDC'03, Maui, Hawaii, 2003.

[2] A. Varga. Reliable algorithms for computing minimal dynamic covers for descriptor systems. Proc. MTNS Symposium, Leuven, Belgium, 2004.

source

DescriptorSystems.ginv — Function

ginv(sys; type = "1-2", mindeg = false, fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol, 
          offset = sqrt(ϵ)) -> (sysinv, info)

Compute for a descriptor system sys = (A-λE,B,C,D) with the transfer function matrix G(λ) a generalized inverse system sysinv = (Ai-λEi,Bi,Ci,Di) with the transfer function matrix Gi(λ) such that two or more of the following Moore-Penrose conditions are satisfied:

     (1) G(λ)*Gi(λ)*G(λ) = G(λ);        
     (2) Gi(λ)*G(λ)*Gi(λ) = Gi(λ);
     (3) G(λ)*Gi(λ) = (G(λ)*Gi(λ))';
     (4) Gi(λ)*G(λ) = (Gi(λ)*G(λ))'.

The desired type of the computed generalized inverse can be specified using the keyword parameter type as follows:

  "1-2"     - for a generalized inverse which satisfies conditions (1) and (2) (default);
  "1-2-3"   - for a generalized inverse which satisfies conditions (1), (2) and (3);
  "1-2-4"   - for a generalized inverse which satisfies conditions (1), (2) and (4);
  "1-2-3-4" - for the Moore-Penrose pseudoinverse, which satisfies all conditions (1)-(4).

The vector poles specified as a keyword argument, can be used to specify the desired eigenvalues alternatively to or jointly with enforcing a desired stability degree sdeg of the poles of the computed generalized inverse.

The keyword argument mindeg can be used to specify the option to determine a minimum order generalized inverse, if mindeg = true, or a particular generalized inverse which has possibly non-minimal order, if mindeg = false (default).

To assess the presence of zeros on the boundary of the stability domain Cs, a boundary offset β can be specified via the keyword parameter offset = β. Accordingly, for a continuous-time setting, the boundary of Cs contains the complex numbers with real parts within the interval [-β,β], while for a discrete-time setting, the boundary of Cs contains the complex numbers with moduli within the interval [1-β,1+β]. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

The returned named tuple info has the components info.nrank, info.nfp, info.fnorm and info.tcond, where: info.nrank is the normal rank of G(λ), info.nfp is the number of freely assignable poles of the inverse Gi(λ), info.fnorm is the maximum of norms of employed feedback gains (also for pole assignment) and info.tcond is the maximum of condition numbers of employed non-orthogonal transformations (see below).

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting, if fast = true, or the more reliable SVD-decompositions, if fast = false. The computation of a minimum order inverse is performed by solving suitable minimum dynamic cover problems. These computations involve using non-orthogonal transformations whose maximal condition number is returned in info.tcond, in conjunction with using feedback gains (also for pole assignment), whose maximal norms are returned in info.fnorm. High values of these quantities indicate a potential loss of numerical stability of computations.

The keyword arguments atol1, atol2 and rtol, specify, respectively, the absolute tolerance for the nonzero elements of A, B, C and D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, E, B, C and D. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the maximum dimension of state, input and output vectors of the system sys. The keyword argument atol can be used to simultaneously set atol1 = atol and atol2 = atol.

Method: The methods proposed in [1] are employed in conjunction with full rank factorizations computed using the approach of [2].

References:

[1] A. Varga. Computing generalized inverse systems using matrix pencil methods. Int. J. of Applied Mathematics and Computer Science, vol. 11, pp. 1055-1068, 2001.

[2] Varga, A. A note on computing range space bases of rational matrices. arXiv:1707.0048, https://arxiv.org/abs/1707.00489, 2017.

source