On recursive computation of coprime factorizations of rational matrices

TL;DR

This paper introduces recursive pole dislocation methods for computing minimal-degree coprime factorizations of rational matrices, ensuring stability and properness through descriptor state-space realizations.

Contribution

It presents novel recursive pole dislocation techniques that enable minimal-degree coprime factorizations with stable and proper factors, applicable to improper rational matrices.

Findings

Successfully computes coprime factorizations with minimal denominator degree.

Ensures numerical reliability in pole placement and factorization.

Provides implementation details and illustrative examples.

Abstract

General computational methods based on descriptor state-space realizations are proposed to compute coprime factorizations of rational matrices with minimum degree denominators. The new methods rely on recursive pole dislocation techniques, which allow to successively place all poles of the factors into a "good" region of the complex plane. The resulting McMillan degree of the denominator factor is equal to the number of poles lying in the complementary "bad" region and therefore is minimal. The developed pole dislocation techniques are instrumental for devising numerically reliable procedures for the computation of coprime factorizations with proper and stable factors of arbitrary improper rational matrices and coprime factorizations with inner denominators. Implementation aspects of the proposed algorithms are discussed and illustrative examples are given.

Equations

$$\text{mathds}C=\text{mathds}{C}_g\cup \text{mathds}{C}_b,\text{mathds}{C}_g\cap \text{mathds}{C}_b=\varnothing ,$$

$$G(\lambda )=M^{-1}(\lambda )N(\lambda ),$$

$$G(\lambda )=N(\lambda )M^{-1}(\lambda ),$$

$$G(\lambda )=C(\lambda E-A{)}^{-1}B+D.$$

$$G(\lambda)=\left[\begin{array}[]{c|c}A-\lambda E&B\\ \hline\cr C&D\end{array}\right].$$

$$\tilde{A}-\lambda \tilde{E}=Q(A-\lambda E)Z,\tilde{B}=QB,\tilde{C}=CZ,$$

$$\widetilde{A}-\lambda\widetilde{E}={\left[\begin{array}[]{cccc}A_{\infty}&\ast&\ast&\ast\\ 0&A_{g}\!-\!\lambda E_{g}&\ast&\ast\\ 0&0&A_{f,b}\!-\!\lambda E_{f,b}&\ast\\ 0&0&0&A_{\infty,b}\!-\!\lambda E_{\infty,b}\end{array}\right]},$$

$$Q_{1}(A-\lambda E)Z_{1}=\left[\begin{array}[]{cc}A_{\infty}&\ast\\ 0&A_{22}-\lambda E_{22}\end{array}\right],$$

$$Q_{2}(A_{22}-\!\lambda E_{22})Z_{2}={\left[\begin{array}[]{cc}A_{f}\!-\!\lambda E_{f}&\ast\\ 0&A_{\infty,b}\!-\!\lambda E_{\infty,b}\end{array}\right]},$$

$$Q_{3}(A_{f}-\!\lambda E_{f})Z_{3}={\left[\begin{array}[]{cc}A_{g}\!-\!\lambda E_{g}&\ast\\ 0&A_{f,b}\!-\!\lambda E_{f,b}\end{array}\right]},$$

$$Q_{1}A_{2}=\left[\begin{array}[]{c}\ast\\ A_{22}\end{array}\right],\quad Q_{1}E_{2}=\left[\begin{array}[]{c}\ast\\ E_{22}\end{array}\right]\,.$$

$$P_{2}Q_{2}(A_{22}-\lambda E_{22})Z_{2}P_{2}={\left[\begin{array}[]{cc}A_{f}\!-\!\lambda E_{f}&\ast\\ 0&A_{\infty,b}\!-\!\lambda E_{\infty,b}\end{array}\right]},$$

$$Q_{3}(A_{f}-\!\lambda E_{f})Z_{3}={\left[\begin{array}[]{cc}A_{g}\!-\!\lambda E_{g}&\ast\\ 0&A_{f,b}\!-\!\lambda E_{f,b}\end{array}\right]}$$

$$M(\lambda )={\tilde{M}}_1(\lambda ){\tilde{M}}_2(\lambda )\cdots {\tilde{M}}_k(\lambda ),$$

$$\left[\begin{array}[]{c}N_{i}(\lambda)\\ M_{i}(\lambda)\end{array}\right]=\left[\begin{array}[]{c}N_{i-1}(\lambda)\\ M_{i-1}(\lambda)\end{array}\right]\widetilde{M}_{i}(\lambda),\quad i=1,\ldots,k\,,$$

$$\left[\begin{array}[]{c}N_{i-1}(\lambda)\\ M_{i-1}(\lambda)\end{array}\right]=\left[\begin{array}[]{cc|c}A_{11}-\lambda E_{11}&A_{12}-\lambda E_{12}&B_{1}\\ 0&A_{22}-\lambda E_{22}&B_{2}\\ \hline\cr C_{N,1}&C_{N,2}&D_{N}\\ C_{M,1}&C_{M,2}&D_{M}\end{array}\right]=:\left[\begin{array}[]{c|c}\widetilde{A}-\lambda\widetilde{E}&\widetilde{B}\\ \hline\cr\\[-9.95845pt] \widetilde{C}_{N}&\widetilde{D}_{N}\\ \widetilde{C}_{M}&\widetilde{D}_{M}\end{array}\right],$$

$${\left[\begin{array}[]{c}N_{i}(\lambda)\\ M_{i}(\lambda)\end{array}\right]=\left[\begin{array}[]{c|c}A_{11}\!-\!\lambda E_{11}&B_{1}\\ \hline\cr C_{N,1}&D_{N}\\ C_{M,1}&D_{M}\end{array}\right]=:\left[\begin{array}[]{c|c}\widetilde{A}\!-\!\lambda\widetilde{E}&\widetilde{B}\\ \hline\cr\\[-9.95845pt] \widetilde{C}_{N}&\widetilde{D}_{N}\\ \widetilde{C}_{M}&\widetilde{D}_{M}\end{array}\right]}.$$

$${\tilde{M}}_i(\lambda )=(A_{22}+B_2{F}_2-\lambda E_{22},B_{2}W,F_2,W),$$

$${\left[\begin{array}[]{c}N_{i}(\lambda)\\ M_{i}(\lambda)\end{array}\right]=\left[\begin{array}[]{cc|c}A_{11}-\lambda E_{11}&A_{12}+B_{1}F_{2}-\lambda E_{12}&B_{1}W\\ 0&A_{22}+B_{2}F_{2}-\lambda E_{22}&B_{2}W\\ \hline\cr C_{N,1}&C_{N,2}+D_{N}F_{2}&D_{N}W\\ C_{M,1}&C_{M,2}+D_{M}F_{2}&D_{M}W\end{array}\right]}.$$

$$\begin{array}[]{lcl}\widetilde{A}&\leftarrow&\widetilde{A}+\widetilde{B}\widetilde{F},\\ \widetilde{B}&\leftarrow&\widetilde{B}W,\\ \widetilde{C}_{N}&\leftarrow&\widetilde{C}_{N}+\widetilde{D}_{N}\widetilde{F},\\ \widetilde{C}_{M}&\leftarrow&\widetilde{C}_{M}+\widetilde{D}_{M}\widetilde{F},\\ \widetilde{D}_{N}&\leftarrow&\widetilde{D}_{N}W,\\ \widetilde{D}_{M}&\leftarrow&\widetilde{D}_{M}W.\end{array}$$

$${\tilde{M}}_i(\lambda )=(\gamma -\lambda \eta ,B_2,F_2,W),$$

$$\begin{array}[]{l}F_{2}=-B_{2}^{T}(B_{2}B_{2}^{T})^{-1}A_{22},\\ W=I-B_{2}^{T}(B_{2}B_{2}^{T})^{-1}B_{2}.\end{array}$$

$${\left[\begin{array}[]{c}N_{i}(\lambda)\\ M_{i}(\lambda)\end{array}\right]=\left[\begin{array}[]{cc|c}A_{11}-\lambda E_{11}&A_{12}+B_{1}F_{2}-\lambda E_{12}&B_{1}W\\ 0&\gamma-\lambda\eta&B_{2}\\ \hline\cr C_{N,1}&C_{N,2}+D_{N}F_{2}&D_{N}W\\ C_{M,1}&C_{M,2}+D_{M}F_{2}&D_{M}W\end{array}\right]}.$$

$$\begin{array}[]{lcl}\widetilde{A}&\leftarrow&\left[\begin{smallmatrix}A_{11}&A_{12}+B_{1}F_{2}\\ 0&\gamma\end{smallmatrix}\right],\\ \\[-9.95845pt] \widetilde{E}&\leftarrow&\left[\begin{smallmatrix}E_{11}&E_{12}\\ 0&\eta\end{smallmatrix}\right],\\ \\[-9.95845pt] \widetilde{B}&\leftarrow&\left[\begin{smallmatrix}B_{1}W\\ B_{2}\end{smallmatrix}\right],\\ \widetilde{C}_{N}&\leftarrow&\big{[}\,C_{N,1}\;\;C_{N,2}+D_{N}F_{2}\,\big{]},\\ \widetilde{C}_{M}&\leftarrow&\big{[}\,C_{M,1}\;\;C_{M,2}+D_{M}F_{2}\,\big{]},\\ \widetilde{D}_{N}&\leftarrow&D_{N}W,\\ \widetilde{D}_{M}&\leftarrow&D_{M}W.\end{array}$$

