Optimal input design for system identification using spectral decomposition

TL;DR

This paper develops a spectral decomposition-based method for designing optimal band-limited inputs with power constraints to identify linear MIMO systems, providing bounds and solutions for a non-convex optimization problem.

Contribution

It introduces a novel spectral decomposition approach for optimal input design in system identification, including bounds and solutions for the non-convex optimization problem.

Findings

Bounds on the optimal input design can be efficiently computed.

The bounds often match, indicating global optimality.

Simulations validate the proposed method.

Abstract

The aim of this paper is to design a band-limited optimal input with power constraints for identifying a linear multi-input multi-output system. It is assumed that the nominal system parameters are specified. The key idea is to use the spectral decomposition theorem and write the power spectrum as $\phi_{u}(j\omega)=\frac{1}{2}H(j\omega)H^*(j\omega)$. The matrix $H(j\omega)$ is expressed in terms of a truncated basis for $\mathcal{L}^2\left(\left[-\omega_{\mbox{cut-off}},\omega_{\mbox{cut-off}}\right]\right)$. With this parameterization, the elements of the Fisher Information Matrix and the power constraints turn out to be homogeneous quadratics in the basis coefficients. The optimality criterion used are the well-known $\mathcal{D}-$optimality, $\mathcal{A}-$optimality, $\mathcal{T}-$optimality and $\mathcal{E}-$optimality. The resulting optimization problem is non-convex in general. A lower bound on the optimum is obtained through a bi-linear formulation of the problem, while an upper bound is obtained through a convex relaxation. These bounds can be computed efficiently as the associated problems are convex. The lower bound is used as a sub-optimal solution, the sub-optimality of which is determined by the difference in the bounds. Interestingly, the bounds match in many instances and thus, the global optimum is achieved. A discussion on the non-convexity of the optimization problem is also presented. Simulations are provided for corroboration.

Figures

Equations

$$y(t)=G(s;\theta )u(t)+e(t),$$

$$\begin{array}[]{l}\displaystyle\hskip 56.9055ptG_{ij}(s;\theta)=\frac{\sum_{k=0}^{m_{ij}}b_{ij}^{(k)}s^{k}}{\sum_{k=0}^{n_{ij}}a_{ij}^{(k)}s^{k}},\\ \displaystyle\theta^{ij}=\left\{b_{ij}^{(0)},~{}\dots,~{}b_{ij}^{(m_{ij})},~{}a_{ij}^{(0)},~{}\dots,~{}a_{ij}^{(n_{ij})}\right\}~{}~{}\&~{}~{}\theta=\bigcup_{\{i,j\}}\theta^{ij}.\end{array}$$

$$\begin{array}[]{l}\displaystyle\left[\mathcal{M}_{\theta}\right]_{l,m}=\intop_{-\infty}^{\infty}\mbox{Tr}\left(\phi_{e}^{-1}\left(j\omega\right)\frac{\partial G(j\omega)}{\partial\theta_{l}}\phi_{u}(j\omega)\frac{\partial G(j\omega)}{\partial\theta_{m}}^{*}\right)d\omega,~{}~{}1\leq l,m\leq N,\\ \displaystyle\hskip 51.21504pt\mbox{where,~{}~{}}\phi_{u}\left(j\omega\right)=\left[\phi_{ij}(j\omega)\right]_{\left\{1\leq i,j\leq r\right\}}\mbox{~{}~{}and~{}~{}}\phi_{e}(j\omega)=\Sigma.\end{array}$$

$${\left[{\mathcal{M}}_{\theta }\right]}_{l,m}=\underset{-\infty }{\overset{\infty }{\int }}\text{mbox}Tr\left(\frac{\partial G(j\omega )}{\partial {\theta }_l}{\varphi }_u(j\omega ){\frac{\partial G(j\omega )}{\partial {\theta }_m}}^{\ast }\right)d\omega ,1\le l,m\le N.$$

$$\begin{array}[]{l}\displaystyle\intop_{-\infty}^{\infty}\mbox{Tr}\left(\phi_{u}(j\omega)\right)d\omega\leq K_{u}\mbox{~{}~{}and~{}~{}}\intop_{-\infty}^{\infty}\mbox{Tr}\left(G(j\omega)\phi_{u}(j\omega)G(j\omega)^{*}\right)d\omega\leq K_{y},\end{array}$$

$$\begin{array}[]{l}\displaystyle\phi^{*}=\mbox{arg~{}}\max_{\phi_{u}}~{}\mbox{log}\left(\mbox{det}\left(\mathcal{M}_{\theta}\right)\right)\mbox{~{}{or}~{}}-\mbox{Tr}\left(\mathcal{M}_{\theta}^{-1}\right)\mbox{~{}{or}~{}}\sigma_{\min}\left(\mathcal{M}_{\theta}\right)\mbox{~{}{or}~{}}\mbox{Tr}\left(\mathcal{M}_{\theta}\right)\\ \mbox{subject to}\\ \displaystyle\hskip 42.67912pt\intop_{-1}^{1}\mbox{Tr}\left(\phi_{u}(j\omega)\right)d\omega\leq K_{u}\\ \displaystyle\hskip 42.67912pt\intop_{-1}^{1}\mbox{Tr}\left(G(j\omega)\phi_{u}(j\omega)G(j\omega)^{*}\right)d\omega\leq K_{y}.\end{array}$$

