# Advanced operations on transfer function matrices

Additive spectral decompositions.`gsdec`

Right nullspace basis of a transfer function matrix.`grnull`

Left nullspace basis of a transfer function matrix.`glnull`

Range space basis of a transfer function matrix.`grange`

Coimage space basis of a transfer function matrix.`gcrange`

Solution of the linear rational matrix equation`grsol`

`G(λ)*X(λ) = F(λ)`

.Solution of the linear rational matrix equation`glsol`

`X(λ)*G(λ) = F(λ)`

.Right minimum dynamic cover of Type 1 based order reduction.`grmcover1`

Left minimum dynamic cover of Type 1 based order reduction.`glmcover1`

Right minimum dynamic cover of Type 2 based order reduction.`grmcover2`

Left minimum dynamic cover of Type 2 based order reduction.`glmcover2`

Generalized inverses.`ginv`

`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`

.

`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.

`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.

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 `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.

`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.

`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 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.

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 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`

and `atol2 = atol`

.

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.

`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.

`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)
```

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

`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 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 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.

`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.

`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 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`

.

`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.

`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 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`

.

`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.

`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.

`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`

.

`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:

`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.