Partial Realization Theory and System Identification Redux

TL;DR

This paper revisits a nonstandard matrix Riccati equation from 20 years ago, deriving it anew within system identification, and demonstrating its universality for solving various analytic interpolation problems in control systems.

Contribution

It provides a new derivation of a nonstandard Riccati equation and shows its broad applicability to general analytic interpolation problems in control theory.

Findings

The Riccati equation is universal for interpolation problems.

The derivation connects realization theory with system identification.

It addresses a question posed by R.E. Kalman in the 1970s.

Abstract

Some twenty years ago we introduced a nonstandard matrix Riccati equation to solve the partial stochastic realization problem. In this paper we provide a new derivation of this equation in the context of system identification. This allows us to show that the nonstandard matrix Riccati equation is universal in the sense that it can be used to solve more general analytic interpolation problems by only changing certain parameters. Such interpolation problems are ubiquitous in systems and control. In this context we also discuss a question posed by R.E. Kalman in beginning of the 1970s.

Equations

$$\{\begin{array}{l}x(t+1)=Ax(t)+Bu(t)\\ y(t)=Cx(t)+Du(t)\end{array},$$

$$y_0,y_1,y_2,\dots ,y_N$$

$$\text{white noise}\stackrel{{\scriptstyle u}}{{\negmedspace}}{\longrightarrow}\framebox{$w(z)$}\negmedspace\stackrel{{\scriptstyle y}}{{\longrightarrow}}$$

$$w(z)=\rho \frac{\sigma (z)}{a(z)},$$

$$\sigma (z)$$

$$a(z)$$

$$f(z)=\frac{1}{2}\frac{b(z)}{a(z)}$$

$$\text{Re}\{f(e^{i\theta })\}=|w(e^{i\theta }){|}^2$$

$$b(z)=z^n+b_1{z}^{n-1}+\cdots +b_n$$

$$a(z)b(z^{-1})+b(z)a(z^{-1})=2{\rho }^2\sigma (z)\sigma (z^{-1}).$$

$$f(z)+f(z^{-1})={\rho }^{2}w(z)w(z^{-1})=:\Phi (z),$$

$$\Phi (z)=\sum _{k=-\infty }^{\infty }{c}_k{z}^{-k}$$

$$c_k=\mathbb{E}\{y(t+k)y(t)\}$$

$$f(z)=\frac{1}{2}+c_1{z}^{-1}+c_2{z}^{-2}+\dots$$

$$c_k=\underset{N\to \infty }{\lim }\frac{1}{N+1}\sum _{t=0}^{N-k}{y}_{t+k}{y}_t;$$

$$\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_d\end{array}\right]=2\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_d\end{array}\right]+\left[\begin{array}{cccc}1& & & \\ 2c_1& 1& & \\ ⋮& ⋮& & \\ 2c_{d-1}& 2c_{d-2}& \dots & 1\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_d\end{array}\right]$$

$$\left[\begin{array}{cccc}{c}_1& c_2& \cdots & c_d\\ c_2& c_3& \cdots & c_{d+1}\\ ⋮& ⋮& \ddots & ⋮\\ c_d& c_{d+1}& \cdots & c_{2d-1}\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_d\end{array}\right]=-\left[\begin{array}{c}{c}_{d+1}\\ c_{n+2}\\ ⋮\\ c_{2d}\end{array}\right]$$

$$d:=\mathrm{deg}f(z)=\text{rank}{H}_{\infty }=\text{rank}{H}_d.$$

$$c_k=\frac{1}{N+1}\sum _{t=0}^{N-k}{y}_{t+k}{y}_t$$

$$T=\left[\begin{array}{cccc}{c}_0& c_1& \cdots & c_n\\ c_1& c_0& \cdots & c_{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ c_n& c_{n-1}& \cdots & c_0\end{array}\right]$$

$$w(z)=\rho \frac{\sigma (z)}{a(z)}$$

$$\sigma =\left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ ⋮\\ {\sigma }_n\end{array}\right],\Gamma =\left[\begin{array}{ccccc}-{\sigma }_1& 1& 0& \cdots & 0\\ -{\sigma }_2& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -{\sigma }_{n-1}& 0& 0& \cdots & 1\\ -{\sigma }_n& 0& 0& \cdots & 0\end{array}\right],h=\left[\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right].$$