$${\tilde{C}}_M=[0{\tilde{C}}_{M,2}],$$

$${\widetilde{A}=\left[\begin{array}[]{cc}{A}_{11}&{A}_{12}\\ 0&{A}_{22}\end{array}\right],\;\widetilde{E}=\left[\begin{array}[]{cc}{E}_{11}&{E}_{12}\\ 0&{E}_{22}\end{array}\right],\;\widetilde{B}=\left[\begin{array}[]{c}{B}_{1}\\ {B}_{2}\end{array}\right],}$$

$$F_2=-B_2^T(Y{E}_{22}^T{)}^{-1},W=I_m,$$

$$A_{22}Y{E}_{22}^T+E_{22}Y{A}_{22}^T-B_2{B}_2^T=0,$$

$$\begin{array}[]{l}F_{2}=-B_{2}^{T}(YA_{22}^{T})^{-1},\\ W^{T}\big{(}I+B_{2}^{T}(E_{22}YE_{22}^{T})^{-1}B_{2}\big{)}W=I,\end{array}$$

$$A_{22}Y{A}_{22}^T-B_2{B}_2^T=E_{22}Y{E}_{22}^T.$$

Full text

††journal: LAA

1 Introduction

Consider a disjunct partition of the complex plane $\mathds{C}$ as

[TABLE]

where we only assume $\infty\not\in\mathds{C}_{g}$. Any rational matrix $G(\lambda)$ can be expressed in a left fractional form

[TABLE]

or in a right fractional form

[TABLE]

where both the denominator factor $M(\lambda)$ and the numerator factor $N(\lambda)$ have only (finite) poles in $\mathds{C}_{g}$. Any method to compute a right factorization can be used to compute a left factorization by simply applying it to the transposed matrix $G^{T}(\lambda)$. Therefore, in this paper we only focus on right factorizations of the form (3).

Of special interest are the so-called coprime factorizations, where the factors satisfy additional conditions. A fractional representation of $G(\lambda)$ of the form (3), with $N(\lambda)$ and $M(\lambda)$ having only poles in $\mathds{C}_{g}$, is a right coprime factorization (RCF) over $\mathds{C}_{g}$, if there exist $U(\lambda)$ and $V(\lambda)$ with poles only in $\mathds{C}_{g}$ which satisfy $U(\lambda)M(\lambda)+V(\lambda)N(\lambda)=I\,.$ Coprime factorizations with minimum-degree denominators satisfy additionally the condition that the number of poles of the resulting $M(\lambda)$ is equal to the number of poles of $G(\lambda)$ in $\mathds{C}_{b}$. When counting poles, both finite and infinite poles are counted with their multiplicities. When determining minimum-degree coprime factorizations over $\mathds{C}_{g}$, all poles of $M(\lambda)$ can be arbitrarily chosen from $\mathds{C}_{g}$. The coprime factorizations over a “good” domain of poles $\mathds{C}_{g}$, also including factorizations with special properties of the denominator factor (e.g., stability, innerness), are important in extending various controller synthesis methods [1] to more general systems or in solving synthesis problems of fault detection filters [2] in the most general setting.

If $G(\lambda)$ represents the transfer-function matrix (TFM) of a (generalized) linear time-invariant system, then the variable $\lambda$ stands either for the complex variable $s$ in the Laplace transform, in the case of a continuous-time system, or for the complex variable $z$ in the $Z$-transform, in the case of a discrete-time system. In this case, $G(\lambda)$ has real coefficients (i.e., all entries are ratios of polynomials with real coefficients) and $\mathds{C}_{g}$ and $\mathds{C}_{b}$ are chosen symmetrically located with respect to the real axis, with $\mathds{C}_{g}$ having at least one point on the real axis. This guarantees the existence of real factorizations, where both factors $N(\lambda)$ and $M(\lambda)$ have real coefficients. Typical choices for $\mathds{C}_{g}$ are the open left half-plane of $\mathds{C}$ for a continuous-time system, or the open unit disc centered in the origin for a discrete-time system. In these cases, the corresponding factorizations are called proper and stable RCFs. An important class of such factorizations is the class of stable and proper RCF with inner denominators, where $M(\lambda)$ is additionally an inner TFM satisfying the all-pass property $M^{\sim}(\lambda)M(\lambda)=I$. Recall that, the conjugate TFM $M^{\sim}(\lambda)$ is defined in a continuous-time setting as $M^{\sim}(s)=M^{T}(-s)$, while in a discrete-time setting $M^{\sim}(z)=M^{T}(1/z)$. A TFM $M(\lambda)$ which satisfies only $M^{\sim}(\lambda)M(\lambda)=I$ is called all-pass. Therefore, an inner TFM is a stable all-pass TFM.

For the computation of right coprime factorizations with minimum-degree denominators a recursive pole-dislocation-based approach can be employed. Inspired by the ideas of Belevitch [3], such a method has been developed in [4] and later refined in [5]. The basic approach can be simply described as an $n_{b}$-step procedure, where $n_{b}$ is the number of poles of $G(\lambda)$ in $\mathds{C}_{b}$. To compute the RCF of $G(\lambda)$, it is possible to find a sequence of $n_{b}$ nonsingular rational matrices $\widetilde{M}_{i}(\lambda)$, $i=1,\ldots,n_{b}$, each of McMillan degree 1, with one (arbitrary) pole in $\mathds{C}_{g}$ and one (fixed) zero in $\mathds{C}_{b}$, such that the sequences $N_{i}(\lambda):=N_{i-1}(\lambda)\widetilde{M}_{i}(\lambda)$ and $M_{i}(\lambda):=M_{i-1}(\lambda)\widetilde{M}_{i}(\lambda)$ for $i=1,\ldots,n_{b}$, initialized with $N_{0}(\lambda)=G(\lambda)$ and $M_{0}(\lambda)=I$, generate the factors $N(\lambda):=N_{n_{b}}(\lambda)$ and $M(\lambda):=M_{n_{b}}(\lambda)$ of the RCF (3). The zero of $\widetilde{M}_{i}(\lambda)$ is chosen to cancel with a pole of $N_{i-1}(\lambda)$ lying in $\mathds{C}_{b}$, such that after $n_{b}$ steps, all poles of $G(\lambda)$ are canceled and dislocated to values in $\mathds{C}_{g}$. This approach can be also employed when additionally imposing that all elementary factors $\widetilde{M}_{i}(\lambda)$ are inner. In the case of complex poles, the above technique leads, in general, to complex factorizations. Therefore, to obtain real factorizations, for each complex conjugate pair of poles a second-degree real factor can be used to cancel simultaneously both poles. Second-degree factors may also be necessary when dislocating a pair of real poles into a pair of complex conjugate poles. General formulas for constructing first- and second-degree factors are given in [4]. In the case of an improper $G(\lambda)$, similar first-degree formulas can be devised for canceling a single infinite pole.

In this paper, we describe general descriptor state-space realization based versions of the above recursive approach, which eliminate the delicate computations involving the manipulation of rational matrices. The recursive pole dislocation techniques underlying the new general algorithms have been developed in the spirit of the approach described in [5] and extend the methods proposed in [6] by addressing the dislocation of infinite poles. Alternative, non-iterative approaches to compute coprime factorizations with inner denominators have been proposed in [7, 8] and involve the solution of generalized Lyapunov equations.

2 Descriptor system based computations

In this section we present some basic notions on descriptor systems, and describe the main reduction and updating techniques employed in the paper.

2.1 Preliminaries on descriptor systems

Let $G(\lambda)$ be a $p\times m$ real rational matrix and consider an $n$-th order descriptor system realization $(A-\lambda E,B,C,D)$, with $A-\lambda E$ an $n\times n$ regular pencil (i.e., $\det(A-\lambda E)\not\equiv 0$), which satisfies

[TABLE]

We will also use the equivalent notation for the TFM (4)

[TABLE]

If $Q,Z\in\mathds{R}^{n\times n}$ are invertible matrices, then the two realizations $(A-\lambda E,B,C,D)$ and $(\widetilde{A}-\lambda\widetilde{E},\widetilde{B},\widetilde{C},D)$, whose matrices are related by a similarity transformation of the form

[TABLE]

have the same TFM $G(\lambda)$.

We recall from [9, 10] some basic notions related to descriptor system realizations. A realization $(A-\lambda E,B,C,D)$ is minimal if it is controllable, observable and has no non-dynamic modes. The poles of $G(\lambda)$ are related to $\Lambda(A-\lambda E)$, the eigenvalues of the pencil $A-\lambda E$ (also known as the generalized eigenvalues of the pair $(A,E)$). For a minimal realization, the finite poles of $G(\lambda)$ are the finite eigenvalues of $A-\lambda E$, while the multiplicities of the infinite eigenvalues of $A-\lambda E$ are in excess one to the multiplicities of infinite poles. A finite eigenvalue $\lambda_{f}$ is controllable if $\mathop{\mathrm{rank}}\,[\,A-\lambda_{f}E\;B\,]=n$, otherwise it is uncontrollable. Similarly, a finite eigenvalue $\lambda_{f}$ is observable if $\mathop{\mathrm{rank}}\left[\begin{smallmatrix}A-\lambda_{f}E\\ C\end{smallmatrix}\right]=n$, otherwise it is unobservable. Infinite controllability requires that $\mathop{\mathrm{rank}}\,[\,E\;B\,]=n$, while infinite observability requires that $\mathop{\mathrm{rank}}\left[\begin{smallmatrix}E\\ C\end{smallmatrix}\right]=n$. The lack of non-dynamic modes can be equivalently expressed as $A\ker(E)\subseteq\mathop{\mathrm{Im}}(E)$.

