Multivariable analytic interpolation with complexity constraints: A modified Riccati approach

TL;DR

This paper introduces a new multivariable analytic interpolation method with complexity constraints, extending previous scalar results, and employs a homotopy continuation algorithm to solve Riccati-type matrix equations.

Contribution

It generalizes scalar interpolation techniques to multivariable cases, addressing new challenges and providing an algorithm for solving associated Riccati equations.

Findings

Developed a new multivariable interpolation method.

Extended scalar results to multivariable systems.

Provided an algorithm based on homotopy continuation for Riccati equations.

Abstract

Analytic interpolation problems with rationality and derivative constraints occur in many applications in systems and control. In this paper we present a new method for the multivariable case, which generalizes our previous results on the scalar case. This turns out to be quite nontrivial, as it poses many new problems. A basic step in the procedure is to solve a Riccati type matrix equation. To this end, an algorithm based on homotopy continuation is provided.

Equations

$$F(e^{i\theta })+F(e^{-i\theta }{)}^{\prime }>0,-\pi \le \theta \le \pi ,$$

$$\frac{1}{k!}{F}^{(k)}(z_j)=W_{jk},$$

$$k=0,\cdots n_j-1,$$

$$n=(\sum _{j=0}^m{n}_j-1).$$

$$F(z)=\frac{1}{2}I+zH(I-zF{)}^{-1}G,$$

$$W:=\left[\begin{array}{ccc}{W}_0& & \\ & \ddots & \\ & & W_m\end{array}\right]$$

$$W_j=\left[\begin{array}{cccc}{W}_{j0}& & & \\ W_{j1}& W_{j0}& & \\ ⋮& \ddots & \ddots & \\ W_{j{n}_{j-1}}& \cdots & W_{j1}& W_{j0}\end{array}\right]$$

$$Z:=\left[\begin{array}{ccc}{Z}_0& & \\ & \ddots & \\ & & Z_m\end{array}\right]$$

$$Z_j=\left[\begin{array}{cccc}{z}_j& & & \\ 1& z_j& & \\ & \ddots & \ddots & \\ & & 1& z_j\end{array}\right]j=0,1,\dots ,m.$$

$$e:=[e_{n_0}^1,e_{n_1}^1,\cdots ,e_{n_m}^1{]}^{\prime },$$

$$S=ZS{Z}^{\ast }+e{e}^{\prime }.$$

$$W(S\otimes I_{\ell })+(S\otimes I_{\ell })W^{\ast }>0$$

$${\Phi }_{+}(z)=\frac{1}{2}I+H(zI-F{)}^{-1}G,$$

$${\Phi }_{+}(e^{i\theta })+{\Phi }_{+}(e^{-i\theta }{)}^{\prime }>0,-\pi \le \theta \le \pi ,$$

$$H=\text{diag}(h_{t_1},h_{t_2},\dots ,h_{t_{\ell }})\in {\mathbb{R}}^{\ell \times n\ell }$$

$$F=J-AH\in {\mathbb{R}}^{n\ell \times n\ell }$$

$$J_{\nu }=\left[\begin{array}{ccccc}0& 1& 0& \dots & 0\\ 0& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & 0\\ 0& 0& 0& \dots & 1\\ 0& 0& 0& \dots & 0\end{array}\right]$$

$$t_1+t_2+\cdots +t_{\ell }=n\ell .$$

$$A(z)=D(z)+\Pi (z)A,$$

$$D(z):=\text{diag}(z^{t_1},z^{t_2},\dots ,z^{t_{\ell }}).$$

$$\Pi (z)(zI-J)=\text{diag}(z^{t_1},z^{t_2},\dots ,z^{t_{\ell }})H.$$

$$\Pi (z)(zI-F)=\Pi (z)(zI-J)+\Pi (z)AH=A(z)H$$

$${\Phi }_{+}(z)=\frac{1}{2}A(z{)}^{-1}B(z),$$

$$B(z)=D(z)+\Pi (z)B$$

$$B=A+2G.$$

$$V(z)V(z^{-1}{)}^{\prime }=\Phi (z):={\Phi }_{+}(z)+{\Phi }_{+}(z^{-1}{)}^{\prime }.$$

