# Factorization of descriptor system transfer function matrices

Right coprime factorization with proper and stable factors.`grcf`

Left coprime factorization with proper and stable factors.`glcf`

Right coprime factorization with inner denominator.`grcfid`

Left coprime factorization with inner denominator.`glcfid`

Normalized right coprime factorization.`gnrcf`

Normalized left coprime factorization.`gnlcf`

Inner-outer/QR-like factorization.`giofac`

Co-outer-co-inner/RQ-like factorization.`goifac`

Right spectral factorization of`grsfg`

`γ^2*I-G'*G`

.Left spectral factorization of`glsfg`

`γ^2*I-G*G'`

.

`DescriptorSystems.grcf`

— Function```
grcf(sys; smarg, sdeg, evals, mindeg = false, mininf = false, fast = true,
atol = 0, atol1 = atol, atol2 = atol, atol3 = atol, rtol = n*ϵ) -> (sysn, sysm)
```

Compute for the descriptor system `sys = (A-λE,B,C,D)`

, the factors `sysn = (An-λEn,Bn,Cn,Dn)`

and `sysm = (Am-λEm,Bm,Cm,Dm)`

of its stable and proper right coprime factorization. If `sys`

, `sysn`

and `sysm`

have the transfer function matrices `G(λ)`

, `N(λ)`

and `M(λ)`

, respectively, then `G(λ) = N(λ)*inv(M(λ))`

, with `N(λ)`

and `M(λ)`

proper and stable transfer function matrices. The resulting matrix pairs `(An,En)`

and `(Am,Em)`

are in (generalized) Schur form. The stability domain `Cs`

of poles is defined by the keyword argument `smarg`

for the stability margin, as follows: for a continuous-time system `sys`

, `Cs`

is the set of complex numbers with real parts at most `smarg < 0`

, while for a discrete-time system `sys`

, `Cs`

is the set of complex numbers with moduli at most `smarg < 1`

(i.e., the interior of a disc of radius `smarg`

centered in the origin). If `smarg`

is missing, then the employed default values are `smarg = -sqrt(eps)`

for a continuous-time system and `smarg = 1-sqrt(eps)`

for a discrete-time system.

The keyword argument `sdeg`

specifies the prescribed stability degree for the assigned eigenvalues of the factors. If both `sdeg`

and `smarg`

are missing, then the employed default values are `sdeg = -0.05`

for a continuous-time system and `sdeg = 0.95`

for a discrete-time system, while if `smarg`

is specified, then `sdeg = smarg`

is used.

The keyword argument `evals`

is a real or complex vector, which contains a set of finite desired eigenvalues for the factors. For a system with real data, `evals`

must be a self-conjugated complex set to ensure that the resulting factors are also real.

If `mindeg = false`

, both factors `sysn`

and `sysm`

have descriptor realizations with the same order and with `An = Am`

, `En = Em`

and `Bn = Bm`

. If `mindeg = true`

, the realization of `sysm`

is minimal. The number of (finite) poles of `sysm`

is equal to the number of unstable finite poles of `sys`

.

If `mininf = false`

, then `An-λEn`

and `Am-λEm`

may have simple infinite eigenvalues. If `mininf = true`

, then `An-λEn`

and `Am-λEm`

have no simple infinite eigenvalues. Note that the removing of simple infinite eigenvalues involves matrix inversions.

The keyword arguments `atol1`

, `atol2`

, `atol3`

, and `rtol`

, specify, respectively, the absolute tolerance for the nonzero elements of `A`

, the absolute tolerance for the nonzero elements of `E`

, the absolute tolerance for the nonzero elements of `B`

, and the relative tolerance for the nonzero elements of `A`

, `E`

and `B`

. The default relative tolerance is `n*ϵ`

, where `ϵ`

is the machine epsilon of the element type of `A`

and `n`

is the order of the system `sys`

. The keyword argument `atol`