Consider a given disjunct partition of the complex plane as in (1). We say the descriptor system (5) is proper $\mathds{C}_{g}$-stable if all finite eigenvalues of $A-\lambda E$ belong to $\mathds{C}_{g}$ and all infinite eigenvalues of $A-\lambda E$ are simple. The descriptor system (5) (or equivalently the pair $(A-\lambda E,B)$) is $\mathds{C}_{b}$-stabilizable if $\mathop{\mathrm{rank}}[\,A-\lambda E\;B\,]=n$ for all finite $\lambda\in\mathds{C}_{b}$. The descriptor system (5) (or equivalently the pair $(A-\lambda E,C)$) is $\mathds{C}_{b}$-detectable if $\mathop{\mathrm{rank}}\left[\begin{smallmatrix}A-\lambda E\\ C\end{smallmatrix}\right]=n$ for all finite $\lambda\in\mathds{C}_{b}$.

2.2 Reduction of the pole pencil to a special form

Using two orthogonal transformation matrices $Q$ and $Z$, it is possible to reduce a regular pencil $A-\lambda E$ to a specially ordered generalized real Schur form (GRSF) $\widetilde{A}-\lambda\widetilde{E}:=Q(A-\lambda E)Z$, with the eigenvalues of the component diagonal blocks split in four diagonal blocks:

[TABLE]

where: $(i)$ $A_{\infty}$ is an $(n-r)\times(n-r)$ invertible (upper triangular) matrix, with $r=\mathop{\mathrm{rank}}E$; the pair $(A_{\infty},0)$ contains all simple infinite eigenvalues of $A-\lambda E$ corresponding to first-order eigenvectors; $(ii)$ $A_{g}$ and $E_{g}$ are $n_{g}\times n_{g}$ matrices, such that the pair $(A_{g},E_{g})$ is in a GRSF (i.e., $A_{g}$ upper quasi-triangular and $E_{g}$ upper triangular) and $\Lambda(A_{g}-\lambda E_{g})\subset\mathds{C}_{g}$; $(iii)$ $A_{f,b}$ and $E_{f,b}$ are $n_{b}^{f}\times n_{b}^{f}$ matrices, such that the pair $(A_{f,b},E_{f,b})$ is in a GRSF and $\Lambda(A_{f,b}-\lambda E_{f,b})$ contains the finite eigenvalues lying in $\mathds{C}_{b}$; and $(iv)$ $A_{\infty,b}$ and $E_{\infty,b}$ are $n_{b}^{\infty}\times n_{b}^{\infty}$ upper triangular matrices and $\Lambda(A_{\infty,b}-\lambda E_{\infty,b})$ contains the rest of infinite eigenvalues. The matrices $E_{g}$, $E_{f,b}$, $A_{\infty,b}$ are invertible, while $E_{\infty,b}$ is nilpotent.

The reduction of a pair $(A,E)$ to the specially ordered GRSF (7) can be performed in three main steps. In the first step, we isolate the infinite eigenvalues having first order eigenvectors using two orthogonal transformation matrices $Q_{1}$ and $Z_{1}$ such that

[TABLE]

where the regularity of $A-\lambda E$ ensures that $A_{\infty}$ is nonsingular and the pencil $A_{22}-\lambda E_{22}$ is regular. This step involves the compression of the columns of $E$ using an orthogonal transformation matrix $Z_{1}$ such that $EZ_{1}=[\,0\;E_{2}\,]$, with $E_{2}$ full column rank $r=\mathop{\mathrm{rank}}E$, and the conformably partitioned $AZ_{1}=[\,A_{1}\;A_{2}\,]$, with $A_{1}$ having full column rank $n-r$. Typically, such a column compression of $E$ is performed using the rank-revealing QR decomposition with column pivoting applied to $E^{T}$ [11]. Then, an orthogonal matrix $Q_{1}$ is determined such that $Q_{1}A_{1}=\left[\begin{smallmatrix}A_{\infty}\\ 0\end{smallmatrix}\right]$, with $A_{\infty}$ an $(n-r)\times(n-r)$ invertible upper triangular matrix. This computation can be done using the standard QR decomposition of a full column rank matrix [11].

At the second step, the pair $(A_{22},E_{22})$ is reduced, using two orthogonal transformation matrices $Q_{2}$ and $Z_{2}$, to the form

[TABLE]

where the pair $(A_{\infty,b},E_{\infty,b})$ has only infinite generalized eigenvalues, $A_{\infty,b}\!-\!\lambda E_{\infty,b}$ is upper triangular, with $E_{\infty,b}$ nilpotent, and all eigenvalues of $A_{f}\!-\!\lambda E_{f}$ are finite. The infinite-finite splitting of eigenvalues can be reliably computed by applying the orthogonal reduction algorithm of [12] to the transposed pencil $A_{22}^{T}-\lambda E^{T}_{22}$. This algorithm determines the orthogonal transformation matrices $Z_{2}^{T}$ and $Q_{2}^{T}$, such that the reduced pencil $Z_{2}^{T}(A_{22}^{T}-\lambda E^{T}_{22})Q_{2}^{T}$ is in a $2\times 2$ block upper triangular form, where the leading block contains the infinite eigenvalues, with the matrices in an upper triangular (staircase) form, while the trailing block contains the finite eigenvalues. The form (8) is obtained by pertransposing the resulting pair,222Pertransposing a square matrix $M$ is the operation to form $PM^{T}\!P$, where $P$ is a permutation matrix with ones down the secondary diagonal. using a permutation matrix $P_{2}$ of appropriate size. This operation interchanges the order of the infinite and finite blocks and preserves the upper triangular shape of the blocks containing the infinite eigenvalues, which form the resulting pencil $A_{\infty,b}\!-\!\lambda E_{\infty,b}$.

At the third step, we reduce the pair $(A_{f},E_{f})$ to an ordered GRSF, using two orthogonal transformation matrices $Q_{3}$ and $Z_{3}$, such that

[TABLE]

where $\Lambda(A_{g}-\lambda E_{g})\subset\mathds{C}_{g}$ and $\Lambda(A_{f,b}-\lambda E_{f,b})\subset\mathds{C}_{b}$. Suitable algorithms (e.g., the QZ algorithm) to reduce a matrix pair to the GRSF and to reorder the eigenvalues by permuting the $1\times 1$ and $2\times 2$ diagonal blocks are described in [11].

The following procedure summarizes the main computational steps of the reduction of a regular pair $(A,E)$, to a pair $(\widetilde{A},\widetilde{E})=(QAZ,QEZ)$ in the specially ordered GRSF (7), by employing the orthogonal transformation matrices $Q$ and $Z$.

Procedure GSORSF.

0:   Compute an orthogonal $Z_{1}$ such that $EZ_{1}=[\,0\;E_{2}\,]$, with $E_{2}$ having full column rank $r=\mathop{\mathrm{rank}}E$; compute the conformably partitioned $AZ_{1}=[\,A_{1}\;A_{2}\,]$, with $A_{1}$ having full column rank $n-r$.

0:   Compute an orthogonal $Q_{1}$ such that $Q_{1}A_{1}=\left[\begin{smallmatrix}A_{\infty}\\ 0\end{smallmatrix}\right]$, with $A_{\infty}$ an $(n-r)\times(n-r)$ invertible upper triangular matrix; compute the conformably partitioned matrices

[TABLE]

0:   Apply the staircase algorithm of [12] to the transposed pencil $A_{22}^{T}-\lambda E_{22}^{T}$ to obtain the orthogonal transformation matrices $Z_{2}^{T}$ and $Q_{2}^{T}$ such that

[TABLE]

where $P_{2}$ is a permutation matrix of appropriate size with ones down on the secondary diagonal, the pair $(A_{\infty,b},E_{\infty,b})$ has only infinite generalized eigenvalues, $A_{\infty,b}\!-\!\lambda E_{\infty,b}$ is upper triangular, with $E_{\infty,b}$ nilpotent, and all eigenvalues of $A_{f}\!-\!\lambda E_{f}$ are finite.

0:   Compute orthogonal $Q_{3}$ and $Z_{3}$ such that

[TABLE]

is in a GRSF, where $\Lambda(A_{g}-\lambda E_{g})\subset\mathds{C}_{g}$ and $\Lambda(A_{f,b}-\lambda E_{f,b})\subset\mathds{C}_{b}$.