$${\varphi }_u(\omega )=\sum _{k=-M}^M{c}_k{e}^{-j\omega k}.$$

$$\begin{array}[]{l}\hskip 85.35826ptH_{m}(j\omega)=W\mathcal{H}\\ \mbox{where,~{}~{}}W=\left[\begin{array}[]{cccc}I_{r\times r}&w_{1}I_{r\times r}&\dots&w_{m}I_{r\times r}\end{array}\right]\\ \mbox{and~{}~{}}\mathcal{H}=\left[\begin{array}[]{cccc}H_{0}&H_{1}&\dots&H_{m}\end{array}\right]^{\top},~{}~{}\forall i~{}H_{i}\in\mathcal{R}^{r\times r}.\end{array}$$

$${\varphi }_u^m(j\omega )=\frac{1}{2}W(j\omega )\mathcal{H}{\mathcal{H}}^{\top }W(j\omega {)}^{\ast }.$$

$$\underset{-\infty }{\overset{\infty }{\int }}\text{mbox}Tr\left(\frac{1}{2}W(j\omega )\mathcal{H}{\mathcal{H}}^{\top }W(j\omega {)}^{\ast }\right)d\omega \le K_u.$$

$$\begin{array}[]{l}\intop_{-\infty}^{\infty}\mbox{Tr}\left(\frac{1}{2}W(j\omega)\mathcal{H}\mathcal{H}^{\top}W(j\omega)^{*}\right)d\omega=\mbox{Tr}\left(\mathcal{H}^{\top}\left(\intop_{-1}^{1}\frac{1}{2}W(j\omega)^{*}W(j\omega)d\omega\right)\mathcal{H}\right)\\ \hskip 156.49014pt=\mbox{Tr}\left(\mathcal{H}^{\top}\mathcal{W}\mathcal{H}\right)\\ \hskip 156.49014pt=h^{\top}P_{I}h,\\ \hskip 56.9055pt\mbox{where,~{}~{}}h=\mbox{vec}(\mathcal{H})\mbox{~{}~{}and~{}~{}}P_{I}=\left[\begin{array}[]{cccc}\mathcal{W}\\ &\ddots\\ &&\mathcal{W}\end{array}\right].\end{array}$$

$$h^{\top }{P}_{I}h\le K_u.$$

$$\underset{-\infty }{\overset{\infty }{\int }}\text{mbox}Tr\left(G(j\omega )W(j\omega )\mathcal{H}{\mathcal{H}}^{\top }W(j\omega {)}^{\ast }G(j\omega {)}^{\ast }\right)d\omega \le K_y-p.$$

$$\begin{array}[]{l}\hskip 142.26378pt\displaystyle h^{\top}P_{O}h\leq K_{y},\mbox{~{}where}\\ \displaystyle P_{O}=\left[\begin{array}[]{cccc}\mathcal{W}_{G}\\ &\ddots\\ &&\mathcal{W}_{G}\end{array}\right]\displaystyle\mbox{~{}~{}and~{}~{}}\mathcal{W}_{G}=\frac{1}{2}\intop_{-1}^{1}W(j\omega)^{*}G(j\omega)^{*}G(j\omega)W(j\omega)d\omega.\end{array}$$

$$\begin{array}[]{l}\displaystyle\left[\mathcal{M}_{\theta}\right]_{l,m}=\intop_{-\infty}^{\infty}\mbox{Tr}\left(\frac{\partial G(j\omega)}{\partial\theta_{l}}\phi_{u}(j\omega)\frac{\partial G(j\omega)}{\partial\theta_{m}}^{*}\right)d\omega=\mbox{Tr}\left(\left(\mathcal{H}^{\top}\tilde{G}_{l}\right)\left(\mathcal{H}^{\top}\tilde{G}_{m}\right)^{\top}\right),\end{array}$$

$$\begin{array}[]{l}\displaystyle\max_{h\in\mathcal{R}^{mr}}~{}\Phi(\mathcal{M}_{\theta}(h))\\ \displaystyle\mbox{subject to}\\ \displaystyle\hskip 28.45274pth^{\top}P_{I}h\leq K_{u}\\ \displaystyle\hskip 28.45274pth^{\top}P_{O}h\leq K_{y}.\end{array}$$

$$\begin{array}[]{l}x^{\top}M_{i,j}x\mbox{~{}~{}is converted to~{}~{}}x^{\top}M_{i,j}y\mbox{~{}~{}for all $i,j$}.\end{array}$$