can be used to simultaneously set `atol1 = atol`

, `atol2 = atol`

, `atol3 = atol`

.

The preliminary separation of finite and infinite eigenvalues of `A-λE`

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`

.

*Method:* The Procedure GRCF from [2] is implemented, which represents an extension of the recursive factorization approach of [1] to cope with infinite eigenvalues. All infinite poles are assigned to finite real values. If `evals`

is missing or does not contain a sufficient number of real values, then a part or all of infinite eigenvalues of `A-λE`

are assigned to the value specified by `sdeg`

. The pairs `(An,En)`

and `(Am,Em)`

result in *generalized Schur form* with both `An`

and `Am`

quasi-upper triangular and `En`

and `Em`

either both upper triangular or both UniformScalings.

*References:*

[1] A. Varga. Computation of coprime factorizations of rational matrices. Linear Algebra and Its Applications, vol. 271, pp.88-115, 1998.

[2] A. Varga. On recursive computation of coprime factorizations of rational matrices. arXiv:1703.07307, https://arxiv.org/abs/1703.07307, 2020. (to appear in Linear Algebra and Its Applications)

`DescriptorSystems.glcf`

— Function```
glcf(sys; smarg, sdeg, evals, mindeg = false, mininf = false, fast = true,
atol = 0, atol1 = atol, atol2 = atol, atol3 = atol, rtol = n*ϵ) -> (sysn, sysm)
```

Compute for the descriptor system `sys = (A-λE,B,C,D)`

, the factors `sysn = (An-λEn,Bn,Cn,Dn)`

and `sysm = (Am-λEm,Bm,Cm,Dm)`

of its stable and proper left coprime factorization. If `sys`

, `sysn`

and `sysm`

have the transfer function matrices `G(λ)`

, `N(λ)`

and `M(λ)`

, respectively, then `G(λ) = inv(M(λ))*N(λ)`

, with `N(λ)`

and `M(λ)`

proper and stable transfer function matrices. The resulting matrix pairs `(An,En)`

and `(Am,Em)`

are in (generalized) Schur form. The stability domain `Cs`

of poles is defined by the keyword argument `smarg`

for the stability margin, as follows: for a continuous-time system `sys`

, `Cs`

is the set of complex numbers with real parts at most `smarg`

, while for a discrete-time system `sys`

, `Cs`

is the set of complex numbers with moduli at most `smarg < 1`

(i.e., the interior of a disc of radius `smarg`

centered in the origin). If `smarg`

is missing, then the employed default values are `smarg = -sqrt(eps)`

for a continuous-time system and `smarg = 1-sqrt(eps)`

for a discrete-time system.

The keyword argument `sdeg`

specifies the prescribed stability degree for the assigned eigenvalues of the factors. If both `sdeg`

and `smarg`

are missing, then the employed default values are `sdeg = -0.05`

for a continuous-time system and `sdeg = 0.95`

for a discrete-time system, while if `smarg`

is specified, then `sdeg = smarg`

is used.

The keyword argument `evals`

is a real or complex vector, which contains a set of finite desired eigenvalues for the factors. For a system with real data, `evals`

must be a self-conjugated complex set to ensure that the resulting factors are also real.

If `mindeg = false`

, both factors `sysn`

and `sysm`

have descriptor realizations with the same order and with `An = Am`

, `En = Em`

and `Cn = Cm`

. If `mindeg = true`

, the realization of `sysm`

is minimal. The number of (finite) poles of `sysm`

is equal to the number of unstable finite poles of `sys`

.

If `mininf = false`

, then `An-λEn`

and `Am-λEm`

may have simple infinite eigenvalues. If `mininf = true`

, then `An-λEn`

and `Am-λEm`

have no simple infinite eigenvalues. Note that the removing of simple infinite eigenvalues involves matrix inversions.

The keyword arguments `atol1`

, `atol2`

, `atol3`

, and `rtol`

, specify, respectively, the absolute tolerance for the nonzero elements of `A`

, the absolute tolerance for the nonzero elements of `E`

, the absolute tolerance for the nonzero elements of `C`

, and the relative tolerance for the nonzero elements of `A`

, `E`

and `C`

. The default relative tolerance is `n*ϵ`

, where `ϵ`

is the machine epsilon of the element type of `A`

and `n`

is the order of the system `sys`

. The keyword argument `atol`

can be used to simultaneously set `atol1 = atol`

, `atol2 = atol`

, `atol3 = atol`

.

The preliminary separation of finite and infinite eigenvalues of `A-λE`

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`