0:   Set $Q=\mathop{\mathrm{diag}}(I_{n-r},Q_{3},I)\mathop{\mathrm{diag}}(I_{n-r},P_{2}Q_{2})Q_{1}$, $Z=Z_{1}\mathop{\mathrm{diag}}(I_{n-r},Z_{2}P_{2})\mathop{\mathrm{diag}}(I_{n-r},Z_{3},I)$ and define $\widetilde{A}=QAZ$ and $\widetilde{E}=QEZ$ from the reduced pencil (7).

The computations performed at Steps 1, 2 and 4 of Procedure GSORSF rely on standard numerically stable algorithms (described, for example, in [11]), as the rank revealing QR decomposition with column pivoting or the QZ algorithm to compute and reorder the GRSF of a matrix pair. The numerical complexity of these steps is $\mathcal{O}(n^{3})$. The finite-infinite separation performed at Step 3, using the numerically stable staircase algorithm of [12], employs a rank revealing QR-decomposition-based staircase reduction technique. The performed orthogonal transformations consist of sequences of Householder or Givens transformations, which are not accumulated but directly applied to the involved submatrices. Using this reduction technique (i.e., without explicitly accumulating the performed orthogonal transformations employed for rank determinations), guarantees the $\mathcal{O}(n^{3})$ computational complexity of this step. Therefore, the overall computational complexity of the Procedure GSORSF is $\mathcal{O}(n^{3})$ as well. The reliability of rank decisions can be improved, using singular value decomposition (SVD)-based rank decisions instead QR-decomposition-based rank decisions. Such an approach has been proposed, for example, in [13]. Due to the need to explicitly accumulate the performed orthogonal transformations, the worst-case complexity of this approach is $\mathcal{O}(n^{4})$. An alternative computational technique of complexity $\mathcal{O}(n^{3})$, with the same level of reliability of rank determinations, has been suggested in [12], and consists in combining QR decomposition based reductions (without column pivoting) with SVD-based rank decisions (performed only for small matrices).

2.3 Descriptor system based updating formulas

The denominator factor $M(\lambda)$ of the RCF can be represented in a product form as

[TABLE]

where each $\widetilde{M}_{i}(\lambda)$ ($i=1,\ldots,k$) is an invertible elementary proper factor having either a real pole or a pair of complex conjugate poles. The recursive computational procedure can be formalized as $k$ successive applications of the updating formula

[TABLE]

initialized with $N_{0}(\lambda)=G(\lambda)$ and $M_{0}(\lambda)=I_{m}$. Then, $N(\lambda)=N_{k}(\lambda)$ and $M(\lambda)=M_{k}(\lambda)$. By this approach, it is automatically achieved that the resulting $M(\lambda)$ has the least achievable McMillan degree $n_{b}$.

We can derive state-space formulas for the efficient implementation of the updating operations in (10). Assume $N_{i-1}(\lambda)$ and $M_{i-1}(\lambda)$ have the descriptor realizations

[TABLE]

where $\Lambda(A_{22}-\lambda E_{22})\subset\mathds{C}_{b}$. We assume that $A_{22}-\lambda E_{22}$ is a $1\times 1$ pencil in the case when $A_{22}-\lambda E_{22}$ has a real or an infinite eigenvalue, or is a $2\times 2$ pencil, in the case when $A_{22}-\lambda E_{22}$ has a pair of complex conjugate eigenvalues. This form automatically results if the pair $(\widetilde{A},\widetilde{E})$ is in the specially ordered GRSF (7). We discuss several cases which can be encountered when performing the updating.

If $B_{2}=0$, then the eigenvalue(s) of $A_{22}-\lambda E_{22}$ is (are) not controllable, and thus can be removed to obtain lower order realizations $N_{i}(\lambda):=N_{i-1}(\lambda)$ and $M_{i}(\lambda):=M_{i-1}(\lambda)$

[TABLE]

If $B_{2}\not=0$ and the pencil $A_{22}-\lambda E_{22}$ has finite eigenvalues (i.e., $E_{22}$ is invertible), then the pair $(A_{22}-\lambda E_{22},B_{2})$ is (finite) controllable and there exists $F_{2}$ such that the eigenvalues of $A_{22}+B_{2}F_{2}-\lambda E_{22}$ can be placed in arbitrary locations in $\mathds{C}_{g}$. Assume that such an $F_{2}$ has been determined and define the elementary factor

[TABLE]

where $W$ is chosen to ensure the invertibility of $\widetilde{M}_{i}(\lambda)$. To compute stable and proper RCFs, the choice $W=I_{m}$ is always possible. However, alternative choices of $W$ are necessary to ensure, for example, that $\widetilde{M}_{i}(\lambda)$ is inner (see Section 3.2). It is easy to check that the updated factors $N_{i}(\lambda)$ and $M_{i}(\lambda)$ in (10) have the realizations

[TABLE]

If we denote $\widetilde{F}=[\,0\;F_{2}\,]$, then the above relations lead to the following updating formulas [6]:

[TABLE]

If $B_{2}\not=0$ and the $1\times 1$ pencil $A_{22}-\lambda E_{22}$ has an infinite eigenvalue (i.e., $E_{22}=0$), then we choose the elementary factor

[TABLE]

where $\eta\not=0$ and $\gamma/\eta$ is an arbitrary real pole in $\mathds{C}_{g}$, $W$ is a projection matrix chosen such $B_{2}W=0$ and $\mathop{\mathrm{rank}}\left[\begin{smallmatrix}B_{2}\\ W\end{smallmatrix}\right]=m$, and $F_{2}$ is chosen such that $B_{2}F_{2}=-A_{22}$ and $\mathop{\mathrm{rank}}[\,F_{2}\;W\,]=m$ (the rank conditions guarantee the invertibility of $\widetilde{M}_{i}(\lambda)$). Straightforward choices of $F_{2}$ and $W$ are, for example,

[TABLE]

By this choice of $\widetilde{M}_{i}(\lambda)$, we made the infinite eigenvalue in the realization of the updated factors $N_{i}(\lambda)$ and $M_{i}(\lambda)$ uncontrollable, and after its elimination, we obtain the realizations

[TABLE]

The above relations lead to the following updating formulas:

[TABLE]

The updating techniques relying on the formulas (14) and (17) ensure that, if the original pair $(\widetilde{A},\widetilde{E})$ was in a GRSF, then the updated pair will have a similar form, possibly with $\widetilde{A}-\lambda\widetilde{E}$ having a $2\times 2$ trailing block which corresponds to two real generalized eigenvalues (to recover the GRSF, such a block can be further split into two $1\times 1$ blocks using an orthogonal similarity transformation).

3 Computation of RCFs

In this section we present the computational procedures to compute proper $\mathds{C}_{g}$-stable RCFs and stable RCFs with inner denominators, and discuss the main numerical features of the proposed algorithms, as generality, numerical reliability and computational efficiency.

3.1 Computation of proper $\mathds{C}_{g}$-stable RCFs

To compute a proper $\mathds{C}_{g}$-stable RCF of an arbitrary rational TFM $G(\lambda)$ with a descriptor system realization $(A-\lambda E,B,C,D)$, we start by reducing the pair $(A,E)$ to the special GRSF (7) by performing a system similarity transformation as in (6) using orthogonal transformation matrices $Q$ and $Z$. Then, by employing the updating techniques described previously, we can dislocate the generalized eigenvalues of the trailing elementary $1\times 1$ or $2\times 2$ blocks to locations in $\mathds{C}_{g}$. Finally, the reordering of the diagonal blocks in the GRSF of the updated pair $(\widetilde{A},\widetilde{E})$, brings in the trailing position a new block, whose generalized eigenvalues (finite or infinite) lie in $\mathds{C}_{b}$. This eigenvalue dislocation process is repeated until all eigenvalues are moved into $\mathds{C}_{g}$.

The following procedure computes for an arbitrary $p\times m$ rational TFM $G(\lambda)$, with a descriptor system realization $(A-\lambda E,B,C,D)$, and for a given disjunct partition of the complex plane $\mathds{C}$ as in (1), a RCF $G(\lambda)=N(\lambda)M^{-1}(\lambda)$, with the resulting factors having proper $\mathds{C}_{g}$-stable descriptor realizations $N(\lambda)=(\widetilde{A}-\lambda\widetilde{E},\widetilde{B},\widetilde{C}_{N},\widetilde{D}_{N})$ and $M(\lambda)=(\widetilde{A}-\lambda\widetilde{E},\widetilde{B},\widetilde{C}_{M},\widetilde{D}_{M})$.

Procedure GRCF.

0:   Compute, using Procedure GSORSF, the orthogonal matrices $Q$ and $Z$ to reduce the pair $(A,E)$ to the specially ordered GRSF $(\widetilde{A},\widetilde{E})$ in (7). Compute $\widetilde{B}:=QB$, $\widetilde{C}_{N}:=CZ$, and set $\widetilde{C}_{M}=0$, $\widetilde{D}_{N}=D$, $\widetilde{D}_{M}=I_{m}$, $q:=r+n_{g}$, $n_{b}=n_{b}^{f}+n_{b}^{\infty}$.

0:   If $n_{b}=0$, Exit.