$$\begin{array}[]{l}\displaystyle\max_{x\in\mathcal{R}^{N}}~{}\Phi(\mathcal{M}(x,y)),~{}~{}\mbox{where~{}}\mathcal{M}_{ij}(x,y)=x^{\top}M_{i,j}y\mbox{~{}and~{}}M_{i,j}\in\mathcal{R}^{N\times N}\\ \hskip 85.35826pt~{}\&~{}\Phi(.)=\mbox{logdet}(.)|\mbox{$\sigma_{\min}$}(.)|\mbox{-trace-inv}(.)|\mbox{trace}(.)\\ \displaystyle\mbox{subject to}\\ \displaystyle\hskip 56.9055ptx^{\top}P_{I}x\leq 1\\ \displaystyle\hskip 56.9055ptx^{\top}P_{O}x\leq 1\\ \\ \displaystyle\mbox{tol}=||x-y||_{2}\\ \displaystyle y=x\end{array}$$

$$\begin{array}[]{l}x^{\prime}Ax=\mbox{Tr}(AX),\mbox{~{}~{}where~{}~{}}X=xx^{\top}.\end{array}$$

$$\begin{array}[]{l}\displaystyle\max_{X\in\mathcal{R}^{N\times N}}~{}\Phi(\mathcal{M}_{\theta}(X)),~{}~{}\mbox{where~{}}\left[\mathcal{M}_{\theta}\right]_{ij}(X)=\mbox{Tr}(A^{ij}X)\\ \displaystyle\mbox{subject to}\\ \displaystyle\hskip 28.45274pt\mbox{trace}(P_{I}X)\leq 1\\ \displaystyle\hskip 28.45274pt\mbox{trace}(P_{O}X)\leq 1\\ \displaystyle\hskip 28.45274ptX\succeq 0\mbox{~{}~{}and~{}~{}}\mbox{Rank}(X)=1\end{array}$$

$$\begin{array}[]{l}\displaystyle\max_{X\in\mathcal{R}^{N\times N}}~{}\Phi(\mathcal{M}(X)),~{}~{}\mbox{where~{}}\mathcal{M}_{ij}(X)=\mbox{trace}(M_{i,j}X),~{}~{}M_{i,j}\in\mathcal{R}^{N\times N}\\ \hskip 85.35826pt~{}\&~{}\Phi(.)=\mbox{logdet}(.)|\mbox{$\sigma_{\min}$}(.)|\mbox{trace-inv}(.)|\mbox{trace}(.)\\ \displaystyle\mbox{subject to}\\ \displaystyle\hskip 56.9055pt\mbox{Tr}(P_{I}X)\leq 1\\ \displaystyle\hskip 56.9055pt\mbox{Tr}(P_{O}X)\leq 1\\ \displaystyle\hskip 56.9055ptX\succeq 0\end{array}$$

$$\begin{array}[]{l}\displaystyle\max_{h\in\mathcal{R}^{mr^{2}}}~{}~{}\Phi\left(\left[\begin{array}[]{cccc}h^{\top}M_{1,1}h&\dots&h^{\top}M_{1,N}h\\ \vdots&&\vdots\\ h^{\top}M_{N,1}h&\dots&h^{\top}M_{N,N}h\\ \end{array}\right]\right)\\ \displaystyle\mbox{subject to}\\ \displaystyle\hskip 56.9055pth^{\top}P_{I}h\leq K_{u}\\ \displaystyle\hskip 56.9055pth^{\top}P_{O}h\leq K_{y},\\ \end{array}$$

$$\begin{array}[]{l}Q(x)=\left[x^{\top}M_{1,1}x~{}\dots~{}x^{\top}M_{i,j}x~{}\dots~{}x^{\top}M_{N,N}x~{}~{}x^{\top}P_{I}x~{}~{}x^{\top}P_{O}x\right]^{\top};~{}~{}i\leq j.\end{array}$$

$$\underset{G}{\underbrace{\left[\begin{array}{cc}\frac{3.2}{s^2+1.6s+1.5}& \frac{3.1}{s^2+1.9s+1}\\ & \\ \frac{5.2}{s^2+0.5s+0.95}& \frac{1.5}{s^2+0.9s+0.3}\end{array}\right]}},\underset{G_{nom}}{\underbrace{\left[\begin{array}{cc}\frac{3.125}{s^2+1.5s+1.563}& \frac{3.188}{s^2+2s+1.063}\\ & \\ \frac{5.313}{s^2+0.5s+1.063}& \frac{1.563}{s^2+s+0.312}\end{array}\right]}}.$$