$$V(z)=H(zI-F{)}^{-1}K+R,$$

$$V(z)=A(z{)}^{-1}\Sigma (z)R,$$

$$\Sigma (z)=D(z)+\Pi (z)\Sigma$$

$$\Sigma =A+K{R}^{-1}.$$

$$K$$

Full text

I Introduction

A common problem in robust control and spectral estimation is to find an $\ell\times\ell$ matrix-valued rational function $F$, analytic in the unit disc $\mathbb{D}=\{z\mid|z|>1\}$, such that

[TABLE]

which also satifies the interpolation condition

[TABLE]

where ′ denotes transposition, $F^{(k)}(z)$ is the $k$th derivative of $F(z)$, and $z_{0},z_{1},\dots,z_{m}$ are distinct points in $\mathbb{D}$. We restrict the complexity of the rational function $F(z)$ by requiring that its McMillan degree be at most $\ell n$, where

[TABLE]

Without loss of generality we may assume that $z_{0}=0$ and $W_{0}=\frac{1}{2}I$. Then $F(z)$ has a realization

[TABLE]

where $H\in\mathbb{R}^{\ell\times\ell n}$, $F\in\mathbb{R}^{\ell n\times\ell n}$, $G\in\mathbb{R}^{\ell n\times\ell}$, the matrix $F$ has all its eighenvalues in $\mathbb{D}$ and $(H,F)$ is an observable pair.

Let $W$ be the $\ell(n+1)\times\ell(n+1)$ matrix

[TABLE]

with

[TABLE]

for each $j=0,1,\dots,m$. Moreover, let $Z$ be the $(n+1)\times(n+1)$ matrix

[TABLE]

with

[TABLE]

Finally define the $n+1$-dimensional column vector

[TABLE]

where $e_{n_{j}}^{1}=[1,0,\cdots,0]\in\mathbb{R}^{n_{j}}$ for each $j=0,1,\dots,m$, and let $S$ be the unique solution of the Lyapunov equation

[TABLE]

Note that the eigenvalues of $Z$ are all located in the open unit disc $\mathbb{D}$.

The problem of determining the interpolant $F(z)$ is an inverse problems which has a solution if and only if

[TABLE]

(see, e.g., [22]), and then there are an infinite number of solutions. We would like to find a parametrization of these solutions.

The special case when $\ell=1$, $m=0$ and $n_{0}=n+1$ is called the rational covariance extension problem and was first formulated by Kalman [1] and then solved in steps in [2, 3, 4, 5], where a complete parameterization in terms of spectral zeros was obtained, and in [6, 7], where a convex optimization approach was introduced. This problem have occurred in many applications in systems and control such as in signal and speech processing [8] and in identification [9]. The case $n_{0}=n_{1}=\dots=n_{m}=1$ and $m=n$ is called the Nevanlinna-Pick interpolation problem with degree constraint and was early considered in robust control [10] and many other applications in systems and control [11, 12]. It was completely parameterized, again in steps, in [13, 14, 15, 16], and a convex optimization approach was introduced in [15, 16]. Since then a large number of papers on the more general scalar problem has appeared [17, 18, 19, 20, 21]. We refer to [9] for further references.

The multivariable case $(\ell>1)$ is much harder, and the nice spectral-zero assignability present in the scalar case appears to be lost or at lease elusive. Restrictive classes of such problems have been considered in large number of papers [2, 22, 23, 24, 25, 26, 28, 27, 29, 30], but the theory remains incomplete, and many problems have been left open.

In [31] we presented a complete parameterization of the problem presented above for the scalar case ($\ell=1$) in terms of a modified Riccati equation, which was first introduced for more restricted classes of interpolation problems in [5] and [32]. As [5] studied the rational covariance extension problem, the modified Riccati equation was named the Covariance Extension Equation (CCE), and we retain this name although the problems now considered are much more general.