0:   Let $(A_{22},E_{22})$ be the last $k\times k$ diagonal blocks of the GRSF of $(\widetilde{A},\widetilde{E})$ (with $k\!=\!1$ or $k\!=\!2$) and let $B_{2}$ be the $k\times m$ matrix formed from the last $k$ rows of $\widetilde{B}$. If $B_{2}=0$, then remove the uncontrollable eigenvalues $\Lambda(A_{22}-\lambda E_{22})$ and redefine $\widetilde{A}$, $\widetilde{E}$, $\widetilde{B}$, $\widetilde{C}_{N}$, $\widetilde{C}_{M}$ according to (12); update $n\leftarrow n-k$, $n_{b}\leftarrow n_{b}-k$ and go to Step 2.

0:   If $E_{22}=0$, compute $F_{2}=-B_{2}^{T}(B_{2}B_{2}^{T})^{-1}A_{22}$ and $W=I_{m}-B_{2}^{T}(B_{2}B_{2}^{T})^{-1}B_{2}$. Choose $\eta=1$ and $\gamma\in\mathds{C}_{g}$ and update $\widetilde{A}$, $\widetilde{E}$, $\widetilde{B}$, $\widetilde{C}_{N}$, $\widetilde{D}_{N}$, $\widetilde{C}_{M}$ and $\widetilde{D}_{M}$ using (17).

0:   If $E_{22}\not=0$, determine $F_{2}$ such that $\Lambda(A_{22}+B_{2}F_{2}-\lambda E_{22})\subset\mathds{C}_{g}$. Set $\widetilde{F}=[\,0\;\;F_{2}\,]$, $W=I_{m}$ and update $\widetilde{A}$, $\widetilde{E}$, $\widetilde{B}$, $\widetilde{C}_{N}$, $\widetilde{D}_{N}$, $\widetilde{C}_{M}$ and $\widetilde{D}_{M}$ using (14).

0:   Compute the orthogonal matrices $\widetilde{Q}$ and $\widetilde{Z}$ to move the last blocks of $(\widetilde{A},\widetilde{E})$ to positions $(q+1,q+1)$ by interchanging the diagonal blocks of the GRSF. Compute $\widetilde{A}\leftarrow\widetilde{Q}\widetilde{A}\widetilde{Z}$, $\widetilde{E}\leftarrow\widetilde{Q}\widetilde{E}\widetilde{Z}$, $\widetilde{B}\leftarrow\widetilde{Q}\widetilde{B}$, $\widetilde{C}_{N}\leftarrow\widetilde{C}_{N}\widetilde{Z}$, $\widetilde{C}_{M}\leftarrow\widetilde{C}_{M}\widetilde{Z}$. Put $q\leftarrow q+k$, $n_{b}\leftarrow n_{b}-k$ and go to Step 2.

With the special GRSF (7) computed at Step 1, the procedure executes at the beginning the Steps 3, 4 and 6 repeatedly, until all controllable infinite eigenvalues are dislocated to finite locations. In the reordering of diagonal blocks performed at Step 6 to recover the special GRSF (7), it is possible to exploit the structure of the blocks containing the infinite generalized eigenvalues of the pair $(\widetilde{A},\widetilde{E})$. Consequently, the swapping of two $1\times 1$ diagonal blocks, where the first one contains an infinite and the second one contains a finite eigenvalue, can be reliably performed by explicitly restoring the zero diagonal elements in $\widetilde{E}$ (these elements are usually blurred by roundoff errors incurred by swapping). This part of the algorithm ensures that the resulting factors are proper and can be interpreted as a more elaborated version of both of the Algorithms PRRCF1 and PRRCF2 presented in [6].

After the completion of the dislocation of infinite eigenvalues, the procedure executes the Steps 3, 5 and 6 repeatedly, until all controllable finite eigenvalues in $\mathds{C}_{b}$ are dislocated to new locations in $\mathds{C}_{g}$. This computation corresponds to Algorithm GRCF-P in [6] and can be interpreted as an extension of the GRSF-based method for pole assignment [14].

The Procedure GRCF is completely general, being applicable regardless the original descriptor realization of $G(\lambda)$ is $\mathds{C}_{b}$-stabilizable or not, is infinite controllable or not. The resulting pair $(\widetilde{A},\widetilde{E})$ is in a particular GRSF with $n-r$ simple infinite eigenvalues in the leading $n-r$ positions (no such block exists if $E$ is invertible). Thus, the elimination of non-dynamic (simple) modes can be done by applying standard state residualization formulas.

A minimal realization of the denominator $M(\lambda)$ can be easily determined. The resulting $\widetilde{C}_{M}$ has always the form

[TABLE]

where the number of columns of $\widetilde{C}_{M,2}$ is equal to the number of controllable generalized eigenvalues of the pair $(A,E)$ lying in $\mathds{C}_{b}$. By partitioning accordingly the resulting $\widetilde{E}$, $\widetilde{A}$ and $\widetilde{B}$

[TABLE]

we obtain $({A}_{22}-\lambda{E}_{22},{B}_{2},\widetilde{C}_{M,2},\widetilde{D}_{M})$, which is a minimal descriptor system realization of $M(\lambda)$. The order of the minimal realization of $M(\lambda)$ has the least possible McMillan degree if all eigenvalues of $A-\lambda E$ in $\mathds{C}_{b}$ are observable (i.e, the pair $(A-\lambda E,C)$ is $\mathds{C}_{b}$-detectable).

3.2 Computation of stable RCFs with inner denominator

For the computation of stable RCFs with inner denominators, a recursive procedure, similar to that of previous section, can be developed. In this case, we use the partition of the complex plane with $\mathds{C}_{g}=\mathds{C}_{s}$ and $\mathds{C}_{b}=\mathds{C}\setminus\mathds{C}_{s}$, where $\mathds{C}_{s}$ is the stability domain for the poles. Recall, that $\mathds{C}_{s}$ is the open left half-plane for a continuous-time system, or the interior of the unit circle centered in the origin, for a discrete-time system. For a given rational matrix $G(\lambda)$, a necessary and sufficient condition for the existence of a stable RCF with inner denominator is that $G(\lambda)$ has no poles in $\partial\mathds{C}_{s}$, the boundary of $\mathds{C}_{s}$. In the continuous-time case, this means that $G(s)$ is proper and has no poles on the imaginary axis, while in the discrete-time case $G(z)$ has no poles on the unit circle centered in the origin, but for the sake of generality, $G(z)$ can be improper (i.e., having infinite poles).

For a descriptor system realization $(A-\lambda E,B,C,D)$ of $G(\lambda)$, sufficient conditions for the existence of a stable RCF with inner denominator can be formulated in terms of the eigenvalues of the pole pencil $A-\lambda E$. In the continuous-time case, we merely require that the pencil $A-sE$ has no finite controllable eigenvalues on the imaginary axis and no higher order controllable infinite eigenvalues. In the discrete-time case, the pencil $A-zE$ has no controllable eigenvalues on the unit circle centered in the origin, but may have arbitrary number of infinite eigenvalues.

In the recursive procedure to compute the RCF with inner denominator, elementary inner factors are employed to dislocate one or a pair of complex conjugate poles at each iteration. Let us assume that at the $i$-th step, we have the matrices $A_{22}$, $E_{22}$ and $B_{2}$ as defined in (11). We encounter two cases when determining the $i$-th elementary factor $\widetilde{M}_{i}(\lambda)$ in (9).

If $B_{2}\not=0$ and $E_{22}$ invertible, then we choose $\widetilde{M}_{i}(\lambda)$ of the form (13), with $F_{2}$ and $W$ determined as follows (see [6, Fact 5]): for a continuous-time system

[TABLE]

where $Y$ is the solution of the generalized continuous-time Lyapunov equation (GCLE)

[TABLE]

while for a discrete-time system $F_{2}$ and $W$ are computed from

[TABLE]

where $Y$ is the solution of the generalized discrete-time Lyapunov equation (GDLE)

[TABLE]

The above choice ensures that the poles of $\widetilde{M}_{i}(\lambda)$ are the reflected generalized eigenvalues of the pair $(A_{22},E_{22})$ with respect to the imaginary axis, in the continuous-time case, or with respect to the unit circle in the origin, in the discrete-time case. For this case, the updating formulas (14) are used.

The second case may appear only for an improper discrete-time system, for which $B_{2}\not=0$, $A_{22}\not=0$, and $E_{22}=0$. The corresponding infinite pole is reflected to the origin, by choosing $\widetilde{M}_{i}(z)$ of the form (15), with $\gamma=0$, $\eta=-A_{22}$, and $F_{2}$ and $W$ determined as in (16). For this case, the updating formulas (17) are used. The formulas (16) underlying the computation of the elementary inner factor $\widetilde{M}_{i}(z)$ have been derived using the dual results of [8, Theorem 3.2].

The following procedure computes for a $p\times m$ rational TFM $G(\lambda)$, with a descriptor system realization $(A-\lambda E,B,C,D)$, a RCF with inner denominator $G(\lambda)=N(\lambda)M^{-1}(\lambda)$, with the resulting factors having proper $\mathds{C}_{s}$-stable descriptor realizations $N(\lambda)=(\widetilde{A}-\lambda\widetilde{E},\widetilde{B},\widetilde{C}_{N},\widetilde{D}_{N})$ and $M(\lambda)=(\widetilde{A}-\lambda\widetilde{E},\widetilde{B},\widetilde{C}_{M},\widetilde{D}_{M})$.