$$\underset{G_{\mathcal{D}}}{\underbrace{\left[\begin{array}{cc}\frac{3.21}{s^2+1.603s+1.497}& \frac{3.105}{s^2+1.897s+0.998}\\ & \\ \frac{5.073}{s^2+0.492s+0.951}& \frac{1.387}{s^2+0.856s+0.284}\end{array}\right]}},\underset{G_{\mathcal{A}}}{\underbrace{\left[\begin{array}{cc}\frac{3.555}{s^2+1.639s+1.501}& \frac{2.742}{s^2+1.863s+0.961}\\ & \\ \frac{5.252}{s^2+0.502s+0.952}& \frac{1.458}{s^2+0.884s+0.295}\end{array}\right]}},$$

$$\underset{G_{\mathcal{T}}}{\underbrace{\left[\begin{array}{cc}\frac{3.197}{s^2+1.606s+1.506}& \frac{3.082}{s^2+1.920s+0.999}\\ & \\ \frac{5.172}{s^2+0.498s+0.951}& \frac{1.425}{s^2+0.873s+0.292}\end{array}\right]}},\underset{G_{\mathcal{E}}}{\underbrace{\left[\begin{array}{cc}\frac{2.787}{s^2+1.502s+1.488}& \frac{3.144}{s^2+1.905s+1.006}\\ & \\ \frac{5.308}{s^2+0.499s+0.952}& \frac{1.473}{s^2+0.887s+0.296}\end{array}\right]}}.$$

$$0=x^{\ast }\mathcal{W}x$$

$$\small H(j\omega)=\left[\begin{array}[]{cc}H_{11}(j\omega)&H_{12}(j\omega)\\ H_{21}(j\omega)&H_{22}(j\omega)\end{array}\right]\mbox{~{}~{}and~{}~{}}\phi(j\omega)=\frac{1}{2}H(j\omega)H^{*}(j\omega)=\left[\begin{array}[]{cc}\phi_{11}(j\omega)&\phi_{12}(j\omega)\\ \phi_{21}(j\omega)&\phi_{22}(j\omega)\end{array}\right].$$

$$\begin{array}[]{l}\phi_{11}(j\omega)=\frac{1}{2}\left(\left\lvert H_{11}(j\omega)\right\rvert^{2}+\left\lvert H_{12}(j\omega)\right\rvert^{2}\right),\phi_{22}(j\omega)=\frac{1}{2}\left(\left\lvert H_{21}(j\omega)\right\rvert^{2}+\left\lvert H_{22}(j\omega)\right\rvert^{2}\right)\\ \hskip 28.45274pt\phi_{12}(j\omega)=\frac{1}{2}\left(H_{11}(j\omega)H_{21}^{*}(j\omega)+H_{12}(j\omega)H_{22}^{*}(j\omega)\right)=\phi_{21}(j\omega)^{*}.\end{array}$$

$${\mathcal{R}}_u(\tau )=\mathbb{E}\left\{u(t)u^{\ast }(t-\tau )\right\}=\frac{1}{2\pi }\underset{-{\omega }_c}{\overset{{\omega }_c}{\int }}{\varphi }_u(j\omega )e^{j\omega \tau }d\omega ,$$

Full text

Keywords: System Identification; Optimal Input Design; Spectral Decomposition; Fisher Information Matrix; Convex Optimization

1 Introduction

In this paper, the problem of designing a band-limited optimal input for identifying a continuous time linear multi-input multi-output (MIMO) system, subject to input and output power constraints, is considered. The linear system considered is given by

[TABLE]

where

[TABLE]

For all $t\in\mathcal{R}^{+}$, the input $u(t)\in\mathcal{R}^{r\times 1}$, the output $y(t)\in\mathcal{R}^{p\times 1}$ and the additive white noise $e(t)\sim\mathcal{N}\left(0,\Sigma\right)$ ($\Sigma\in\mathcal{S}^{+}_{p}$, the space of positive definite $p\times p$ matrices). It is assumed that the linear system under consideration is stable and the order of elements of the matrix $G(s;\theta)$ are specified. The aim of system identification, in general, is to estimate the set of parameters $\theta$, the nominal values of which are specified, by exciting the system with inputs and recording the outputs. The schematic in Figure 1 shows the general system identification framework.

The quality of an unbiased estimator is judged by a suitable functional defined on the covariance matrix of the estimated system. It turns out that the covariance matrix of an unbiased estimator of a linear system is a matrix multiple of the inverse of the Fisher Information Matrix (FIM). Suppose that the set of parameters $\theta$ is indexed and its cardinality is $N$, then the FIM for (1) is defined as (see Kumar et al. (2015)):

[TABLE]

For simplicity of analysis, it is assumed that $\Sigma$ is the identity matrix. This assumption simplifies the expression of FIM to:

[TABLE]

Since the FIM depends on the input power spectrum, it is obvious that the covariance matrix of the estimated system would depend on the input chosen for exciting the system. Therefore, the goal is to determine a power spectrum $\phi_{u}(j\omega)$ which is optimal with respect one of the following commonly used functionals defined as :

• •

$\mathcal{D}$-Optimality which maximizes the determinant of $\mathcal{M}_{\theta}$ (Sanchez et al. (2012); Kumar et al. (2015)),