.

*Method:* The dual of Procedure GRCF from [2] is used, which represents an extension of the recursive factorization approach of [1] to cope with infinite poles. All infinite eigenvalues are assigned to finite real values. If `evals`

is missing or does not contain a sufficient number of real values, then a part or all of infinite eigenvalues of `A-λE`

are assigned to the value specified by `sdeg`

. The pairs `(An,En)`

and `(Am,Em)`

result in *generalized Schur form* with both `An`

and `Am`

quasi-upper triangular and `En`

and `Em`

either both upper triangular or both UniformScalings.

*References:*

[1] A. Varga. Computation of coprime factorizations of rational matrices. Linear Algebra and Its Applications, vol. 271, pp.88-115, 1998.

[2] A. Varga. On recursive computation of coprime factorizations of rational matrices. arXiv:1703.07307, https://arxiv.org/abs/1703.07307, 2020. (to appear in Linear Algebra and Its Applications)

`DescriptorSystems.grcfid`

— Function```
grcfid(sys; mindeg = false, mininf = false, fast = true, offset = sqrt(ϵ),
atol = 0, atol1 = atol, atol2 = atol, atol3 = atol, rtol = n*ϵ) -> (sysni, sysmi)
```

Compute for the descriptor system `sys = (A-λE,B,C,D)`

, the factors `sysni = (Ani-λEni,Bni,Cni,Dni)`

and `sysmi = (Ami-λEmi,Bmi,Cmi,Dmi)`

of its right coprime factorization with inner denominator. If `sys`

, `sysni`

and `sysmi`

have the transfer function matrices `G(λ)`

, `N(λ)`

and `M(λ)`

, respectively, then `G(λ) = N(λ)*inv(M(λ))`

, with `N(λ)`

and `M(λ)`

proper and stable transfer function matrices and the denominator factor `M(λ)`

inner. The resulting matrix pairs `(Ani,Eni)`

and `(Ami,Emi)`

are in (generalized) Schur form. The system `sys`

must not have poles on the boundary of the stability domain `Cs`

. In terms of eigenvalues, this requires for a continuous-time system, that `A-λE`

must not have controllable eigenvalues on the imaginary axis (excepting simple infinite eigenvalues), while for a discrete-time system, `A-λE`

must not have controllable eigenvalues on the unit circle centered in the origin.

To assess the presence of poles on the boundary of `Cs`

, a boundary offset `β`

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, then 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.

If `mindeg = false`

, both factors `sysni`

and `sysmi`

have descriptor realizations with the same order and with `Ani = Ami`

, `Eni = Emi`

and `Bni = Bmi`

. If `mindeg = true`

, the realization of `sysmi`

is minimal. The number of (finite) poles of `sysmi`

is equal to the number of unstable finite poles of `sys`

.

If `mininf = false`

, then `Ani-λEni`

and `Ami-λEmi`

may have simple infinite eigenvalues. If `mininf = true`

, then `Ani-λEni`

and `Ami-λEmi`

have no simple infinite eigenvalues. Note that the removing of simple infinite eigenvalues involves matrix inversions.

The keyword arguments `atol1`

, `atol2`

, `atol3`

, and `rtol`

, specify, respectively, the absolute tolerance for the nonzero elements of `A`