Procedure GRCFID.

0:   Compute, using Procedure GSORSF, the orthogonal matrices $Q$ and $Z$ to reduce the pair $(A,E)$ to the specially ordered GRSF $(\widetilde{A},\widetilde{E})$ in (7). Compute $\widetilde{B}:=QB$, $\widetilde{C}_{N}:=CZ$, and set $\widetilde{C}_{M}=0$, $\widetilde{D}_{N}=D$, $\widetilde{D}_{M}=I_{m}$, $q:=r+n_{g}$, $n_{b}=n_{b}^{f}+n_{b}^{\infty}$.

0:   If $n_{b}=0$, Exit.

0:   Let $(A_{22},E_{22})$ be the last $k\times k$ diagonal blocks of the GRSF of $(\widetilde{A},\widetilde{E})$ (with $k\!=\!1$ or $k\!=\!2$) and let $B_{2}$ be the $k\times m$ matrix formed from the last $k$ rows of $\widetilde{B}$. If $B_{2}=0$, then remove the uncontrollable eigenvalues $\Lambda(A_{22}-\lambda E_{22})$ and redefine $\widetilde{A}$, $\widetilde{E}$, $\widetilde{B}$, $\widetilde{C}_{N}$, $\widetilde{C}_{M}$ according to (12); update $n\leftarrow n-k$, $n_{b}\leftarrow n_{b}-k$ and go to Step 2.

0:   If $E_{22}=0$, then if the descriptor system is continuous-time, Exit (no solution exists); else, set $\gamma=0$, $\eta=-A_{22}$, compute $F_{2}$ and $W$ according to (16), and update $\widetilde{A}$, $\widetilde{E}$, $\widetilde{B}$, $\widetilde{C}_{N}$, $\widetilde{D}_{N}$, $\widetilde{C}_{M}$ and $\widetilde{D}_{M}$ using (17).

0:   If $E_{22}\not=0$, then if $\Lambda(A_{22}-\lambda E_{22})\subset\partial\mathds{C}_{s}$, Exit (no solution exists); else, determine $F_{2}$ and $W$ as in (20) for a continuous-time system, or as in (22) for a discrete-time system, set $\widetilde{F}=[\,0\;F_{2}\,]$, and update $\widetilde{A}$, $\widetilde{E}$, $\widetilde{B}$, $\widetilde{C}_{N}$, $\widetilde{D}_{N}$, $\widetilde{C}_{M}$ and $\widetilde{D}_{M}$ using (14).

0:   Compute the orthogonal matrices $\widetilde{Q}$ and $\widetilde{Z}$ to move the last blocks of $(\widetilde{A},\widetilde{E})$ to positions $(q+1,q+1)$ by interchanging the diagonal blocks of the GRSF. Compute $\widetilde{A}\leftarrow\widetilde{Q}\widetilde{A}\widetilde{Z}$, $\widetilde{E}\leftarrow\widetilde{Q}\widetilde{E}\widetilde{Z}$, $\widetilde{B}\leftarrow\widetilde{Q}\widetilde{B}$, $\widetilde{C}_{N}\leftarrow\widetilde{C}_{N}\widetilde{Z}$, $\widetilde{C}_{M}\leftarrow\widetilde{C}_{M}\widetilde{Z}$. Put $q\leftarrow q+k$, $n_{b}\leftarrow n_{b}-k$ and go to Step 2.

The Procedure GRCFID is completely general and is able to compute a stable RCF with inner denominator, provided the existence conditions formulated in terms of the poles of $G(\lambda)$ are fulfilled. Note that uncontrollable eigenvalues in $\mathds{C}_{b}$ or at infinity are automatically deflated. In contrast, the alternative, non-iterative approaches proposed in [7, 8], are only applicable to stabilizable descriptor realizations. The resulting pair $(\widetilde{A},\widetilde{E})$ is in a particular GRSF with $n-r$ simple infinite eigenvalues in the leading $n-r$ positions, which allows the immediate elimination of non-dynamic (simple) modes. A minimal realization of the denominator $M(\lambda)$ can be similarly determined as in the case of Procedure GRCF and has the least McMillan degree, provided the descriptor realization is $\mathds{C}_{b}$-detectable.

Remark 1. In the particular case, when $G(\lambda)$ is a square all-pass and anti-stable TFM (i.e., with all poles in $\mathds{C}\setminus\mathds{C}_{s}$ and all zeros in $\mathds{C}_{s}$), a RCF with inner denominator $G(\lambda)=N(\lambda)M^{-1}(\lambda)$ has the factors $N(\lambda)=I$ and $M(\lambda)=G^{\sim}(\lambda)$. It follows, that all eigenvalues in the resulting realization of $N(\lambda)$ computed by the Procedure GRCFID will be unobservable. More generally, if the given TFM $G(\lambda)$ contains an unstable all-pass factor with, say $k$ unstable poles, then the resulting descriptor system realization of $N(\lambda)$, computed by the Procedure GRCFID, will contain $k$ unobservable eigenvalues, which are precisely the stable zeros of the all-pass factor. See Example 2 presented in Section 5, with a discrete-time TFM $G(z)$ containing an unstable all-pass factor.

3.3 Numerical aspects

As already mentioned, the Procedures GRCF and GRCFID are completely general, being able to compute RCFs of rational matrices, independently of the properties of the underlying descriptor system realization. For the computation of RCFs with minimum degree denominator factors, the $\mathds{C}_{b}$-detectability condition is sufficient to guarantee the least achievable McMillan degrees of the denominators. The right coprimeness of the computed factors follows from the full (column) rank property of $\left[\begin{smallmatrix}N(\lambda)\\ M(\lambda)\end{smallmatrix}\right]$ for all $\lambda\in\mathds{C}_{g}$ [15].

The Procedure GRCF is essentially a recursive pole assignment algorithm which exploits and maintains the GRSF of the pair $(\widetilde{A},\widetilde{E})$. It represents a specialization of the generalized pole assignment algorithm of [14] for descriptor systems or the pole assignment algorithm of [16] for standard systems (i.e., $E=I$). These algorithms are generally considered as satisfactory computational methods, as long as the partial feedback gains $F_{2}$ computed at Steps 4 or 5, have moderate sizes, as—for example, $\|F_{2}\|\leq\kappa\|A\|/\|B\|$, with say $\kappa<100$. A careful implementation of the computation of $F_{2}$ in conjunction with optimal choices of the poles to be assigned can significantly contribute to the reduction of the size of $F_{2}$. Unfortunately, the above restriction on the partial feedback gains cannot be always fulfilled if large gains are necessary to move all poles into $\mathds{C}_{g}$. This may occur either if the poles in the “bad” region $\mathds{C}_{b}$ are too far from those in the “good” region $\mathds{C}_{g}$ or if these poles are weekly controllable (e.g., small value of $\|B_{2}\|$). Nevertheless, the Procedure GRCF can be still considered a numerically reliable algorithm, since the above norm condition on $F_{2}$ can be easily checked at each computational step and, therefore, the potential loss of numerical reliability can be easily detected.

The numerical properties of Procedure GRCFID are similar to those of Procedure GRCF, as long as the partial feedback gains $\|F_{2}\|$ at Steps 4 and 5 are reasonably small. However, this condition for numerical reliability may not always be fulfilled due to the lack of any freedom in assigning the poles. Recall that the unstable poles are reflected into symmetrical positions with respect to $\partial\mathds{C}_{s}$, and this may occasionally require large gains.

The main computation burden in the proposed algorithms is the computation of the specially ordered GRSF in (7) at Step 1 and the preservation of this form using eigenvalue reordering techniques at Step 6. With the use of standard algorithms (i.e., the QR decomposition based on the Householder method, the QZ algorithm to compute and reorder the GRSF), see [11], and of the staircase algorithm of [12], the overall numerical complexity is $\mathcal{O}(n^{3})$.

4 Implementation issues

In this section we discuss several implementation issues, which are instrumental for an efficient and robust implementation of the proposed computational procedures. The discussed issues are related to the choice of the poles to be assigned at each iteration in the Procedure GRCF, the use of enhanced accuracy (square-root) techniques to compute the elementary inner factors in Procedure GRCFID, the elimination of the non-dynamic modes of the resulting factors and the enforcement of the upper quasi-triangular shape of the resulting system matrices when computing left coprime factorizations.

4.1 Pole selection and assignment

The selection of appropriate poles to be assigned at each iteration can significantly influence the overall numerical behavior of the Procedure GRCF. The strategies to select poles depend on the concrete definition of the “good” region $\mathds{C}_{g}$ (and also of its complement $\mathds{C}_{b}=\mathds{C}\setminus\mathds{C}_{g}$), and target either a stabilization oriented factorization or a pole allocation oriented factorization.

4.1.1 Stabilization oriented factorization

A relevant measure to characterize the desired dynamics of the factors is the stability degree, say $\alpha$, of the poles of the resulting factors. For a continuous-time system, $\alpha<0$ and $\alpha$ represents the maximum admissible real part of poles. Consequently, $\mathds{C}_{g}$ can be defined as

[TABLE]