• •

$\mathcal{A}$-Optimality which minimizes the trace of $\mathcal{M}^{-1}_{\theta}$ (Aoki and Staley (1970); Nahi and DE Wallis (1969)),

• •

$\mathcal{T}$-Optimality which maximizes the trace of $\mathcal{M}_{\theta}$ (Srivastava and Anderson (1974)), and

• •

$\mathcal{E}$-Optimality which maximizes the minimum eigenvalue of $\mathcal{M}_{\theta}$ (Cheng (1980); Constantine (1981)).

It is also natural to consider power constraints on the input signals and the output signals. The input power constraint may arise due to limitations of the excitation system. If the linear system has peak responses at some frequencies, the output power constraint would limit possible damages during excitation. The expressions for the input and the output power constraints are given by:

[TABLE]

where $K_{u}$ and $K_{y}$ are specified. It is assumed that the elements of the transfer function matrix are strictly proper. This assumption holds good in many practical systems and allows for the integral limits to be limited to a band of frequencies around zero. With the aforementioned definitions and constraints, the mathematical problem at hand is outlined in Figure 2.

Optimal input design for linear system identification has been studied for a long time and is still an active area of research (see Mehra (1974); Annergren et al. (2017)). Although not comprehensive, this paragraph would touch upon some of the seminal pieces of work in this area and discuss in brief, the methodologies used therein. In each of Barenthin et al. (2008); Kumar et al. (2015), the power spectrum is parameterized as

[TABLE]

The optimal input design problem is cast as a convex optimization problem in terms of the covariance matrices $\left\{c_{k}\right\}_{k=0}^{M}$. In particular, the Kalman-Yakubovich-Popov lemma and the Hamburger’s moment condition aid the formulation. The former ensures the positive definiteness of the power spectrum while the latter guarantees the existence of a valid solution. The idea of optimal input design has many applications. Ljung et al. (2011) is an excellent article which discusses the applications of system identification in communication systems, sensor networks and machine learning algorithms. In addition, automotive systems and chemical engineering processes also find the need for system identification and optimal input design (Alberer et al. (2011); Kumar et al. (2015)).

The way this paper differs from the earlier references is that the spectral density matrix is parameterized in terms of $H(j\omega)$, such that $\phi_{u}(j\omega)=\frac{1}{2}H(j\omega)H^{*}(j\omega)$. The matrix $H(j\omega)$ is expressed as a linear combination of a suitable truncated basis of $\mathcal{L}^{2}\left(-\omega_{\mbox{\small cut-off}},\omega_{\mbox{\small cut-off}}\right)$ and the coefficients become the optimization variables. Observe that this idea readily enforces positivity of $\phi_{u}(j\omega)$, which otherwise requires the Fourier basis and the Kalman-Yakubovich-Popov lemma. This also implies that the proposed method is independent of the choice of the basis. Indeed, the resulting optimization problem in this paper is non-convex in general. The optimization problem has an interesting structure which has been studied in the literature and has wide applications in signal processing and control systems. Moreover, solving this class of optimization problems still remains open.

The paper is organized as follows. The second section sets up the optimization problem, presents methods to compute lower and upper bounds, and discusses aspects related to the convexity of the problem. The third section presents simulation results, i.e., the input and the output waveforms, and the estimated parameters of the system. The concluding remarks are presented in the last section followed by references. The appendices to the paper present a couple of theorems and discuss the methodology of generating inputs from the optimal power spectrum.

2 The proposed method

This section discusses the formulation of the optimization problem, the non-convex nature of it and a methodology for obtaining a sub-optimal solution. As mentioned earlier, the input power spectrum is limited to a pre-defined band given by $D=[-\omega_{c},\omega_{c}]$. Now, each element of $\phi_{u}(j\omega)$ must belong to the space $C^{\infty}(D)$. Spectral decomposition theorem states that one can find $H(j\omega)$ belonging to $C^{1}(D)$ such that $\phi_{u}(j\omega)=H(j\omega)H^{*}(j\omega)$. With the knowledge that $C^{1}(D)$ is dense in $\mathcal{L}^{2}(D)$, one can as well approximate elements $H(j\omega)$ through a linear combination of a truncated basis of $\mathcal{L}^{2}(D)$. Therefore, each element of $H(j\omega)$ is assumed to the space of square integrable functions over $D$ denoted by $\mathcal{L}^{2}(D)$ (see Nair (2001)).

2.1 Setting up the optimization problem

As mentioned in the beginning of this section, is it assumed that each element of $H(j\omega)$ belongs to $\mathcal{L}^{2}(\mathcal{D})$. Thus, $H(j\omega)$ can be approximated with as a linear combination of a truncated basis for $\mathcal{L}^{2}(\mathcal{D})$. Suppose the truncated basis elements are given by the set $\left\{w_{1},w_{2},\dots,w_{m}\right\}$. Then $H(j\omega)$ is expressed as:

[TABLE]