, the absolute tolerance for the nonzero elements of `E`

, the absolute tolerance for the nonzero elements of `B`

, and the relative tolerance for the nonzero elements of `A`

, `E`

and `B`

. The default relative tolerance is `n*ϵ`

, where `ϵ`

is the working machine epsilon. The keyword argument `atol`

can be used to simultaneously set `atol1 = atol`

, `atol2 = atol`

, `atol3 = atol`

.

The preliminary separation of finite and infinite eigenvalues of `A-λE`

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`

.

*Method:* An extension of the recursive factorization approach of [1] is used (see [2] for details). The pairs `(Ani,Eni)`

and `(Ami,Emi)`

result in *generalized Schur form* with both `Ani`

and `Ami`

quasi-upper triangular and `Eni`

and `Emi`

either both upper triangular or both UniformScalings.

*References:*

[1] A. Varga. Computation of coprime factorizations of rational matrices. Linear Algebra and Its Applications, vol. 271, pp.88-115, 1998.

[2] A. Varga. On recursive computation of coprime factorizations of rational matrices. arXiv:1703.07307, https://arxiv.org/abs/1703.07307, 2020. (to appear in Linear Algebra and Its Applications)

`DescriptorSystems.glcfid`

— Function```
glcfid(sys; mindeg = false, mininf = false, fast = true, offset = sqrt(ϵ),
atol = 0, atol1 = atol, atol2 = atol, atol3 = atol, rtol = n*ϵ) -> (sysni, sysmi)
```

Compute for the descriptor system `sys = (A-λE,B,C,D)`

, the factors `sysni = (Ani-λEni,Bni,Cni,Dni)`

and `sysmi = (Ami-λEmi,Bmi,Cmi,Dmi)`

of its left coprime factorization with inner denominator. If `sys`

, `sysni`

and `sysmi`

have the transfer function matrices `G(λ)`

, `N(λ)`

and `M(λ)`

, respectively, then `G(λ) = inv(M(λ))*N(λ)`

, with `N(λ)`

and `M(λ)`

proper and stable transfer function matrices and the denominator factor `M(λ)`

inner. The resulting matrix pairs `(Ani,Eni)`

and `(Ami,Emi)`

are in Schur forms. The system `sys`

must not have poles on the boundary of the stability domain `Cs`

. In terms of eigenvalues, this requires for a continuous-time system, that `A-λE`

must not have controllable eigenvalues on the imaginary axis (excepting simple infinite eigenvalues), while for a discrete-time system, `A-λE`

must not have controllable eigenvalues on the unit circle centered in the origin.

To assess the presence of poles on the boundary of `Cs`

, a boundary offset `β`

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, then 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.

If `mindeg = false`

, both factors `sysni`

and `sysmi`

have descriptor realizations with the same order and with `Ani = Ami`

, `Eni = Emi`

and `Cni = Cmi`

. If `mindeg = true`

, the realization of `sysmi`

is minimal. The number of (finite) poles of `sysmi`

is equal to the number of unstable finite poles of `sys`

.

If `mininf = false`

, then `Ani-λEni`

and `Ami-λEmi`

may have simple infinite eigenvalues. If `mininf = true`

, then `Ani-λEni`

and `Ami-λEmi`

have no simple infinite eigenvalues. Note that the removing of simple infinite eigenvalues involves matrix inversions.

The keyword arguments `atol1`

, `atol2`

, `atol3`

, and `rtol`

, specify, respectively, the absolute tolerance for the nonzero elements of `A`

, the absolute tolerance for the nonzero elements of `E`

, the absolute tolerance for the nonzero elements of `C`

, and the relative tolerance for the nonzero elements of `A`

, `E`

and `C`

. The default relative tolerance is `n*ϵ`

, where `ϵ`

is the working machine epsilon. The keyword argument `atol`

can be used to simultaneously set `atol1 = atol`

, `atol2 = atol`

, `atol3 = atol`

.

`A-λE`

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`