In the present paper we take a first step in generalizing the results in [31] to the multivariable case $(\ell>1)$. In Section II we provide the basic tools for the multivariable problem. To describe our ultimate goal we provide in Section III a brief review of the scalar results in [31], and then in Section IV we develop the multivariable case in the same spirit. In Section V we present our main results and an algorithm based on homotopy continuation in the style of [33, 31]. The results fall somewhat short of what the scalar case promises, and, given some results in [27], we suspect that this is due to problems introduced by the nontrivial Jordan structure of the multivariable case. In SectionVI-C we provide some simulations to illustrate this and also an example of model reduction. Finally, in Section VII we give some concluding remarks and suggestions for future research.

II Preliminaries

Defining $\Phi_{+}(z):=F(z^{-1})$ we have

[TABLE]

which has all its poles in the unit disc $\mathbb{D}$. In view of (1)

[TABLE]

and hence $\Phi_{+}(z)$ is (strictly) positive real [9, Chapter 6]. By a coordinate transformation $(H,F,G)\to(HT^{-1},TFT^{-1},TG)$ we can choose $(H,F)$ in the observer canonical form

[TABLE]

with $h_{\nu}:=(1,0,\dots,0)\in\mathbb{R}^{\nu}$, and

[TABLE]

where $J:=\text{diag}(J_{t_{1}},J_{t_{2}},\dots,J_{t_{\ell}})$ with $J_{\nu}$ the $\nu\times\nu$ shift matrix

[TABLE]

and $A\in\mathbb{R}^{n\ell\times\ell}$. The numbers $t_{1},t_{2},\dots,t_{\ell}$ are the observability indices of $\Phi_{+}(z)$, and

[TABLE]

Next define $\Pi(z):=\text{diag}(\pi_{t_{1}}(z),\pi_{t_{2}}(z),\dots,\pi_{t_{\ell}}(z))$, where $\pi_{\nu}(z)=(z^{\nu-1},\dots,z,1)$, and the $\ell\times\ell$ matrix polynomial

[TABLE]

where

[TABLE]

Lemma 1

$H(zI-F)^{-1}=A(z)^{-1}\Pi(z)$**

Proof:

First note that

[TABLE]

Then

[TABLE]

as claimed. ∎

Consequently

[TABLE]

where

[TABLE]

with

[TABLE]

Moreover let $V(z)$ be the minumum-phase spectral factor of

[TABLE]

We know [9, Chapter 6] that $V(z)$ has a realization of the form

[TABLE]

which, by Lemma 1, can be written

[TABLE]

where

[TABLE]

with

[TABLE]

From stochastic realization theory [9, Chapter 6] we have

[TABLE]

where $P$ is the unique minimum solution of the algebraic Riccati equation

[TABLE]

Now, from (13), (22) and (23) we have

[TABLE]

where, in view of (21),

[TABLE]

Hence

[TABLE]

Since $F=\Gamma+KR^{-1}H$ and $G-\Gamma PH^{\prime}=KR^{-1}$, (24) can be written

[TABLE]

where we have also used (23). Inserting $KR^{-1}=G-\Gamma PH^{\prime}$ we have

[TABLE]

III A review of the scalar case

To motivate our approach to the multivariable problem presented in Section I we shall briefly review some results on the scalar case ($\ell=1$) presented in [31]. To stress the fact that the matrices $H,G,A,B$ and $\Sigma$ are now $n$-vectors we shall here denote them $h,g,a,b$ and $\sigma$ instead, and the scalar $R$ will be denoted $\rho$.

Introducing the interpolation conditions (2) into the calculation we obtain

[TABLE]

where $P$ is the unique solution of the algebraic Riccati equation (27). Elimination $g$ we then have the modified Riccati equation

[TABLE]

where $u\in\mathbb{R}^{n}$ and $U\in\mathbb{R}^{n\times n}$ are given by

[TABLE]

with

[TABLE]

We showed in [31] that there is a map sending $W$ to $u$ which is a diffeomorphism and that there is a linear map $L$ such that $U=Lu$.

Let $\mathcal{S}_{n}$ be the space of Schur polynomials (i.e., polynomials with all zeros in the open unit disc $\mathbb{D}$) of the form

[TABLE]