$$\begin{array}{rl}& \frac{z^n}{z^n+c_1{z}^{n-1}+\cdots +c_n}\\ & =1-u_1{z}^{-1}-u_2{z}^{-2}-u_3{z}^{-3}-\dots \end{array}$$

$$u=\left[\begin{array}{c}{u}_1\\ u_2\\ ⋮\\ u_n\end{array}\right],U=\left[\begin{array}{ccccc}0& & & & \\ u_1& 0& & & \\ u_2& u_1& & & \\ ⋮& ⋮& \ddots & & \\ u_{n-1}& u_{n-2}& \cdots & u_1& 0\end{array}\right].$$

$$g(P)=u+U\sigma +U\Gamma Ph.$$

$$P=\Gamma (P-Ph{h}^{\prime }P){\Gamma }^{\prime }+g(P)g(P{)}^{\prime },$$

$$\begin{array}{rl}& a=(I-U)(\Gamma Ph+\sigma )-u,\\ & \rho =\sqrt{1-h^{\prime }Ph}.\end{array}$$

$$f(z)=\frac{1}{2}+h^{\prime }(zI-F{)}^{-1}g,$$

$$F=\left[\begin{array}{ccccc}-a_1& 1& 0& \cdots & 0\\ -a_2& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -a_{n-1}& 0& 0& \cdots & 1\\ -a_n& 0& 0& \cdots & 0\end{array}\right]=J-a{h}^{\prime },$$

$$g=\frac{1}{2}(b-a),$$

$$\frac{b(z)}{a(z)}=1+2h^{\prime }(zI-F{)}^{-1}g,$$

Full text

I INTRODUCTION

A classical basic problem in system identification is to estimate the (unknown) constant coefficient matrices $A,B,C,D$ in a stable linear (SISO) stochastic system

[TABLE]

driven by white noise $\{u(t)\}_{t\in\mathbb{Z}}$, from a record of observations

[TABLE]

of the output process $\{y(t)\}_{t\in\mathbb{Z}}$, which is stationary in steady state. Modulo possible unobservable and/or unreachable modes, choice of coordinates and placement of zeros, this is equivalent to finding a shaping filter

[TABLE]

with a minimum-phase transfer function

[TABLE]

where

[TABLE]

are Schur polynomials, i.e., polynomials with all its roots in the open unit disc, and $\rho$ is a positive number. Then the rational function

[TABLE]

satisfying

[TABLE]

is positive real (see, e.g., [1]), and a simple calculation shows that

[TABLE]

is a Schur polynomial whose coefficients can be determined from the linear system of equations corresponding to the relation

[TABLE]

In fact, (7) is equivalent to

[TABLE]

where

[TABLE]

is the power spectral density of the stationary output process $y$ and

[TABLE]

are the covariance lags [1]. Then,

[TABLE]

is analytic in the complement of the unit disc in the complex plane and maps to the right half plane, estabishing the positive-real property.

Now, if we had an infinite observation record (2), i.e, $N=\infty$, then we would have an infinite sequence $(c_{0},c_{1},c_{2},\dots)$ of covariance lags that could be determined from the ergodic limits

[TABLE]

see, e.g., [1]. Then, identifying coefficients of powers of $z$ in $b(z)=2f(z)a(z)$ as done in [2], we obtain

[TABLE]

for nonnegative powers and

[TABLE]

for negative powers. The coefficient matrix in (14) is a Hankel matrix that we denote $H_{d}$. By Kronecker’s theorem [3],

[TABLE]

Hence, $f(z)$, and hence also $w(z)$, can be determined from finite sequence $(c_{0},c_{1},\dots,c_{n})$ of covariance lags, where $n:=2d$.

Therefore it may seem that only a finite covariance record $(c_{0},c_{1},\dots,c_{n})$ is needed, and this also a basic assumption in the early work on so-called subspace identification [4, 5, 6], where in general the biased ergodic estimates