.

*Method:* An extension of the recursive factorization approach of [1] is used to the dual system (see [2] for details). The pairs `(Ani,Eni)`

and `(Ami,Emi)`

result in *generalized Schur form* with both `Ani`

and `Ami`

quasi-upper triangular and `Eni`

and `Emi`

either both upper triangular or both UniformScalings.

*References:*

`DescriptorSystems.gnrcf`

— Function```
gnrcf(sys; fast = true, ss = false,
atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (sysn, sysm)
```

Compute for the descriptor system `sys = (A-λE,B,C,D)`

, the factors `sysn = (An-λEn,Bn,Cn,Dn)`

and `sysm = (An-λEn,Bn,Cm,Dm)`

of its normalized right coprime factorization. If `sys`

, `sysn`

and `sysm`

have the transfer function matrices `G(λ)`

, `N(λ)`

and `M(λ)`

, respectively, then `G(λ) = N(λ)*inv(M(λ))`

, with `N(λ)`

and `M(λ)`

proper and stable transfer function matrices and `[N(λ);M(λ)]`

inner. The resulting `En = I`

if `ss = true`

.

Pencil reduction algorithms are employed which perform rank decisions based on rank revealing QR-decompositions with column pivoting if `fast = true`

or the more reliable SVD-decompositions if `fast = false`

.

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 machine epsilon of the element type of `A`

and `n`

is the order of the system `sys`

. The keyword argument `atol`

can be used to simultaneously set `atol1 = atol`

, `atol2 = atol`

.

*Method:* Pencil reduction algorithms are employed to determine the inner range space `R(λ)`

of the transfer function matrix `[G(λ); I]`

using the method described in [1], which is based on the reduction algorithm of [2]. Then the factors `N(λ)`

and `M(λ)`

result from the partitioning of `R(λ)`

as `R(λ) = [N(λ);M(λ)]`

.

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

`DescriptorSystems.gnlcf`

— Function```
gnlcf(sys; fast = true, ss = false,
atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (sysn, sysm)
```

Compute for the descriptor system `sys = (A-λE,B,C,D)`

, the factors `sysn = (An-λEn,Bn,Cn,Dn)`

and `sysm = (An-λEn,Bm,Cn,Dm)`

of its normalized right coprime factorization. If `sys`

, `sysn`

and `sysm`

have the transfer function matrices `G(λ)`

, `N(λ)`

and `M(λ)`

, respectively, then `G(λ) = inv(M(λ))*N(λ)`

, with `N(λ)`

and `M(λ)`

proper and stable transfer function matrices and `[N(λ) M(λ)]`

coinner. The resulting `En = I`

if `ss = true`

.

Pencil reduction algorithms are employed which perform rank decisions based on rank revealing QR-decompositions with column pivoting if `fast = true`

or the more reliable SVD-decompositions if `fast = false`

.

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 machine epsilon of the element type of `A`

and `n`

is the order of the system `sys`

. The keyword argument `atol`

can be used to simultaneously set `atol1 = atol`

, `atol2 = atol`

.

*Method:* Pencil reduction algorithms are employed to determine the coinner coimage space `R(λ)`

of the transfer function matrix `[G(λ) I]`

using the dual of method described in [1], which is based on the reduction algorithm of [2]. Then the factors `N(λ)`

and `M(λ)`

result from the partitioning of `R(λ)`

as `R(λ) = [N(λ) M(λ)]`

.

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

`DescriptorSystems.giofac`

— Function```
giofac(sys; atol = 0, atol1 = atol, atol2 = atol, atol3 = atol, rtol,
fast = true, minphase = true, offset = sqrt(ϵ)) -> (sysi, syso, info)
```

Compute for the descriptor system `sys = (A-λE,B,C,D)`

with the transfer function matrix `G(λ)`

, the square inner factor `sysi = (Ai-λEi,Bi,Ci,Di)`