The matrix $\mathcal{H}$ comprises of the coefficients which need to be determined by solving an optimization problem. As a consequence, the input power spectral density can be expressed as:

[TABLE]

Using (9), the constraint on the total input power is given by:

[TABLE]

Interchanging the integral operator and the trace operator, one obtains

[TABLE]

Since $\mathcal{W}$ is a positive definite matrix by (A), $P_{I}$ is also positive definite. In short, the total input power constraint given by (10) is equivalent to

[TABLE]

Note that $h\in\mathcal{R}^{mr^{2}}$. Similarly the total output power constraint given by

[TABLE]

can be rewritten as

[TABLE]

Note that, as in (A), it can be shown that $\mathcal{W}_{G}$ is positive definite and so is $P_{O}$. Finally, the FIM $\mathcal{M}_{\theta}$ can be rewritten as

[TABLE]

where $\displaystyle\tilde{G}_{l}=\intop_{-\infty}^{\infty}\frac{\partial G(j\omega)}{\partial\theta_{l}}W$. In other words, every element of the FIM is a quadratic in the coefficients $h$. Now assuming $\displaystyle\Phi(.)=\mbox{logdet}(.)$ or \mbox{\sigma_{\min}}(.) or $-\mbox{trace-inv}(.)$ or $\mbox{trace}(.)$, the optimization problem (6) is equivalent to

[TABLE]

2.2 Finding a suboptimal solution

To the best of the authors’ knowledge, the optimization problem (16) does not have an efficient algorithm. The problem, due to its quadratic structure, renders itself to bi-linear forms, as well as semi-definite convex relaxations. The method proposed here depends on these derived forms and has two steps in it: (i) find a lower bound to (16) and (ii) find an upper bound to (16). The first step, which obtains a lower bound to the optimization problem through a bi-linear form, is iterative and in each iteration, a convex problem is solved. The procedure is asymptotically convergent and produces a feasible solution. The second step find an upper bound to the optimal value of (16) through a convex relaxation. Hence, it can be computed efficiently with semi-definite programming. The solution corresponding to the lower bound is a used as a sub-optimal solution to (16). The sub-optimality of the solution is determined by the difference between the upper and lower bounds.

Lower bound: To obtain a lower bound, the quadratic form in each of the elements of $\mathcal{M}_{\theta}$ is converted into a bi-linear form. The input and output power constraints, which are convex, are retained as they are. In other words,

[TABLE]

The bi-linear optimization problem is described, for example, in Konno (1976). Note that once one of the variables ($x$ or $y$) is fixed to some vector, (16) becomes a convex optimization problem in the other variable. Starting from an arbitrary feasible point $y=y_{0}$, the convex optimization problem in $x$ is solved to obtain $x^{*}$. For the next iteration, $y$ is fixed to the vector the $x^{*}$. The resulting vector of each iteration acts as a lower bound to (16) since it is a feasible solution. Note that with each iteration, the value of the cost function can only increase and hence, the lower bound becomes tighter. The process asymptotically converge to a value but however, a discussion on the rate of convergence is out of the scope of this paper. The algorithm, also termed as the ‘hill-climbing’ procedure, is outlined in Figure 3.

Upper bound: The upper bound is determined by a standard semi-definite convex relaxation of (16). It can be shown that any homogeneous quadratic $x^{\prime}Ax$ can also be written as (see Boyd and Vandenberghe (2004))

[TABLE]

This transformation implies that (16) is equivalent to

[TABLE]

Dropping the Rank-1 in (20) results in a convex semi-definite optimization problem shown in Figure 4. Being a relaxation of the original optimization problem, the solution to (20) gives an upper bound to (16). Moreover, the bound can be obtained efficiently as the associated problem is convex. However, unlike the lower bound algorithm, the solution matrix $X^{*}$ need not be have a rank of 1, in which case deriving a feasible solution for (16) from it may not be possible.

Remark: It might be necessary to have more power constraints in the optimization problem. In particular, it might be necessary to constrain input and output powers at the individual ports, in addition to the total power constraints. Moreover, the power constraints may also be non-convex, particularly when the constraint is to deliver a minimum amount of input or output power. In all the cases mentioned here, the algorithm for the upper bound does not change, except for the increase in the number of constraints. The lower bound algorithm would need a few changes. Firstly, the constraints are also changed to their corresponding bi-linear forms. Secondly, the step “$y=x$” in Figure 3 needs to be replaced by “$y=\frac{x_{k}+x_{k-1}}{2}$”. In other words, $y$ is replaced by the average of the resulting vectors of the last two iterates.

2.3 A discussion on quadratic maps

This subsection discusses the convexity properties of the optimization problem given in (16). As mentioned earlier, every element of the FIM $\mathcal{M}_{\theta}(h)$ is a homogeneous quadratic in $h$. For clarity, the optimization problem is written in an expanded form as:

[TABLE]