and let $\mathcal{P}_{n}$ be the $2n$-dimensional space of pairs $(a,b)\in\mathcal{S}_{n}\times\mathcal{S}_{n}$ such that $b(z)/a(z)$ is positive real. Moreover, for each $\sigma\in\mathcal{S}_{n}$, let $\mathcal{P}_{n}(\sigma)$ be the submanifold of $\mathcal{P}_{n}$ for which

[TABLE]

holds, where $\rho^{2}$ is the appropriate normalizing factor. It was shown in [34] that $\{\mathcal{P}_{n}(\sigma)\mid\sigma\in\mathcal{S}_{n}\}$ is a foliation of $\mathcal{P}_{n}$, i.e., a family of smooth nonintersecting submanifolds, called leaves, which together cover $\mathcal{P}_{n}$. Moreover, for any polynomial (31), let $a_{*}(z)=z^{n}a(z^{-1})$ be the reversed polynomial of $a(z)$. Finally, let $\mathcal{W}_{+}$ be the space of all $W$ such that the generalized Pick matrix $WS+SW^{*}$ is positive definite, where $S$ the unique solution of the Lyapunov equation (10) .

The following result was proved in [31].

Theorem 2

Let $\ell=1$. For each $(\sigma,W)\in\mathcal{S}_{n}\times\mathcal{W}_{+}$, the modified Riccati equation (28) has a unique positive definite solution $P$ such that $hPh^{\prime}<1$, and the problem to find a rational function $b_{*}(z)/a_{*}(z)$ satisfying the interpolation conditions (2) and the positivity condition (32) has a unique solution given by

[TABLE]

In fact, the map sending $(a,b)\in\mathcal{P}_{n}(\sigma)$ to $W\in\mathcal{W}_{+}$ is a diffeomorphism. Finally, the degree of the interpolant equals the rank of $P$.

Consequently, by (32), for each $\sigma\in\mathcal{S}_{n}$ there is a unique interpolant $b_{*}(z)/a_{*}(z)$ with the prescribed properties such that

[TABLE]

Hence

[TABLE]

is the correspondning spectral factor.

In [31] we solved (28) by homotopy continuation by taking $u(\lambda)=\lambda u$ with $\lambda$ varying from [math] to $1$. We showed that this provides an efficient and robust algorithm for analytic interpolation with degree constraint that can handle situations which are difficult with the optimization approach, especially when system poles are close to the unit circle.

IV The matrix case

Next we turn to the general multivariable case and introduce the interpolation condition (2) in the matrix setting of Section II.

Lemma 3

Let the matrices $W$ and $Z$ be given by (5) and (7), respectively. Then the interpolation condition (2) can be written

[TABLE]

where $\otimes$ denotes Kronecker product.

Proof:

Since $F(z)$ is analytic in $\mathbb{D}$, it has the representation

[TABLE]

for all $z\in\mathbb{D}$, where $C_{0}=\frac{1}{2}I_{\ell}$. A straight-forward but tedious calculation, omitted here for lack of space, yields

[TABLE]

where $W_{j}$, defined by (6), is given by (2). Then (34) follows from (7) and (5). ∎

Let $A_{*}(z)$ be the reversed matrix polynomial

[TABLE]

where $D(z)$ is given by (16), and let $B_{*}(z)$ be defined in the same way in terms of $B(z)$. Then

[TABLE]

and the interpolation condition (34) can be written

[TABLE]

Moreover, let the $\ell\times\ell n$ matrices $N_{1},N_{2},\dots,N_{t}$ be defined by

[TABLE]

where $t$ is the largest observability index. Then

[TABLE]

where $A_{k}=N_{k}A$. For later use we observe that

[TABLE]

where $e_{j}^{k}$ is a $1\times j$ row vector with the k:th element 1 and the others 0 whenever $k\leq j$, and a zero row vector of dimension $1\times j$ when $k>j$.

Next we reformulate (37) as

[TABLE]

where $M$ is the $\ell(n+1)\times\ell(n+1)(t+1)$ matrix

[TABLE]

In view of (18), (40) can be written

[TABLE]

or, equivalently,

[TABLE]