with the transfer function matrix `Gi(λ)`

and the minimum-phase quasi-outer factor or the full row rank factor `syso = (Ao-λEo,Bo,Co,Do)`

with the transfer function matrix `Go(λ)`

such that

` G(λ) = Gi[:,1:r](λ)*Go(λ) (*),`

where `r`

is the normal rank of `G(λ)`

. The resulting proper and stable inner factor satisfies `Gi'(λ)*Gi(λ) = I`

. If `sys`

is stable (proper), then the resulting `syso`

is stable (proper). The resulting factor `Go(λ)`

has full row rank `r`

. Depending on the selected factorization option, if `minphase = true`

, then `Go(λ)`

is minimum phase, excepting possibly zeros on the boundary of the appropriate stability domain `Cs`

, or if `minphase = false`

, then `Go(λ)`

contains all zeros of `G(λ)`

, in which case (*) is the extended QR-like factorization of `G(λ)`

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

on the boundary of `Cs`

, and `info.niuz`

is the number of infinite zeros of `syso`

. `info.nfuz`

is set to `missing`

if `minphase = false`

.

*Note:* `syso`

may generally contain a *free inner factor*, which can be eliminated by removing the finite unobservable eigenvalues.

The keyword arguments `atol1`

, `atol2`

, `atol3`

, and `rtol`

, specify, respectively, the absolute tolerance for the nonzero elements of `A`

and `B`

, the absolute tolerance for the nonzero elements of `E`

, the absolute tolerance for the nonzero elements of `C`

and `D`

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

, `atol3 = 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:* For a continuous-time system, the factorization algorithm of [1] is used, while for a discrete-time system, the factorization algorithm of [1] is used.

*References:*

[1] C. Oara and A. Varga. Computation of the general inner-outer and spectral factorizations. IEEE Trans. Autom. Control, vol. 45, pp. 2307-2325, 2000.

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

`DescriptorSystems.goifac`

— Function```
goifac(sys; atol = 0, atol1 = atol, atol2 = atol, atol3 = atol, rtol,
fast = true, minphase = true, offset = sqrt(ϵ)) -> (sysi, syso, info)
```

Compute for the descriptor system `sys = (A-λE,B,C,D)`

with the transfer function matrix `G(λ)`

, the square inner factor `sysi = (Ai-λEi,Bi,Ci,Di)`

with the transfer function matrix `Gi(λ)`

and the minimum-phase quasi-outer factor or the full column rank factor `syso = (Ao-λEo,Bo,Co,Do)`

with the transfer function matrix `Go(λ)`

such that

` G(λ) = Go(λ)*Gi[1:r,:](λ) (*),`

where `r`

is the normal rank of `G(λ)`

. The resulting proper and stable inner factor satisfies `Gi'(λ)*Gi(λ) = I`

. If `sys`

is stable (proper), then the resulting `syso`

is stable (proper). The resulting factor `Go(λ)`

has full column rank `r`

. Depending on the selected factorization option, if `minphase = true`

, then `Go(λ)`

is minimum phase, excepting possibly zeros on the boundary of the appropriate stability domain `Cs`

, or if `minphase = false`

, then `Go(λ)`

contains all zeros of `G(λ)`

, in which case (*) is the extended RQ-like factorization of `G(λ)`

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

on the boundary of `Cs`

, and `info.niuz`

is the number of infinite zeros of `syso`

. `info.nfuz`

is set to `missing`

if `minphase = false`

.

*Note:* `syso`

may generally contain a *free inner factor*, which can be eliminated by removing the finite unobservable eigenvalues.

The keyword arguments `atol1`

, `atol2`

, `atol3`

, and `rtol`

, specify, respectively, the absolute tolerance for the nonzero elements of `A`

and `C`

, the absolute tolerance for the nonzero elements of `E`

, the absolute tolerance for the nonzero elements of `B`