where $M_{i,j}$’s are symmetric matrix of dimension $mr\times mr$. Moreover, by construction, the symmetric matrices are such that the FIM in the argument of $\Phi$ is positive semidefinite for all $h\in\mathcal{R}^{mr}$. The convexity of this optimization problem depends on the convexity of the image of $\mathcal{R}^{mr}$ under the quadratic map $\displaystyle Q:\mathcal{R}^{mr^{2}}\rightarrow\mathcal{R}^{\frac{N(N+1)}{2}+2}$ given by:

[TABLE]

It is assumed here that the total number of coefficients given by $mr$ is much larger than the number of parameters to be estimated given by $N$. Quadratic maps and the convexity of its image have been studied earlier in literature and some of the important results are pointed out here.

Ramana and Goldman (1994) were one of the first few to study, in detail, the convexity of images under quadratic transformations. They give necessary and sufficient condition for characterizing such maps. However, the condition are not checkable for general cases. Polyak (1998) proved that the image of the unit disk under any “definite” quadratic map $Q:\mathcal{R}^{n}\rightarrow\mathcal{R}^{3}$ is convex for $n\geq 3$. Sheriff (2013) presented results for “stable convexity” and discussed its applications in analyzing coupled systems. Dymarsky (2014) showed the convexity of the image of an open ball in $\mathcal{R}^{n}$ under quadratic maps. Since the dimension of the range space in this paper is much larger than 3, these results cannot be applied here. Beck (2009) presented some startling generalizations of earlier results on the subject. One of the results can be applied to obtain a convex relaxation of the optimization problem in this paper, but only at the cost of a large number of redundant variables. Hiriart-Urruty and Torki (2002) provided an excellent overview on the convexity aspects of the images of quadratic maps in their paper. These references point out that quadratic maps with convex images are rather exceptions than the rule.

2.4 Input generation and system identification

The methodology for generating inputs real inputs with the optimal power spectrum obtained from the optimization procedure is outlined in Appendix B. The methodology essentially discretizes the power spectrum, concatenates it with a random phase and performs inverse Fourier transform to obtain the inputs. The system is excited with the the designed optimal inputs and the outputs are recorded. The system identification is done using the Instrument Variables method, which is an unbiased estimator of the system parameters. For details on the system identification algorithm, the reader is directed to Söderström and Stoica (2002).

3 Simulations

In most practical applications, a system is fabricated with a desired/nominal set of parameters. But fabrication is seldom perfect and therefore, the actual parameters vary to a small extent around the desired locations. For optimal input design, the nominal system parameters are used for generating the optimal inputs, while the system with actual parameters generates the outputs when actuated with the optimal inputs. For simulations, a 2-input and 2-output system was chosen. The nominal transfer function matrix and the actual transfer function matrix are given by:

[TABLE]

The total input power and the total output power were constrained to a maximum of 1000 units. Chebyshev’s polynomials up to the $13^{\rm{th}}$ order were considered for parameterizing $H(j\omega)$. It was observed that lower bounds and the upper bounds matched in the $\mathcal{D}$, $\mathcal{A}$ and $\mathcal{E}$ optimality criterion, while they were off by less than $0.1\%$ in case of $\mathcal{T}$ optimality. It has been observed that the lower and upper bounds match in several instances. However, there is no mathematical proof for the observation yet. The optimal power spectrums obtained using the proposed method have been depicted in Figure 5. The inputs generated according to the optimal spectrums corresponding to the four different optimality criterion are shown in Figures 6, 7, 8. In addition, the corresponding outputs, corrupted with uncorrelated additive Gaussian white noise vectors, are also shown in Figures 6, 7, 8. The standard deviation of the measurement noise is assumed to be equal to one-tenth of the minimum of the root-mean-squared values of the outputs. The system identification was done using Instrument Variables method and the results corresponding to the different optimality criteria are given below.

[TABLE]

[TABLE]

4 Conclusion

In this paper, a distinct paradigm to design a band-limited optimal inputs for identifying a linear continuous time MIMO system was presented. Unlike conventional schemes where the power spectral density is parametrized with covariance matrices, this paper used spectral decomposition. The input spectra $\phi_{u}(j\omega)$ was factored as $\frac{1}{2}H(j\omega)H^{*}(j\omega)$ and $H(j\omega)$ was approximated by a truncated basis for $\mathcal{L}^{2}([-\omega_{c},\omega_{c}])$, where $\omega_{c}$ is the cut-off frequency. The coefficients of the finite linear combination thus became the optimization parameters. The parameterization ensured the positivity of the power spectral density matrix, irrespective of the basis chosen for $\mathcal{L}^{2}(\mathcal{D})$. The optimization problem of determining the optimal coefficients for $\mathcal{D},~{}\mathcal{A}$, $\mathcal{T}$ or $\mathcal{E}$ optimality criteria turned out to be a non-standard optimization problem. A discussion on the convexity aspect of the problem was also presented. A lower bound for the optimal value was obtained by first converting the optimization problem into a bi-linear form and then solving it using an iterative procedure. In each of the iterations, a convex problem was solved which made the algorithm computationally efficient. An upper bound for the optimal value was obtained by first converting the optimization problem into a rank constrained semi-definite problem. The Rank-1 constraint was then dropped and the resulting convex semi-definite program was solved using standard routines. If the upper bound turned out to be a Rank-1 solution, the vector defining its column space was chosen as the optimal solution; else the solution corresponding to the lower bound was chosen. The sub-optimality of the final solution was determined by the difference in the bounds. Interestingly, the bounds matched in many instances of simulations. With the optimal spectral density, the corresponding inputs were generated. The linear MIMO system was stimulated with these inputs and the outputs, corrupted with measurement noise, were recorded. The Instrument Variables method was used to estimate the system parameters. The estimated system was observed to be in close agreement with the original system.