For a discrete-time system, $0\leq\alpha<1$ and $\alpha$ represents the maximum admissible magnitude of poles. Therefore, $\mathds{C}_{g}$ can be defined as

[TABLE]

For a given stability degree $\alpha$, the assignment of poles of the elementary factors at Steps 4 and 5 of Procedure GRCF can be done in such a way, that partial feedback gains $F_{2}$ with reduced size result.

In the case of an $1\times 1$ block $A_{22}-\lambda E_{22}$ with a finite real eigenvalue, the pole to be assigned can always be chosen $\gamma=\alpha$. This choice leads at Step 5 to a least norm of $F_{2}$, which can be computed by solving

[TABLE]

for the least-norm solution $F_{2}$. If $B_{2}$ has the RQ decomposition

[TABLE]

where $\sigma\not=0$ and $V=[\,V_{1}\;V_{2}\,]$ is orthogonal, with $V_{1}$ an $m\times 1$ matrix, then $F_{2}$ can be explicitly computed as

[TABLE]

In the case of an $1\times 1$ block $A_{22}-\lambda E_{22}$ with an infinite eigenvalue (i.e., $E_{22}=0$), we can use the RQ decomposition (26) to determine the feedback gain $F_{2}$ and the projection matrix $W$ in (16) as

[TABLE]

In the case of a $2\times 2$ block corresponding to a pair of complex conjugate eigenvalues $\mu\pm\mathrm{i}\tau$, to reduce the size of the partial feedback $F_{2}$ it is desirable to perform a minimum shifting of these poles to their new locations. This can be achieved by choosing the poles to be assigned $\{\gamma_{1},\gamma_{2}\}$ as $\alpha\pm\mathrm{i}\tau$ in the continuous-time case and $\alpha(\mu\pm\mathrm{i}\tau)/\sqrt{\mu^{2}+\tau^{2}}$ in the discrete-time case. The computation of $F_{2}$ is done depending on the rank of $B_{2}$.

If $\mathop{\mathrm{rank}}B_{2}=2$ and $\Theta$ is $2\times 2$ matrix such that the eigenvalues of $E_{22}^{-1}\Theta$ are the eigenvalues to be assigned, then $F_{2}$ can be computed by solving

[TABLE]

for the least-norm solution $F_{2}$. This solution can be computed using the SVD of $B_{2}$

[TABLE]

where $\Sigma$ is a $2\times 2$ invertible diagonal matrix, and $U$ and $V$ are orthogonal matrices. If we partition $V$ as $V=[\,V_{1}\;V_{2}\,]$, with $V_{1}$ an $m\times 2$ matrix, then $F_{2}$ results as

[TABLE]

The choice of $\Theta$ is not unique and an optimal choice would be one for which the (Frobenius) norm of $F_{2}$ is minimum. For a pair of complex conjugate eigenvalues $\{\gamma_{1},\gamma_{2}\}$, the matrix $\Theta$ can be parameterized in terms of two parameters $\theta_{1}$ and $\theta_{2}$ as

[TABLE]

If we denote with $F_{2}(\theta_{1},\theta_{2})$ the expression (30) evaluated for the above $\Theta(\theta_{1},\theta_{2})$, then the optimal values of the free parameters $\theta_{1}$ and $\theta_{2}$ can be determined by minimizing $\|F_{2}(\theta_{1},\theta_{2})\|_{F}^{2}$.

If $\mathop{\mathrm{rank}}B_{2}=1$, then Procedure A of [16] can be applied to the standard controllable pair $(E_{22}^{-1}A_{22},E_{22}^{-1}B_{2})$ to determine $F_{2}$ which assigns the selected pair of complex conjugate poles. In this case, we compute the SVD of $E_{22}^{-1}B_{2}$ as

[TABLE]

where \Sigma=\left[\begin{array}[]{c}\sigma\\ 0\end{array}\right], is a $2\times 1$ matrix with $\sigma\not=0$, and $U$ and $V$ are orthogonal matrices. We compute

[TABLE]

where $\widetilde{\alpha}_{21}\not=0$ is guaranteed by the controllability of the pair $(A_{22}-\lambda E_{22},B_{2})$. If we partition $V$ as $V=[\,V_{1}\;V_{2}\,]$, with $V_{1}$ an $m\times 1$ matrix, then the partial feedback $F_{2}$ is obtained as

[TABLE]

where $\widetilde{F}=[\ \varphi_{1}\;\varphi_{2}\,]$, with $\varphi_{1}$ and $\varphi_{2}$ computed as

[TABLE]

4.1.2 Pole assignment oriented factorization

Let $\alpha$ be a desired stability degree for the poles of the resulting factors and define the “good” region $\mathds{C}_{g}$ as in (24) or (25), depending on the system type, continuous- or discrete-time, respectively. Simultaneously, let $\Gamma$ be a symmetric set of complex values such that $\Gamma\subset\mathds{C}_{s}$, which contains the desired poles to be assigned for the two factors. In this setting, the initial reduction of the pair $(A,E)$ to the special GRSF (7) ensures, that all eigenvalues of $A-\lambda E$ lying in $\mathds{C}_{g}$ are preserved in the descriptor system realizations of the resulting factors and only the eigenvalues lying in the “bad” region $\mathds{C}_{b}$ or at infinity are assigned to the values specified in $\Gamma$. We will assume that the number of specified poles in $\Gamma$ is greater than the sum of the number of controllable (higher order) infinite eigenvalues and the number of controllable finite eigenvalues in $\mathds{C}_{b}$. Otherwise, $\Gamma$ is extended “on the fly” with poles chosen from $\mathds{C}_{g}$, in accordance with the strategy used for the stabilization oriented factorization. In determining the partial feedback $F_{2}$ at Steps 4 and 5 of Procedure GRCF, we encounter several cases according to the size of the last block $A_{22}-\lambda E_{22}$ and the chosen poles to be assigned. After each successful partial pole assignment, $\Gamma$ is updated by removing the already assigned poles.

If $A_{22}-\lambda E_{22}$ is an $1\times 1$ block with a finite eigenvalue, then a real pole $\gamma\in\Gamma$ is selected, which is the nearest one to the eigenvalue of $A_{22}-\lambda E_{22}$. This pole is then assigned, using (26) and (27). In the case, when there is no real pole available in $\Gamma$, but $n_{b}>1$, then a $2\times 2$ block can be formed, either by adjoining an adjacent $1\times 1$ block or by interchanging the last two blocks of the GRSF to bring a $2\times 2$ block in the last position. Then, a pair of complex conjugate poles from $\Gamma$ can be assigned (see below). If $n_{b}=1$, then $\gamma=\alpha$ is assigned.

If $A_{22}-\lambda E_{22}$ has an infinite eigenvalue (i.e., $E_{22}=0$), then $\gamma$ is selected as the nearest real pole in $\Gamma$ to the boundary $\partial\mathds{C}_{s}$ and the formulas for $F_{2}$ and $W$ in (28) are used. If no real pole is available in $\Gamma$, then $\gamma=\alpha$ is chosen. It is possible to use more complicated schemes, as for example, adjoining two adjacent blocks with infinite eigenvalues, or a finite $1\times 1$ block and an infinite block, and assigning a pair of complex conjugate poles using suitable updating formulas (still to be developed). The chosen strategy to assign only real poles for all infinite poles appears to be the simplest one, and therefore, well suited for implementation purposes. A trivial workaround to assign a pair of complex conjugate poles to replace two infinite poles, is to assign first two real poles for the two last $1\times 1$ blocks, and then to update the resulting $2\times 2$ block by assigning a pair of complex conjugate poles (see below).

If $A_{22}-\lambda E_{22}$ is a $2\times 2$ block, then we choose the poles to be assigned $\{\gamma_{1},\gamma_{2}\}$, as the nearest ones in $\Gamma$ to the eigenvalues of $A_{22}-\lambda E_{22}$. If $A_{22}-\lambda E_{22}$ has complex conjugate eigenvalues, then we assign either two real poles or, preferably, a pair of complex conjugate poles, while if $A_{22}-\lambda E_{22}$ resulted by adjoining two $1\times 1$ diagonal blocks, then a pair of complex conjugate poles is assigned. Depending on the rank of the corresponding $B_{2}$, we employ for the computation of $F_{2}$ either the formula (30) if $\mathop{\mathrm{rank}}B_{2}=2$ or the formula (32) if $\mathop{\mathrm{rank}}B_{2}=1$.

If $A_{22}-\lambda E_{22}$ resulted by adjoining two $1\times 1$ diagonal blocks and $\mathop{\mathrm{rank}}B_{2}=1$, it may happen that the resulting pair $(A_{22}-\lambda E_{22},B_{2})$ is not controllable (because the real eigenvalue in the leading $1\times 1$ diagonal block of $A_{22}-\lambda E_{22}$ may be uncontrollable). The lack of controllability can be easily detected, because $\widetilde{\alpha}_{21}=0$ in (31). The uncontrollable eigenvalue $\widetilde{\alpha}_{22}$ can be eliminated by performing an orthogonal similarity transformation on the pair $(A_{22}-\lambda E_{22},B_{2})$ to obtain

[TABLE]