and `D`

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

, `atol3 = 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:* For a continuous-time system, the dual system is formed and the factorization algorithm of [1] is used, while for a discrete-time system, the factorization algorithm of [1] is used.

*References:*

[1] C. Oara and A. Varga. Computation of the general inner-outer and spectral factorizations. IEEE Trans. Autom. Control, vol. 45, pp. 2307–2325, 2000.

`DescriptorSystems.grsfg`

— Function```
sysf = grsfg(sys, γ; fast = true, stabilize = true, offset = β,
atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ)
```

Compute for the descriptor system `sys = (A-λE,B,C,D)`

with the transfer function matrix `G(λ)`

and ${\small γ > \|G(λ)\|_∞}$, the minimum-phase right spectral factor `sysf = (Af-λEf,Bf,Cf,Df)`

with the transfer-function matrix `F(λ)`

, such that `F(λ)'*F(λ) = γ^2*I-G(λ)'*G(λ)`

. If `stabilize = true`

(the default), a preliminary stabilization of `sys`

is performed. In this case, `sys`

must not have poles on the imaginary-axis in the continuous-time case or on the unit circle in the discrete-time case. If `stabilize = false`

, then no preliminary stabilization is performed. In this case, `sys`

must be stable.

To assess the presence of poles on the boundary of the stability domain `Cs`

, a boundary offset `β`

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, then 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.

If `stabilize = true`

, a preliminary separation of finite and infinite eigenvalues of `A-λE`

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`

.

The keyword arguments `atol1`

, `atol2`

, `atol3`

, and `rtol`

, specify, respectively, the absolute tolerance for the nonzero elements of `A`

, the absolute tolerance for the nonzero elements of `E`

, the absolute tolerance for the nonzero elements of `C`

, and the relative tolerance for the nonzero elements of `A`

, `E`

and `C`

. The default relative tolerance is `n*ϵ`

, where `ϵ`

is the working machine epsilon. The keyword argument `atol`

can be used to simultaneously set `atol1 = atol`

, `atol2 = atol`

, `atol3 = atol`

.

*Method:* Extensions of the factorization approaches of [1] are used.

*References:*

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

`DescriptorSystems.glsfg`

— Function```
sysf = glsfg(sys, γ; fast = true, stabilize = true, offset = β,
atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ)
```

Compute for the descriptor system `sys = (A-λE,B,C,D)`

with the transfer function matrix `G(λ)`

and ${\small γ > \|G(λ)\|_∞}$, the minimum-phase right spectral factor `sysf = (Af-λEf,Bf,Cf,Df)`

with the transfer-function matrix `F(λ)`

, such that `F(λ)*F(λ)' = γ^2*I-G(λ)*G(λ)'`

. If `stabilize = true`

(the default), a preliminary stabilization of `sys`

is performed. In this case, `sys`

must not have poles on the imaginary-axis in the continuous-time case or on the unit circle in the discrete-time case. If `stabilize = false`

, then no preliminary stabilization is performed. In this case, `sys`

must be stable.

To assess the presence of poles on the boundary of the stability domain `Cs`

, a boundary offset `β`

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, then 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.

If `stabilize = true`

, a preliminary separation of finite and infinite eigenvalues of `A-λE`

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`

.

The keyword arguments `atol1`

, `atol2`

, `atol3`

, and `rtol`

, specify, respectively, the absolute tolerance for the nonzero elements of `A`

, the absolute tolerance for the nonzero elements of `E`

, the absolute tolerance for the nonzero elements of `C`

, and the relative tolerance for the nonzero elements of `A`

, `E`

and `C`

. The default relative tolerance is `n*ϵ`

, where `ϵ`

is the working machine epsilon. The keyword argument `atol`

can be used to simultaneously set `atol1 = atol`

, `atol2 = atol`

, `atol3 = atol`

.

*Method:* Extensions of the factorization approaches of [1] are used.

*References:*

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