where $Q_{k}:=A_{k}+G_{k}$, $k=1,2,\dots,t$, and

[TABLE]

where

[TABLE]

for $j=0,1,\dots,m$.

Using the rule $(A\otimes B)(C\otimes D)=(AC)\oplus(BD)$, valid for arbitrary matrices of appropriate dimensions, (41) takes the form

[TABLE]

Then multiplying both sides from the right by $(e\otimes I_{\ell})$ and observing that

[TABLE]

we obtain

[TABLE]

where $V$ is the $\ell(n+1)\times\ell t$ matrix

[TABLE]

and $\hat{T}$ is the $\ell(n+1)\times\ell$ matrix

[TABLE]

that is

[TABLE]

Therefore, since $G_{k}=N_{k}G$ and $Q_{k}=N_{k}Q$, we have

[TABLE]

where $N$ is given by (39). Now, $VN$ is an $\ell(n+1)\times\ell n$ matrix in which the top $\ell$ rows are zero, since $z_{0}=0$. Therefore it takes the form

[TABLE]

For the moment assuming that the square matrix $L$ is nonsingular – we shall later see that this is not always true – $VN$ has a psuedo-inverse $(VN)^{\dagger}$, and hence (48) yields

[TABLE]

Since, by definition, $Q=A+G$, (26) yields

[TABLE]

where $u:=(VN)^{\dagger}\hat{T}$ and $U:\mathbb{R}^{\ell n\times\ell}\to\mathbb{R}^{\ell n\times\ell}$ is the linear operator

[TABLE]

Inserting (51) into (27) we obtain the modified Riccati equation

[TABLE]

It was first introduced in [5] for the scalar case $\ell=1$ and for the special case of covariance extension. Therefore it has been called the Covariance Extension Equation (CEE).

V Main results

Next we generalize the results of Section III to the general multivariable problem, which is considerably more difficult. Therefore several key questions will be left unanswered at this time. Nevertheless the theory in its present (preliminary) form does give an workable algorithm for large classes of problems.

V-A Basic results

Now redefine $\mathcal{S}_{n}$ to be the class of $\ell\times\ell$ matrix polynomials (15) such that $\det A(z)$ has all its zeros in the open unit disc $\mathbb{D}$. Clearly $\mathcal{S}_{n}$ consists of subclasses with different Jordan structure $J$ defined via (13). In each such subclass $D(z)$ and $\Pi(z)$ in (15), as well as $N_{1},N_{2},\dots,N_{t}$ in (38), are the same. Let $\mathcal{W}_{+}$ we the values in (2) that satisfy the geralized Pick condition (11).

Lemma 4

Let the $\ell n\times\ell n$ matrix $L$ be defined by (49). Then $L$ is nonsingular if and only if all observability indices are the same, i.e., $t_{1}=t_{2}=\dots=t_{\ell}=n$.

Proof:

Let us order the observability indices as $t_{1}\geq t_{2}\geq\dots\geq t_{\ell}$ and set $t:=t_{1}$. Then, by (14), $t\geq n$. Since $(Z,e)$ is a reachable pair,

[TABLE]

First assume that $t=n$. Then, since $\text{rank}(A\otimes B)=\text{rank}(A)\text{rank}(B)$,

[TABLE]

has rank $\ell n$, and $N\in\mathbb{R}^{n\ell\times n\ell}$ given by

[TABLE]

also has rank $\ell n$. Consequently Sylverster’s inequality,

[TABLE]

(see, e.g., [9, p.741]) implies that $VN$ has rank $\ell n$, and hence $L$ is nonsingular.

Next assume that $t>n$. Then the first $t$ columns of $N$, now given by (39) can be written $I_{t}\otimes(e^{1}_{\ell})^{\prime}$, so the first $t$ columns of $VN$ form the matrix

[TABLE]

which in view of (53) has rank $n<t$. Hence the columns of $VN$ are linearly dependent, and thus $L$ is singular. ∎

In the present matrix case, the relation (32) reads

[TABLE]