As for future directions of work, several interesting questions can be posed. The authors would like to point to a few: (i) The choice of basis for $\mathcal{L}^{2}(\mathcal{D})$ and its effect on the obtained optimality, and (ii) the ratio of the lower bound to the upper bound, in general.

Appendix A Theorems and proofs

Theorem. Suppose $\mathcal{W}=\frac{1}{2}\intop_{-1}^{1}W(j\omega)^{*}W(j\omega)d\omega$. Then $\mathcal{W}$ is an $r(m+1)\times r(m+1)$ positive definite matrix.

Proof.

It is clear that $\mathcal{W}$ is positive semidefinite. Suppose for some $x\neq 0$

[TABLE]

∎

Theorem. Suppose $\Phi$ is a matrix function defined on $[-\pi,\pi]$ such that (i) $\Phi(\theta)\in\mathbb{C}^{q\times q}$ is a Hermitian nonnegative definite matrix for all $\theta\in[-\pi,\pi]$, (ii) $\Phi(\theta)=\Phi^{\top}(-\theta)$ and (iii) $\Phi$ is integrable and admits a Fourier series. Then there exists a matrix function $\Psi$ defined on $[-\pi,\pi]$ such that $\Phi(\theta)=\Psi(\theta)\Psi^{*}(\theta)$. Moreover, $\Psi$ admits a Fourier series.

Proof.

See Wilson [1972]. ∎

Appendix B Design of inputs from optimal power spectrum

An algorithm for generating a $u(t)$ from $H(j\omega)$, which is obtained from the proposed method, is discussed in this section. For brevity, the case for generating a $2\times 1$ input is illustrated. The extension for generating a $r\times 1$ input is presented in Figure 10. Conventionally, the input is generated as a sample path of a stochastic process with the desired power spectral density. The same principle is adopted in this paper as well. To elaborate further, let

[TABLE]

Naturally,

[TABLE]

The goal is then to determine a real-valued random vector $u(t)=\left[u_{1}(t),u_{2}(t)\right]^{\top}$ such that

[TABLE]

where $\mathbb{E}\{.\}$ denotes the expected value. To that end, consider the one sided frequency functions $i.e.$ defined for $\omega\geq 0$,

[TABLE]

where the families of random variables $\left\{\Phi_{1}^{\omega},\omega\geq 0\right\}$ and $\left\{\Phi_{2}^{\omega},\omega\geq 0\right\}$ are such that

[TABLE]

For the input signals have to be real-valued, extend $U_{1}$ and $U_{2}$ to the negative imaginary axis as:

[TABLE]

Define the inputs as:

[TABLE]

The claim is that the power spectral density of $\left[u_{1}(t),u_{2}(t)\right]$ is equal to $\phi_{u}(j\omega)$. For proving the claim, note that by using Fubini’s Theorem Nair [2001],

[TABLE]

From (27) and (28), it can be deduced that

[TABLE]

Therefore,

[TABLE]

On similar lines, it follows that

[TABLE]

It is now easy to note from the relation (26) applied to (34) that the power spectral density of $[u_{1}(t)u_{2}(t)]^{\top}$ is equal to $\phi(j\omega)=H(j\omega)H^{*}(j\omega)$. Now, a closed form expression for the integrals in (30) is difficult to obtain and therefore it needs to be evaluated numerically by discretization of the frequency band. To this end, choose $\Phi_{1}^{\omega}$ and $\Phi_{2}^{\omega}$ to be independent Gaussian white noise vectors. Since $U_{1}(j\omega)$ is zero outside the interval $\left[-\omega_{c},\omega_{c}\right]$,

[TABLE]

With a positive integer $N$ such that $\frac{1}{N}\ll 1$, a Riemann sum approximation of (35) would be

[TABLE]

In a similar way, $\tilde{u_{2}}(t)$ can also be defined. The signal $\tilde{u}(t)=[\tilde{u_{1}}(t)~{}\tilde{u_{2}}(t)]$ is then used for simulations and computations. For a general $r$-input system, the procedure has been shown in Figure 10.