with $\widetilde{B}_{2}\not=0$. Here, $U$ is the orthogonal matrix from the SVD of $B_{2}$ in (29), while $V$ is the orthogonal matrix from the RQ decomposition of $UE_{22}$ (i.e., $UE_{22}V$ is upper triangular). After the separation of the uncontrollable $1\times 1$ block $\alpha_{22}-\lambda\eta_{22}$, the updating formulas (12) are used to remove this block.

4.2 Square-root based computation of elementary inner factors

In the Procedure GRCFID, the computation of each elementary inner factor $\widetilde{M}_{i}(\lambda)$ as in (13) involves the solution of either of the GCLE (21) or of the GDLE (23), from which the partial feedback $F_{2}$ and feedthrough matrix $W$ are determined. In each case, the solution $Y$ is positive definite, thus it can be expressed as $Y=SS^{T}$, where $S$ is an upper triangular matrix, also called (improperly) the square-root of $Y$. Methods to solve Lyapunov equations directly for the square-root factor of the solution have been proposed in [17] for $E_{22}=I$ and in [18] for a general $E_{22}$. These methods are provably more accurate than methods which computes $Y$, because they avoid explicitly forming the matrix product $B_{2}B_{2}^{T}$. These methods are applicable provided $\Lambda(A_{22}-\lambda E_{22})\subset\mathds{C}_{s}$ and $B_{2}B_{2}^{T}$ stays instead $-B_{2}B_{2}^{T}$ in (21) or (23). To cope with these restrictions, instead of solving (21), we can solve

[TABLE]

where $\Lambda(-A_{22}-\lambda E_{22})\subset\mathds{C}_{s}$, and instead solving (23), we can solve

[TABLE]

where $\Lambda(E_{22}-\lambda A_{22})\subset\mathds{C}_{s}$. Thus, both above Lyapunov equations can be solved directly for the square-root factor $S$.

The upper triangular shape of $S$ and the diagonal form of $E_{22}$ can be exploited in evaluating the partial feedback $F_{2}$ in (20) as

[TABLE]

Moreover, for an $1\times 1$ block, we can use that $A_{22}+B_{2}F_{2}=-A_{22}$. Similarly, $F_{2}$ in (22) can be computed as

[TABLE]

and for an $1\times 1$ block, we can use that $A_{22}+B_{2}F_{2}=E_{22}^{2}/A_{22}$. The feedthrough term $W$ in (22) can be computed as $W=R^{-1}$, with $R$ an upper triangular (Cholesky) factor satisfying

[TABLE]

where $X=S^{-1}E_{22}^{-1}B_{2}$ is a $k\times m$ matrix, with $k=1$ or $k=2$. To compute $R$, the initial Cholesky factor $I_{m}$ can be updated, by performing, $k$ times rank-1 changes, thus avoiding to form explicitly the product $X^{T}X$ (see [11][Section 6.5.1] for suitable QR factorization based updating techniques).

4.3 Removing non-dynamic modes

At the end of both Procedure GRCF and Procedure GRCFID, the matrices of the resulting realizations of the factors can be partitioned in the form

[TABLE]

where $A_{\infty}$ is a $(n-r)\times(n-r)$ upper triangular matrix, with $r=\mathop{\mathrm{rank}}E$, and the pair $(A_{\infty},0)$ contains the $n-r$ simple infinite eigenvalues (i.e., the non-dynamic modes), while the pair $(\widetilde{A}_{22},\widetilde{E}_{22})$ is in a GRSF with all eigenvalues in $\mathds{C}_{g}$. Using state residualization formulas, the $n-r$ simple infinite eigenvalues can be eliminated to obtain a reduced order descriptor system realization of the form

[TABLE]

where

[TABLE]

4.4 Computation of left coprime factorizations

If we apply Procedure GRCF or Procedure GRCFID to the dual realization $(A^{T}-\lambda E^{T},C^{T},B^{T},D^{T})$ of $G^{T}(\lambda)$, we obtain the RCF $G^{T}(\lambda)=N(\lambda)M^{-1}(\lambda)$, from which a left coprime factorization of $G(\lambda)$ results as G(\lambda)=\big{(}M^{T}(\lambda)\big{)}^{-1}N^{T}(\lambda). The realizations of $N(\lambda)$ and $M(\lambda)$ have the form (33) with the pair $(\widetilde{A},\widetilde{E})$ in an upper GRSF. This form is advantageous in many applications. For example, the eigenvalues of $\widetilde{A}-\lambda\widetilde{E}$ can be determined practically at no computational cost, and can be further reordered by preserving the upper quasi-triangular shape (GRSF) of the pair $(\widetilde{A},\widetilde{E})$. The dual realizations of $N^{T}(\lambda)$ and $M^{T}(\lambda)$ have the pair \big{(}\widetilde{A}^{T},\widetilde{E}^{T}\big{)} in a lower GRSF, which has some inconveniences for numerical computations. For example, the computation of eigenvalues using standard tools (e.g., MATLAB333MATLAB® is a registered trademark of The Mathworks, Inc.) involves performing the whole QZ algorithm [11] to obtain an upper GRSF. A simple trick can be used to obtain the dual realizations with an upper GRSF. This comes down to form the realizations of $N^{T}(\lambda)$ and $M^{T}(\lambda)$ as

[TABLE]

where $P$ is the permutation matrix of appropriate size with ones down on the secondary diagonal. Note that, in the above realization, the simple infinite eigenvalues are now located in the trailing position of the pair \big{(}P\widetilde{A}^{T}P,P\widetilde{E}^{T}P\big{)}.

5 Numerical examples

The proposed factorization procedures to compute RCFs of rational matrices have been implemented as MATLAB functions and belong to the free software collection of Descriptor Systems Tools (DSTOOLS) [19]. To illustrates the capabilities of the proposed computational algorithms, we present two simple examples computed with the functions grcf and grcfid available in DSTOOLS, which implements the Procedure GRCF and Procedure GRCFID, respectively. Both functions call the function gsorsf to compute the specially ordered GRSF. This function implements the Procedure GSORSF, but employs at Step 3, instead the method of [12], the Algorithm 3.2.1 of [20], for which an implementation is available in the SLICOT library [21].

Example 1

Consider the continuous-time improper TFM

[TABLE]

with the corresponding minimal realization given by

[TABLE]

$G(s)$ has the following set of poles: $\{-1,0,\infty,\infty\}$. To compute a stable and proper RCF of $G(s)$, we employed the pole assignment oriented factorization, with a stability degree of $\alpha=-1$ and the desired set of poles $\Gamma=\{-1,-2,-3\}$. With the option to eliminate the non-dynamic modes, the function grcf computes the following descriptor system realization of the factors (converted automatically to symbolic expressions)

[TABLE]

which correspond to the rational matrices

[TABLE]

The McMillan degree of $M(s)$ is three, thus the least possible one.

Example 2

Consider the discrete-time improper TFM

[TABLE]

with the corresponding minimal realization given by

[TABLE]

$G(z)$ has the following set of poles: $\{2,0,\infty,\infty\}$ and therefore, a RCF with inner denominator exists. With the option to eliminate the non-dynamic modes, the function grcfid computes the following descriptor system realization of the factors (converted automatically to symbolic expressions)

[TABLE]

which correspond to the rational matrices

[TABLE]

The McMillan degree of $M(z)$ is three, thus the least possible one. Interestingly, the McMillan degree of $N(z)$ is only two, because two unobservable eigenvalues in 0 have been removed. These eigenvalues are the zeros of the (improper) all-pass factor $\mathop{\mathrm{diag}}(z^{2},1)$ with two infinite poles, which is contained in $G(z)$.

6 Conclusions

In this paper we proposed two numerically reliable algorithms to compute right coprime factorizations of rational matrices using descriptor system based realizations. For the computation of a proper $\mathds{C}_{g}$-stable RCF, no particular properties of the underlying realizations need to be assumed. In particular, no $\mathds{C}_{b}$-stabilizability is required, because uncontrollable eigenvalues in $\mathds{C}_{b}$ are automatically removed. To determine factorizations with minimum degree denominator factors, the weakest requirement is that all controllable eigenvalues lying in $\mathds{C}_{b}$ are also observable. This requirement is fulfilled, for example, if the underlying descriptor system realization is minimal, or only observable, or $\mathds{C}_{b}$-detectable.

For the computation of a stable RCF with inner denominator, the condition for the lack of poles in $\partial\mathds{C}_{s}$ (the boundary of the stability domain) imposes the same condition on the eigenvalues of the pole pencil, provided the descriptor realization is minimal. However, the minimality condition on the descriptor system realization can be relaxed, by only requiring the lack of controllable eigenvalues in $\partial\mathds{C}_{s}$ (because the uncontrollable eigenvalues in $\partial\mathds{C}_{s}$ or $\mathds{C}_{b}$ are automatically removed). The minimum degree requirement for $M(\lambda)$ imposes additionally the lack of unobservable eigenvalues in $\mathds{C}_{s}$.

In the light of the above considerations, both the Procedure GRCF, in conjunction with the pole selection and assignment schemes presented in Section 4.1, as well as the Procedure GRCFID, in conjunction with the square-root-based computation of elementary inner factors described in Section 4.2, can be considered completely satisfactory numerical algorithms, which fulfill the standard requirements formulated in [22] for generality, numerical reliability and computational efficiency.