Let $\mathcal{P}_{n}$ be the space of pairs $(A,B)\in\mathcal{S}_{n}\times\mathcal{S}_{n}$ such that $A(z)^{-1}B(z)$ is positive real. Then the problem at hand is to find, for each $\Sigma\in\mathcal{S}_{n}$, a pair $(A,B)\in\mathcal{P}_{n}$ such that (55) and (2) hold.

Theorem 5

Given $(\Sigma,W)\in\mathcal{S}_{n}\times\mathcal{W}_{+}$, where $\Sigma(z)$ has all it observability indicies equal. Then to any positive definite solution $P$ of the Covariance Extension Equation (52) such that $HPH^{\prime}<I$, there corresponds a unique analytic interpolant (36), where $A(z)$ and $B(z)$ have the same Jordan structure as $\Sigma(z)$, the matrices $A$ and $B$ are given by

[TABLE]

and $A(z)$ and $B(z)$ satisfy (55) with

[TABLE]

Finally,

[TABLE]

Proof:

The theorem follows from the derivation above. For the details of the proof of (58) we refer to [5]. ∎

Let us stress that these results are considerably weaker than the corresponding theorem for the scalar case reviewed in Section III. In fact, Theorem 5 does not guarantee that there exists a unique solution to (52). In fact, if there were two solutions to (52), there would be two interpolants, a unique one for each solution $P$. Moreover, the condition on the observability indices restricts the classes of Jordan structures that are feasible.

However, Theorem 5 can be combined with other partial results on existence and uniqueness. There are multivariable problems for which we already know that there is a unique solution to the interpolation problem, and then existence and uniquenss of a solution to (52) will follow. A case in point is when $\Sigma(z)=\sigma(z)I$, where $\sigma(z)$ is a stable scalar polynomial [22, 25], in which case the the observability indices are all equal, as required in Theorem 5. In this case the analytic interpolation problem will have a unique solution, and thus, tracing the calculations in Section II backwards, so will (52). The same is true when $\Sigma(z)=\sigma(z)C$, where $C$ is full rank [25].

On the other hand, in recent years there have been a number of results [25, 26, 27, 29, 30] on the question of existence and uniqueness of the multivariate analytic interpolation problem, mostly for the covariance extension problem ($m=0,n_{0}=n+1$), but there are so far only partial results and for special structures of the prior (in our case $\Sigma(z)$). Especially the question of uniqueness has proven elusive. Perhaps, as suggested in [27], this is due to the Jordan structure, and this could be the reason for the condition on the observability indices required in Theorem 5. In any case, as long as our algorithm delivers a solution to the Covariance Extension Equation, we will have a solution to the analytic interpolation problem, unique or not. An advantage of our method is that (58) can be used for model reduction, as will be illustrated in Section VI-C.

V-B Solving CEE by homotopy continuation

We shall provide an algorithm for solving (52) based on homotopy continuation. We assume from now on that $t:=t_{1}=t_{2}=\dots,t_{\ell}=n$. Whenever this algoritm delivers a solution $P$, the interpolant is obtained via (56).

When $u=0$, $\hat{T}=0$, and hence $U=0$. Then the modified Riccati equation (52) becomes $P=\Gamma(P-PH^{\prime}HP)\Gamma^{\prime}$, which has the solution $P=0$. We would like to make a continuous deformation of $u$ to go from this trivial solution to the solution of (52), so we choose $u(\lambda)=\lambda u$ with $\lambda\in[0,1]$. The corresponding deformation of $U$ is $\lambda U$, and $T$ is deformed to $\lambda T$. The value matrix (5) will vary as

[TABLE]

and we need to ascertain that $W(\lambda)$ still satisfies (11) along the whole trajectory.

Lemma 6

Suppose that $W\in\mathcal{W}_{+}$. Then $W(\lambda)\in\mathcal{W}_{+}$ for all $\lambda\in[0,1]$.

Proof:

By (11) we want to show that

[TABLE]

for $E:=S\otimes I_{\ell}$. We have

[TABLE]

which we know is positive definite for $\lambda=1$. However,

[TABLE]

and and therefore (59) holds a fortiori. ∎

Now, note that equation (55) can be written as

[TABLE]

where

[TABLE]

Lemma 7