[TABLE]

were used to insure that the corresponding Toeplitz matrix is positive definite, as required. However, as we pointed out in [7], this is incorrect and may lead to an $f(z)$ that is not positive real; also see [1, Chapter 13]. This is due to the possible difference between algebraic degree and positive degree (Section II-B). The classical subspace identification procedures were based on solving the deterministic partial realization problem [3, 9] rather than the stochastic one, namely the rational covariance extension problem, which we shall return to in Section II.

The focus of this paper will be on a certain nonstandard Riccati equation, called the Covariance Extension Equation (CEE), which solves the rational covariance extension problem. It will be introduced in Section II-C. Section II as a whole is devoted to background material for the paper. In addition, in Section II-D we consider a question posed by Kalman and give a partial answer. In Section III we provide a new derivation of CEE that will allow us to establish that versions of CEE can be used to solve more general analytic interpolation problems encountered in systems and control. Such an example will be given in Section IV. We conclude with Section V, where we some future directions of research are discussed.

II BACKGROUND

II-A Rational covariance Extension problem

The rational covariance extension problem, first formulated by Kalman in [8], can be stated in the following way. Given a positive covariance sequence $(c_{0},c_{1},\dots,c_{n})$, i.e., a sequence $(c_{0},c_{1},\dots,c_{n})$ with the property that the Toeplitz matrix

[TABLE]

is positive definite. find an infinite extension $c_{n+1},c_{n+2},\dots$ such that the function (12) is a rational positive-real function of degree at most $n$. This is the proper partial realization problem connected to the system identification problem posed in the introduction. If we ignore the condition that $f(z)$ be positive real, we have a deterministic partial realization problem equivalent to Padé approximation [9]. This is the problem solved in the original papers on subspace identification.

The following theorem provides a smooth complete parameterization of the set set of solutions to rational covariance extension problem.

Theorem 1

Let $(c_{0},c_{1},\dots,c_{n})$ be a positive covariance sequence. Then, given any Schur polynomial (4), there is one and only one Schur polynomial (5) and $\rho>0$ such that

[TABLE]

is a shaping filter for $(c_{0},c_{1},\dots,c_{n})$. The mapping from $\sigma$ to $(a,\rho)$ is diffeomorphism.

The existence part of Theorem 1 was proved in [10] (also see [11]) and the rest of the theorem in [12]. Note that $\sigma(z)$ and $a(z)$ may have common roots, so the degree of $w(z)$ might be less than $n$.

For each parameter $\sigma$ there is a convex optimization problem solving for $(a,\rho)$, which first appeared in [13] (also see [14, 15]), but here we shall consider a different method of solution described in Section II-C.

II-B Algebraic and positive degree

The algebraic degree of $(c_{0},c_{1},\dots,c_{n})$ is the minimal degree of a deterministic partial realization of $(c_{0},c_{1},\dots,c_{n})$. It is related to the rank of a Hankel matrix [9] and has the generic value $[\tfrac{n}{2}]$, i.e., sequences $(c_{0},c_{1},\dots,c_{n})$ that fail to have this algebraic degree live on a thin (lower-dimensional) subset, and there is no nonempty open set of vectors $c=(c_{0},c_{1},\dots,c_{n})^{\prime}\in\mathbb{R}^{n+1}$ having an algebraic degree different from $[\tfrac{n}{2}]$.

The positive degree of $(c_{0},c_{1},\dots,c_{n})$ is the minimum degree of any solution to the corresponding rational covariance extension problem. By contrast, it does not have a generic value, as seen from the following theorem proved in [2].

Theorem 2

For each $p$ such that $[n/2]\leq p\leq n$, there is a nonempty open set in $\mathbb{R}^{n+1}$ of sequences $(c_{0},c_{1},\dots,c_{n})$ for which $p$ is the positive degree. The maximal positive degree is $n$.

Since the original subspace identification algoritms are based on Hankel factorization, they produce solutions with the algebraic degree, which may not coincide with the required positive degree. In [16] we presented a system that produces data for which there is massive failure of the basic subspace identification algorithms. This led to quite a lot of work devising ad hoc fixes that are now included in the codes.