Let $N_{n}$ be defined by (54). Then

[TABLE]

where $A_{n}=N_{n}A$, $B_{n}=N_{n}B$ and $\Sigma_{n}=N_{n}\Sigma$.

Proof:

From (56) and (23) we have

[TABLE]

Since $e^{n}_{n}J_{n}=0$ and hence $N_{n}J=0$, the statement of the lemma follows. ∎

Applying Lemma 7 and deleting the zero row in (60), it can be reduced to $n\ell\times\ell$ equations

[TABLE]

Therefore, introducing the $n\ell\times\ell$ matrix

[TABLE]

we use the homotopy

[TABLE]

where

[TABLE]

depend on $(p,\lambda)$. Then the problem reduces to solving the differential equation

[TABLE]

[35], which has the solution $\hat{p}(\lambda)$ for $0\leq\lambda\leq 1$. The solution of (52) is then obtained by finding the unique solution of the Lyapunov equation

[TABLE]

VI Some simple illustrative examples

Next we provide a few simulations that illustrate the theory.

VI-A Example 1

We consider a problem with the interpolation constraints

[TABLE]

where $n=2,\ell=2$. This yields

[TABLE]

Taking $\Sigma(z)$ of the form

[TABLE]

we have $t_{1}=t_{2}=2$, and thus

[TABLE]

where clearly $L$ is nonsingular as required. Then choosing

[TABLE]

the matrix polynomial (67) is stable, and we can use our algorithm to obtain

[TABLE]

where

[TABLE]

If instead we choose $\Sigma(z)$ of the form

[TABLE]

then $t=t_{1}=3,t_{2}=1$, and thus

[TABLE]

where $L$ is singular as anticipated by Lemma 4. Hence the algorithm cannot be used.

VI-B Example 2

Next consider the multivariable covariance extension problem

[TABLE]

where

[TABLE]

In this case

[TABLE]

By Lemma 4, we need to choose $t_{1}=t_{2}=2$ and this yields

[TABLE]

where $L$ is nonsingular. Choosing the $\Sigma(z)$ in (68), our algoritm delivers the solution (69) with

[TABLE]

Again the choice $t=t_{1}=3,t_{2}=1$ of obsevability indices does not work (Lemma 4).

VI-C Model reduction

Consider a system with a $2\times 2$ transfer function

[TABLE]

of dimension $10$ and with observability indices $t_{1}=t_{2}=5$, where $A(z)$ and $\Sigma(z)$ are given by (15) repectively (20) with

[TABLE]

Here $\Sigma(z)=\sigma(z)I_{2}$, were

[TABLE]

We pass (normalized) white noise though the system

to obtain the output $y_{0},y_{2},\cdots,y_{N}$, and from this output data we estimate the $2\times 2$ matrix valued covariance sequence

[TABLE]

Then we solve the problem (2) with $\ell=2$, $m=0$, $n_{0}=6$, and $W_{0k}=\hat{C}_{k}$ for $k=0,2,\dots,5$. The modified Riccati equation (52) has a solution $P$ with eigenvalues

[TABLE]

The first four eigenvalues are very small, so we can reduce the degree of this system from ten to six by choosing the first four covariance lags $\hat{C}_{0},\hat{C}_{1},\hat{C}_{2},\hat{C}_{3}$ and removing zeros of $\Sigma(z)$ at $-0.2$ and $0.6$. The reduced-order system will have observability indices $t_{1}=t_{2}=3$. The singular values of the true system (70) are shown the multivariable Bode-type plot in Fig. 1 together with those of the estimated systems of degree ten and six, respectively. As can be seen, the reduced-order system is a good approximation.

VII Concluding remarks

We have extended our previous results [31] for the scalar case to the matrix case. However, multivariable versions of analytic interpolation with rationality constraints have been marred by difficulties to establish existence and, in particular, uniqueness in the various parameterizations [2, 22, 23, 24, 25, 26, 27, 29, 30], and we have encountered similar difficulties here. Our approach attacks these problems from a different angle and might put new light on these challenges. Therefore future research efforts will be directed towards settling these intriguing open questions in the context of the modified Riccati equation (52).