II-C Covariance extension equation

Given the Schur polynomial (4), we define

[TABLE]

Moreover, we represent the covariance data in terms of the first $n$ coefficients in the expansion

[TABLE]

about infinity and define

[TABLE]

Finally define the the function $g:\mathbb{R}^{n\times n}\to\mathbb{R}^{n}$given by

[TABLE]

Then the following nonstandard Riccati equation

[TABLE]

where ′ denotes transposition, was called the Covariance Extension Equation (CEE) in [2]. It parameterizes the solutions to the rational covariance extension problem in terms of the covariance data and the numerator polynomial $\sigma(z)$ corresponding to desired spectral zeros. In fact, the following theorem was proved in [2].

Theorem 3

Let $(c_{0},c_{1},\dots,c_{n})$ be a positive covariance sequence. Then, for each Schur polynomial (4), there is a unique symmetric solution $P$ of CEE satisfying $h^{\prime}Ph<1$. Moreover, for each $\sigma$ there is a unique shaping filter (3) for $(c_{0},c_{1},\dots,c_{n})$, where $a(z)$ and $\rho$ are given in terms of the corresponding $P$ by

[TABLE]

Here $a:=(a_{1},a_{2},\dots,a_{n})^{\prime}$. The degree of $w(z)$ equals the rank of $P$.

In [17, 18] we presented a homotopy continuation algorithm for solving CEE.

II-D Kalman’s question

In view of the fact that the algebraic degree can be determined from the rank of a Hankel matrix of the covariance data [3, 9], in 1972 Kalman [19] posed the question whether there is a similar matrix-rank criterion for determining the positive degree. In view of Theorem 2, this would seem impossible. The closest we have found in this direction is the following result, which follows from Theorem 3 and is reported in [2].

Proposition 4

Let $P(\sigma)$ be the unique solution of (20). Then the positive degree of the covariance sequence $(c_{0},c_{1},\dots,c_{n})$ equals the minimum of $\text{rank}\,P(\sigma)$ over all Schur polynomials $\sigma(z)$.

III A NEW DERIVATION OF THE COVARIANCE EXTENSION EQUATION

We present a new derivation of CEE which is divided into two separate steps, namely one that just imposes the condition that $f(z)$ be positive real and rational of degree at most $n$, and a second one imposing the interpolation condition that $f(z)$ should match the $n+1$ first covariances $(c_{0},c_{1},\dots,c_{n})$. In this way, we see that this nonstandard matrix Riccati equation is universal in the sense that it can be used to solve more general analytic interpolation problems only changing certain parameters. In fact, the first step remains the same in this more general context.

III-A Stochastic realization

The rational function $f(z)$ defined by (6) has a realization

[TABLE]

where

[TABLE]

$J$ is the upward shift matrix, and $g$ is an $n$-vector to be determined. Note that this need not be a minimal realization, as there could be cancellations of common zeros of $a(z)$ and $b(z)$.

Lemma 5

The vector $g$ in (22) is given by

[TABLE]

where $a:=(a_{1},a_{2},\dots,a_{n})^{\prime}$ and $b:=(b_{1},b_{2},\dots,b_{n})^{\prime}$ are the $n$-vectors of coefficients in the polynomials $a(z)$ and $b(z)$, respectively, in (6).

Proof:

From (6) and (22) we have

[TABLE]

to which we apply the matrix inversion lemma (Appendix) to obtain

[TABLE]

Hence, since $b(z)$ is the denominator polynomial, we must have $a+2g=b$, from which (24) follows. ∎

In the same way, $w(z)$, given by (3), has a realization

[TABLE]

for some $n$-vector $k$. Since $w(z)$ is a minimum-phase spectral factor of $f(z)+f(z^{-1})$, from stochastic realization theory [1, Chapter 6] we have that

[TABLE]

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

[TABLE]

Lemma 6

The vectors $g$ and $k$ in (22) and (25) are given by

[TABLE]

where $P$ is the minimal solution of (27).

Proof:

Applying the matrix inversion lemma (Appendix) to

[TABLE]

yields

[TABLE]

However, since the denominator is $\sigma(z)$, we must have

[TABLE]

and hence (28b) follows. Moreover, $\Gamma=F-\rho^{-1}k$. This together with (28b) yields

[TABLE]

from which (28a) follows. ∎

Lemma 7

The minimal solution of (27) is also the minimal solution of

[TABLE]

Proof:

Note that (27) can be written

[TABLE]

from which we have

[TABLE]

which, in view of (28), in turn yields (29). ∎

III-B Interpolation condition

Next we introduce the interpolation condition that $f(z)$ matches the first $n+1$ covariance lags $c_{0},c_{1},\dots,c_{n}$. To this end, we identify coefficients of nonnegative powers in $b(z)=2f(z)a(z)$ with $f(z)$ given by (12). This yields

[TABLE]

where

[TABLE]

However, by Lemma 5, $b=a+2g$, which inserted into (30) yields $g+a=c+Ca$. Combining this with (28a) we have $a=C^{-1}(\Gamma Ph+\sigma-c)$, and hence (28a) can be written

[TABLE]

Note that here $C$ is an $n\times n$ matrix and not a $d\times d$ matrix with $d$ being the algebraic degree of $c_{0},c_{1},\dots,c_{n}$ as in (13). By Theorem 2, $n$ is an upper bound of the positive degree.

Finally we show that (32) is equivalent to (19). To this end, we first identify negative powers of $z$ in

[TABLE]

obtained from (17), to obtain

[TABLE]

from which we have

[TABLE]

Consequently, $C^{-1}c=u$ and $I-C^{-1}=U$, so (19) follows from (19).

IV GENERALIZATION TO ANALYTIC INTERPOLATION

Next we show that the nonstandard Riccati equation (20) is universal in the sense that it holds for more general analytic interpolation problems by merely redefining the parameters $u$ and $U$. We shall demonstrate this for a Nevanlinna-Pick interpolation problem with rationality constraints.

Given distinct point $z_{0},z_{1},\dots,z_{n}$ in the complement of the unit disc of the complex plane and points $c_{0},c_{1},\dots,c_{n}$ in the open right half-plane, find a rational positive real function $f(z)$ of degree at most $n$ satisfying the interpolation condition

[TABLE]

It is convenient to chose the points in conjugate pairs to ensure that $f$ is a real function.

Clearly all calculations in Section III-A remain intact, and it on remains to enforce the interpolation condition, which we may write as

[TABLE]

or, equivalently, as

[TABLE]

Since the points $z_{0},z_{1},\dots,z_{n}$ are distinct, the Vandermonde matrix (35b) is nonsingular, and hence we have

[TABLE]

Therefore, by Lemma 5,

[TABLE]

Consequently, in view of (28a),

[TABLE]

Assuming that $I+T$ is nonsingular, for the moment as a technical condition to be more carefully investigated in the future, we define

[TABLE]

where for clarity we have added an index to each identity matrix to indicate dimension. Then

[TABLE]

which has the same form as (19). In combination with (28a) and (26), this also yields

[TABLE]

which is of the same form as (21).

We have thus demonstrated that the Covariance extension equation can be used also in this case after changing the definition of the interpolation parameters $(u,U)$. Hence the algorithms using homotopy continuation presented in [17, 18] could also be used here.

V CONCLUSIONS

In this paper we provide a new derivation of a nonstandard Riccati equation (CEE) for rational covariance extension that separates the part that only depends on the positivity and rationality constraints and the part that depends on the interpolation condition. In this way we see that the structure of CEE remains intact, and only certain interpolation parameters need to be modified when treating more general analytic interpolation problems with rationality constraints.

It should be possible to generalize this framework to the MIMO case. This, together with some remaining numerical and technical issues, will the topic of a future paper.

APPENDIX

For ease of reference, we here reproduce the well-known matrix inversion lemma. Provided all inverses exist, the formula

[TABLE]

holds for otherwise arbitrary matrices of compatible dimensions. This is seen by direct computation.