Bayesian Input Design for Linear Dynamical Model Discrimination

Abstract

A Bayesian design of the input signal for linear dynamical model discrimination has been proposed. The discrimination task is formulated as an estimation problem, where the estimated parameter indexes particular models. As the mutual information between the parameter and model output is difficult to calculate, its lower bound has been used as a utility function. The lower bound is then maximized under the signal energy constraint. Selection between two models and the small energy limit are analyzed first. The solution of these tasks is given by the eigenvector of a certain Hermitian matrix. Next, the large energy limit is discussed. It is proved that almost all (in the sense of the Lebesgue measure) high energy signals generate the maximum available information, provided that the impulse responses of the models are different. The first illustrative example shows that the optimal signal can significantly reduce error probability, compared to the commonly-used step or square signals. In the second example, Bayesian design is compared with classical average D-optimal design. It is shown that the Bayesian design is superior to D-optimal design, at least in this example. Some extensions of the method beyond linear and Gaussian models are briefly discussed.

1. Introduction

Discrimination of various dynamical models of the same process has a wide area of applications, especially in multiple-model fault detection and isolation [1,2,3,4] and, in many other estimation and control problems [5,6,7], it is necessary to choose the most likely dynamical model from a finite set. The discrimination task can be formulated as an estimation problem, where the estimated parameter $\theta$ indexes particular models. The problem can also be considered as a finite-dimensional approximation of more general identification tasks [8]. As the error probability or variance of the estimator of $\theta$ usually depends on the input signal, it is important to select a signal that minimizes error probability or maximizes a utility function that encodes the purpose of the experiment. Selection of an input signal that maximizes a utility function is strongly related to optimal experimental design [9].

Experimental design methods can be divided into classical and Bayesian. The classical methods, also called optimal experimental design, typically use various functionals of the Fisher information matrix as a utility function. These methods are widely described in the literature and work well if the model is linear in its parameters (see [8,10,11,12] and the review article [9]). Unfortunately, in typical identification tasks, the solution of model equation and the covariance depends non-linearly on $\theta$, even if the model equation is linear. This implies that the information matrix and the utility function depend on the parameter $\theta$ to be estimated. Therefore, only locally-optimal design can be obtained [13]. To obtain more robust methods, an averaging over the prior parameter distribution or minimax design [8,14] (Section 6.1) are commonly used, but these methods are not fully Bayesian.

Bayesian optimal design uses the utility function, a functional of the posterior distribution (see [15,16] and the review articles [13,17,18]). The most commonly used utility functions are mutual information between parameters and model output, Kullback-Leibler divergence between the prior and posterior distributions, and the determinant of the posterior covariance matrix [13,16]. In contrast to classical methods, in Bayesian design the utility function does not depend on the parameters to be estimated. Hence, the method can cope with non-linear problems. The utility function, which is suitable for model discrimination, is the error probability of the MAP estimator of $\theta$ [19]. Such a utility function is generally difficult to calculate (see [19]), but the result of Feder & Merhav [20] implies that the error probability of the MAP estimator is upper-bounded by some decreasing function of mutual information between $\theta$ and the output of the system. Hence, the maximization of mutual information creates the possibility of reducing the error probability, provided that appropriate estimator is used. However, the most serious problem that inhibits the development of this idea is great computational complexity in calculating the mutual information.

The main contribution of this article is a fully-Bayesian (in the terminology of [13]) method for finding an input signal that maximizes the mutual information between $\theta$ and the system output. Maximization of information or, equivalently, maximization of the output entropy has been proposed by many authors (see, e.g., [13,15,17,18,21,22,23]), but the mutual information is very hard to compute and the problem is often intractable. To overcome this serious difficulty, instead of mutual information the lower bound, given by Kolchinsky & Tracey [24], has been used. This is a pairwise-distance based entropy estimator and it it useful here, since it is differentiable, tight, and asymptotically reaches the maximum possible information (see [24] (Sections 3.2, 4, 6)). Maximization of such a lower bound, under the signal energy (i.e., the square of the signal norm) constraints, is much simpler, gives satisfactory solutions, and allows for practical implementation of the idea of maximizing information. This is illustrated with examples. Moreover, it is shown that, for certain cases, this problem reduces to a solution of a certain eigenproblem.

The article is organized as follows. In Section 2, the estimation task is formulated and the upper bound of the error probability and the lower bound of the mutual information are given. In Section 2.1, a selection between two models is discussed and an exact solution is given. Design of input signals with small energy, which is required in some applications, is described in Section 2.2. In Section 2.3, the large energy limit is discussed. An application to linear dynamical systems with unknown parameters is given in Section 3. An example of finding the most likely model among three stochastic models with different structures is given in Section 4. Comparison with classical D-optimal design is performed in Section 5. The article ends with conclusions and references.

2. Maximization of Mutual Information between the System Output and Parameter

Let us consider a family of linear models

$$Y=F_{\theta }U+Z,$$

(1)

where $\theta \in 1,2,...,r$, $Y,Z\in R^{n_Y}$, and $U\in R^{n_U}$. The matrices $F_{\theta }$ are bounded. The parameter $\theta$ is unknown. The prior distribution of $\theta$ is given by

$$P(\theta =i)=p_{0,i},i=1,...,r.$$

(2)

The random variable Z is conditionally normal (i.e., $p\left(Z\right|\theta )=N(Z,0,S_{\theta })$), where the covariance matrices $S_{\theta }$ are given a priori and $S_{\theta }>0$, for all $\theta$. The variable U is called the input signal. In all formulas below, the input signal U is a deterministic variable. The set of admissible signals is given by

$${\mathbf{S}}_{\varrho }=\{U\in R^{n_U};U^{T}U⩽\varrho \}.$$

(3)

Under these assumptions, and after applying Bayes rule:

$$\begin{array}{c}p\left(Y\right|U)=\sum _{\theta =1}^r{p}_{0,\theta }N(Y,F_{\theta }U,S_{\theta }),\hfill \end{array}$$

(4)

$$\begin{array}{c}p\left(Y\right|\theta ,U)=N(Y,F_{\theta }U,S_{\theta }),\hfill \end{array}$$

(5)

$$\begin{array}{c}p\left(\theta \right|Y,U)=\frac{p_{0,\theta }N(Y,F_{\theta }U,S_{\theta })}{{\displaystyle \sum _{j=1}^r}{p}_{0,j}N(Y,F_{j}U,S_j)}.\hfill \end{array}$$

(6)

The entropies of Y and $\theta$ and the conditional entropies are defined as

$$\begin{array}{c}H\left(\theta \right)=-\sum _{\theta =1}^r{p}_{0,\theta }ln{p}_{0,\theta },\hfill \end{array}$$

(7)

$$\begin{array}{c}H\left(\theta \right|Y,U)=-\int p\left(Y\right|U)\left(\sum _{\theta =1}^{r}p\left(\theta \right|Y,U)lnp\left(\theta \right|Y,U)\right)dY,\hfill \end{array}$$

(8)

$$\begin{array}{c}H\left(Y\right|U)=-\int p(Y\left|U\right)lnp\left(Y\right|U)dY,\hfill \end{array}$$

(9)

$$\begin{array}{c}H\left(Y\right|\theta )=\frac{1}{2}\sum _{\theta =1}^r{p}_{0,\theta }ln\left({\left(2\pi e\right)}^{n_y}\left|S_{\theta }\right|\right).\hfill \end{array}$$

(10)

The mutual information between $\theta$ and Y is defined as (see [25] (pp. 19, 250))

$$I(Y;\theta |U)=H(\theta \left|U\right)-H\left(\theta \right|Y,U)=H(Y\left|U\right)-H\left(Y\right|\theta ,U).$$

As $H\left(\theta \right|U)=H(\theta )$ and $H\left(Y\right|\theta ,U)=H(Y\left|\theta \right)$, then $I(Y;\theta |U)$ is given by

$$I(Y;\theta |U)=H(\theta )-H(\theta |Y,U)=H\left(Y\right|U)-H(Y\left|\theta \right).$$

(11)

The MAP estimator of $\theta$ is defined as

$$\widehat{\theta }(Y,U)=\arg \underset{\theta \in \{1,...,r\}}{max}p\left(\theta \right|Y,U).$$

The error probability of $\widehat{\theta }$ is given by (see [20])

$$P_e\left(U\right)=1-\int \left(\underset{\theta \in \{1,...,r\}}{max}p\left(\theta \right|Y,U)\right)p\left(Y\right|U)dY.$$

It follows from Fano’s inequality ([25] (p. 38)), that $P_e$ is lower bounded by an increasing function in $H\left(\theta \right|Y,U)$. Feder & Merhav [20] proved that $2P_e\left(U\right)⩽H\left(\theta \right|Y,U){log}_{2}e$. As $H\left(\theta \right|Y,U)=H(\theta )-I(Y;\theta \left|U\right)$ and $H\left(\theta \right)$ does not depend on U, then the maximization of $I(Y;\theta |U)$ creates the possibility of reducing $P_e$, and the optimal signal is given by

$$U^{*}\left(\varrho \right)=\arg \underset{U\in {\mathbf{S}}_{\varrho }}{max}I(Y;\theta |U).$$

(12)

To overcome the problems associated with the calculation of $I(Y;\theta |U)$, we will use its lower bound.

Lemma 1. 

(Information bounds). For all $U\in R^{n_U}$,

$$I_l\left(U\right)⩽I(Y;\theta |U)⩽H\left(\theta \right),$$

(13)

where

$$\begin{array}{c}{I}_l\left(U\right)=-\sum _{i=1}^r{p}_{0,i}ln\left(\sum _{j=1}^r{p}_{0,j}{e}^{-D_{i,j}\left(U\right)}\right),\hfill \end{array}$$

(14)

$$\begin{array}{c}{D}_{i,j}\left(U\right)=\frac{1}{4}{U}^T{Q}_{i,j}U+\frac{1}{2}ln|\frac{1}{2}(S_i+S_j)|-\frac{1}{4}ln|S_i\left|\right|S_j|,\mathrm{and}\hfill \end{array}$$

(15)

$$\begin{array}{c}{Q}_{i,j}={(F_i-F_j)}^T{(S_i+S_j)}^{-1}(F_i-F_j).\hfill \end{array}$$

(16)

Proof. 

According to (4), $p\left(Y\right|U)$ is finite Gaussian mixture. For such mixtures, the information bounds are known. A detailed proof, based on Chernoff $\alpha$-divergence, is given in [24] (Section 4). □

Lemma 2.

Let $\widehat{\theta }(Y,U)=\arg {max}_{\theta \in \{1,...,r\}}p\left(\theta \right|Y,U)$ be the MAP estimator of θ, and let $P_e\left(U\right)$ denote its error probability. There exists a continuous, increasing, and concave function $f:[0,H\left(\theta \right)]\to [0,1-r^{-1}]$, such that

$$P_e\left(U\right)⩽f(H\left(\theta \right)-I_l\left(U\right))⩽\frac{1}{2}(H\left(\theta \right)-I_l\left(U\right)){log}_{2}e.$$

(17)

Proof. 

Feder & Merhav [20] (see Theorem 1 and Equation (14)) proved that there exists an increasing, continuous, and convex function $\varphi :[0,1-r^{-1}]\to [0,H\left(\theta \right){log}_{2}e]$, such that

$$2P_e\left(U\right)⩽\varphi \left(P_e\left(U\right)\right)⩽H\left(\theta \right|Y,U){log}_{2}e.$$

(18)

As $H\left(\theta \right|Y,U)=H(\theta )-I(Y;\theta \left|U\right)$ and $I_l\left(U\right)⩽I(Y;\theta |U)$, then $\varphi \left(P_e\left(U\right)\right)⩽(H\left(\theta \right)-I_l\left(U\right)){log}_{2}e$. The function $g={\varphi }^{-1}$ is increasing, continuous, concave, and it follows from (18) that $2g\left(\eta \right)⩽\eta$. Hence, $P_e\left(U\right)⩽{\varphi }^{-1}\left((H\left(\theta \right)-I_l\left(U\right)){log}_{2}e\right)=g\left((H\left(\theta \right)-I_l\left(U\right)){log}_{2}e\right)⩽\frac{1}{2}(H\left(\theta \right)-I_l\left(U\right)){log}_{2}e$. Taking $f\left(\eta \right)=g\left(\eta {log}_{2}e\right)$ we obtain the result. □

Now, the approximate solution of (12) is given by

$$U^{*}\left(\varrho \right)=\arg \underset{U\in {\mathbf{S}}_{\varrho }}{max}{I}_l\left(U\right).$$

(19)

As $I_l$ is smooth and ${\mathbf{S}}_{\varrho }$ is compact, (19) is well-defined.

2.1. Selection between Two Models

Suppose that $\theta$ takes only two values, 1 and 2, with prior probabilities $p_{0,1}$ and $p_{0,2}=1-p_{0,1}$, respectively. It’s easy to check, by direct calculation, that

$$e^{-I_l\left(U\right)}={(p_{0,1}+p_{0,2}{e}^{-D_{1,2}\left(U\right)})}^{p_{0,1}}{(p_{0,1}{e}^{-D_{1,2}\left(U\right)}+p_{0,2})}^{p_{0,2}}.$$

(20)

Equation (20) implies that the maximization of $I_l$ is equivalent to the maximization of $D_{1,2}$. On the basis of (15), we have the following optimization: task

$$\underset{U^{T}U⩽\varrho }{max}{U}^T{Q}_{1,2}U.$$

(21)

The solution of (21) is the eigenvector of $Q_{1,2}$ corresponding to its largest eigenvalue; that is,

$$Q_{1,2}{U}^{*}={\lambda }_{max}\left(Q_{1,2}\right)U^{*},\left|\right|U^{*}{\left|\right|}^2=\varrho .$$

(22)

2.2. Small Energy Limit

In many practical applications, the energy of an excitation signal must be small. The second order Taylor expansion of (14) gives

$$I_l\left(U\right)=\frac{1}{4}{U}^{T}QU-\sum _{i=1}^r{p}_{0,i}ln\sum _{j=1}^r{\alpha }_{i,j}+{o\left(\right|\left|U\right||}^2),$$

(23)

where

$$Q=\sum _{i=1}^r{p}_{0,i}\left(\frac{\sum _{j=1}^r{\alpha }_{i,j}{Q}_{i,j}}{\sum _{j=1}^r{\alpha }_{i,j}}\right),\mathrm{and}$$

(24)

$${\alpha }_{i,j}=p_{0,j}{e}^{-D_{i,j}\left(0\right)}.$$

(25)

If the value of $\varrho$ (see (3)) is small, then the last term in (23) can be omitted. As the second term does not depend on U, we have the following optimization task:

$$\underset{U^{T}U⩽\varrho }{max}{U}^{T}QU.$$

(26)

The solution of (26) is the eigenvector of Q corresponding to its largest eigenvalue.

2.3. Large Energy Limit

We will investigate asymptotic behaviour of $I(Y;\theta |U)$ when $\left|\right|U\left|\right|\to \infty$. On the basis of Lemma 1, the condition

$$\underset{i\ne j}{min}{U}^T{Q}_{i,j}U>0$$

(27)

guarantees that ${lim}_{\varrho \to \infty }I(Y;\theta |\varrho U)=H\left(\theta \right)$. It is also possible that ${lim}_{\varrho \to \infty }I(Y;\theta |\varrho U)<H\left(\theta \right)$, for some U. Such signals are weakly informative and they cannot generate the maximum information, even if their amplitude tends to infinity. Let ${\mathbf{S}}_1$ denote the unit ball in $R^{n_U}$ and let $\mu$ be the Lebesgue measure on ${\mathbf{S}}_1$. The set of weakly informative signals is defined as

$$\Omega =\{U\in {\mathbf{S}}_1:\underset{\varrho \to \infty }{lim}I(Y;\theta |\varrho U)<H\left(\theta \right)\}.$$

(28)

Theorem 1.

$\mu (\Omega )=0$ if and only if $F_i\ne F_j$, for all $i\ne j$.

Proof. 

$\Leftarrow :$ On the basis of Lemma 1 and (28), $\Omega ={\cup }_{i\ne j}{\Omega }_{i,j}$, where ${\Omega }_{i,j}=\{\xi \in {\mathbf{S}}_1:{\xi }^T{Q}_{i,j}\xi =0\}$. Since $S_i+S_j$ is positive-definite and $F_i\ne F_j$ then, on the basis of (16), the matrix $Q_{i,j}$ has at least one positive eigenvalue. Hence, $\mu \left({\Omega }_{i,j}\right)=0$ and $\mu (\Omega )={\displaystyle \sum _{i\ne j}}\mu \left({\Omega }_{i,j}\right)=0$. $\Rightarrow :$ The condition $\mu (\Omega )=0$ implies that $\mu \left({\Omega }_{i,j}\right)=0$. Hence, $Q_{i,j}$ has at least one positive eigenvalue, which is possible only if $F_i\ne F_j$. □

As a conclusion, we have the following result.

Theorem 2.

Let $\widehat{\theta }(Y,U)=\arg {max}_{\theta \in \{1,...,r\}}p\left(\theta \right|Y,U)$, be the MAP estimator of θ and let $P_e\left(U\right)$ denote its error probability. If $F_i\ne F_j$ for all $i\ne j$, then, for any $ϵ>0$, there exists a number $\varrho >0$ and a signal $U\in {\mathbf{S}}_{\varrho }$, such that $P_e\left(U\right)<ϵ$.

Proof. 

By the assumption, and from Theorem 1, the set ${\mathbf{S}}_1\setminus \Omega$ is non-empty. If $U\in {\mathbf{S}}_1\setminus \Omega$, then ${min}_{i\ne j}{U}^T{Q}_{i,j}U>0$ and, from Lemma 1, we get ${lim}_{\varrho \to \infty }{I}_l\left(\varrho U\right)={lim}_{\varrho \to \infty }I(Y;\theta |\varrho U)=H\left(\theta \right)$. Now, Lemma 2 implies that $2P_e\left(U\right)⩽(H\left(\theta \right)-I_l\left(\varrho U\right)){log}_{2}e<ϵ$, for sufficiently large $\varrho$. □

3. Application to Linear Dynamical Systems

Consider, now, the family of linear systems

$$\begin{array}{c}{x}_{k+1}=A_{\theta }{x}_k+B_{\theta }{u}_k+G_{\theta }{w}_k,k=0,1,2,...,N-1,\hfill \end{array}$$

(29)

$$\begin{array}{c}{y}_k=C_{\theta }{x}_k+D_{\theta }{v}_k,k=1,2,...,N,\hfill \end{array}$$

(30)

where the prior distribution of $\theta$ is given by (2) and $x_k\in R^n,y_k\in R^m,w_k\in R^{n_w},v_k\in R^m$, $w_k\sim N(0,I_{n_w}),$ and $v_k\sim N(0,I_m)$. The variables $w_0,...,w_{N-1},v_1,...,v_N$ are mutually independent. The initial condition is zero. The solution of (29) with initial condition $x_0=0$ has the form

$$x_k=\sum _{i=0}^{k-1}{A}_{\theta }^{k-i-1}{B}_{\theta }{u}_i+\sum _{i=0}^{k-1}{A}_{\theta }^{k-i-1}{G}_{\theta }{w}_i.$$

(31)

If we denote $X=\mathrm{col}(x_1,...,x_N)$, $Y=\mathrm{col}(y_1,...,y_N)$, $U=\mathrm{col}(u_0,...,u_{N-1})$, $W=\mathrm{col}(w_0,...,w_{N-1})$, and $V=\mathrm{col}(v_1,...,v_N)$, then (31) and (30) can be rewritten as

$$\begin{array}{c}X={\mathcal{B}}_{\theta }U+{\mathcal{G}}_{\theta }W,\mathrm{and}\hfill \end{array}$$

(32)

$$\begin{array}{c}Y={\mathcal{C}}_{\theta }X+{\mathcal{D}}_{\theta }V,\hfill \end{array}$$

(33)

where the matrices ${\mathcal{B}}_{\theta }$, ${\mathcal{G}}_{\theta }$, ${\mathcal{C}}_{\theta }$, and ${\mathcal{D}}_{\theta }$ follow forms (30) and (31). The variables W and V are independent, where $W\sim N(0,I_{N{n}_w})$ and $V\sim N(0,I_{Nm})$. Substituting (32) into (33) we get Equation (1), where $F_{\theta }={\mathcal{C}}_{\theta }{\mathcal{B}}_{\theta }$, $Z={\mathcal{C}}_{\theta }{\mathcal{G}}_{\theta }W+{\mathcal{D}}_{\theta }V$. The conditional density of Z has the form $p\left(Z\right|\theta )=N(Z,0,S_{\theta })$, where the covariance matrix is given by

$$S_{\theta }={\mathcal{D}}_{\theta }{\mathcal{D}}_{\theta }^T+{\mathcal{C}}_{\theta }{\mathcal{G}}_{\theta }{\mathcal{G}}_{\theta }^T{\mathcal{C}}_{\theta }^T.$$

(34)

Hence, the results of Section 2 can be applied to the dynamical system (29) and (30).

4. Example

In some fault detection and isolation problems [3,4], there is a need to determine which of the known models of the process is the most adequate. It is, therefore, important to find a signal that emphasizes the differences between these various models. As an example of this type of problem, let us consider three stochastic continuous-time models:

$$dx=(A_{\theta }x+B_{\theta }u)dt+G_{\theta }dw,$$

(35)

where $\theta \in \{1,2,3\}$, $x\left(t\right)\in R^{\theta }$, $x\left(0\right)=0$, $u\left(t\right),w\left(t\right)\in R$, w is a standard Wiener process, and

$$\begin{array}{c}{A}_1=-1,B_1=1,G_1=0.05,\hfill \end{array}$$

(36)

$$\begin{array}{c}{A}_2=\left[\begin{array}{cc}0& 1\\ -3& -2.5\end{array}\right],B_2=\left[\begin{array}{c}0\\ 3\end{array}\right],G_2=\left[\begin{array}{c}0\\ 0.05\end{array}\right],\hfill \end{array}$$

(37)

$$\begin{array}{c}{A}_3=\left[\begin{array}{ccc}0& 1& 0\\ -3& -3.5& 1\\ 0& 0& -10\end{array}\right],B_3=\left[\begin{array}{c}0\\ 0\\ 30\end{array}\right],G_3=\left[\begin{array}{c}0\\ 0\\ 0.05\end{array}\right].\hfill \end{array}$$

(38)

The step responses of these models are similar and they are difficult to experimentally distinguish from each other if the noise level is significant. The observation equation has the form

$$y_k=x_1\left(t_k\right)+0.05v_k,k=1,2,...,N,$$

(39)

where $v_k\sim N(0,1)$, $t_k=k{T}_0$, $T_0=0.1$, and $x_1$ is the first component of $x\left(t\right)$. If $x_k=x\left(t_k\right)$ and $u\left(t\right)=u_k,t\in [t_{k-1},t_k)$, then, after discretization, the state $x_k$ and the output $y_k$ are described by (29) and (30), with appropriate matrices $A_{\theta }$, $B_{\theta }$, $C_{\theta }$, $G_{\theta }$, and $D_{\theta }$. The matrices $F_{\theta }$ and $S_{\theta }$ are calculated by using (31)–(34). Let us observe that, although the orders of the systems are different, the size of both $F_{\theta }$ and $S_{\theta }$ is always $N\times N$. We are interested in the maximization of $I_l\left(U\right)$. The solutions of (19) and (26) with a uniform prior, $\varrho =N$, and $N=200$ steps, are shown in the upper part of Figure 1. The step responses and the optimal responses are shown in the bottom part of Figure 1.

Let us observe that, in contrast to the step signal, the optimal signal clearly distinguishes the systems—although the energy of all input signals was the same and equal to N.

Let $U_{st},U_{sq}\in \partial {\mathbf{S}}_1$ denote the normalized step and square signal with period of three, respectively, and let $U^{*}\left(\varrho \right)$ denote the optimal signal. To check the validity of the results, the error probabilities $P_e\left(\varrho U_{st}\right)$, $P_e\left(\varrho U_{sq}\right)$, and $P_e\left(U^{*}\left(\varrho \right)\right)$ were estimated by Monte Carlo simulation with ${10}^6$ trials and $N=50$ steps. The results are shown in Figure 2. It was observed that the optimal signal gives an error probability several thousand times smaller than the step or square signal with the same energy. The second observation is that $P_e\left(\varrho U_{sq}\right)$ initially increased with $\varrho$. To explain this, let us note that Inequality (17) does not guarantee that $P_e$ is decreasing function of $\varrho$. Hence, it is possible that $P_e$ increases in certain directions, although Theorem 2 guarantees that $P_e$ tends to zero, provided that signal norm tends to infinity.

5. Comparison with the Average D-Optimal Design

Classical methods of signal design for parameter identification use various functionals of the Fisher information matrix as a utility function. One of the most popular is D-optimal design, which consists of finding a signal that maximizes the determinant of the information matrix (see [10,11,12] and the review article [9]). These methods are well-suited to models that are linear in their parameters. Unfortunately, in typical identification and discrimination tasks, the output is a non-linear function of the parameters and the information matrix depends on unknown parameters to be identified. One of the possibilities for avoiding this problem is the averaging of the utility function over the prior parameter distribution. This method is called average D-optimal design (see [14,26] and [9] (Sections 5.3.5 and 6), for details). The Bayesian design, described in the previous sections, will be compared with the average D-optimal design. To that end, let us consider a finite family of linear models (see also [12] (pp. 91–93))

$$y_k=\frac{b_{\theta }{z}^{-1}}{1-a_{\theta }{z}^{-1}}{u}_k+{\sigma }_v{v}_k,$$

(40)

where $\theta \in \{1,2,3,4\}$, $a_{\theta }=0.6+0.1(\theta -1)$, $b_{\theta }=1-a_{\theta }$, ${\sigma }_v=0.1$, and $v_k\sim N(0,1)$. The prior distribution of $\theta$ is uniform (i.e., $p_{0,\theta }=0.25$). The state space representation of (40) has the form

$$\begin{array}{c}{x}_{k+1}=a_{\theta }{x}_k+b_{\theta }{u}_k,\hfill \end{array}$$

(41)

$$\begin{array}{c}{y}_k=x_k+{\sigma }_v{v}_k,\hfill \end{array}$$

(42)

which is consistent with (29) and (30). The Fisher information matrix is given by

$$M_F(\theta ,U)=\frac{1}{N{\sigma }_v^2}\sum _{k=1}^N{d}_k{d}_k^T,$$

(43)

where $d_k={({\xi }_k,{\eta }_k)}^T$ and ${\xi }_k=\frac{\partial y_k}{\partial a_{\theta }}$, ${\eta }_k=\frac{\partial y_k}{\partial b_{\theta }}$, denote the sensitivity of the output $y_k$ to changes in parameters a and b, respectively. The derivatives ${\xi }_k$ and ${\eta }_k$ fulfil the sensitivity equations

$$\begin{array}{c}{\xi }_k=2a_{\theta }{\xi }_{k-1}-a_{\theta }^2{\xi }_{k-2}+b_{\theta }{u}_{k-2},\hfill \end{array}$$

(44)

$$\begin{array}{c}{\eta }_k=a_{\theta }{\eta }_{k-1}+u_{k-1},k=1,2,...,N,\hfill \end{array}$$

(45)

with zero initial conditions. The average D-optimal design consists in finding a signal U that maximizes the expectation of the determinant of the information matrix (see [9] (Sections 5.3.5 and 6), [11,14] (Chapter 6), and [12] for details). Hence, the utility function to be maximized has the form

$$J\left(U\right)=\sum _{\theta =1}^4{p}_{0,\theta }\left|M_F(\theta ,U)\right|,$$

(46)

with the energy constraints given by (3). Maximization of the utility function (46) has been performed for various signal energies and the error probability of the MAP estimator was estimated by Monte Carlo with ${10}^5$ trials. The same procedure was repeated using Bayesian design for (41) and (42). The results are shown in Figure 3. The error rate of Bayesian method is significantly smaller when compared to D-optimal design, at least in this example. In particular, the signal shown in the upper-right part of Figure 3 gives an error probability approximately three times smaller than that of D-optimal signal, although the energy of both signals was the same.

6. Possible Extensions of the Results

In this section, we will briefly discuss some possible extensions of the results to an infinite set of parameters and beyond linear and Gaussian models.

6.1. Non-Linear Models

Although the article refers to linear models, it is possible to extend the results to non-linear models of the form

$$Y=F_{\theta }\left(U\right)+Z,$$

(47)

where the conditional density of variable Z is given by

$$p\left(Z\right|\theta ,U)=N(Z,0,S_{\theta }\left(U\right))$$

(48)

and $S_{\theta }\left(U\right)>0$, for all $U\in U_{ad}$, $\theta \in \{1,...,r\}$. Under these assumptions, the density of Y still remains a Gaussian mixture and the information lower bound takes the form

$$I_l\left(U\right)=-\sum _{i=1}^r{p}_{0,i}ln\left(\sum _{j=1}^r{p}_{0,j}{e}^{-D_{i,j}\left(U\right)}\right),$$

(49)

where

$$\begin{array}{c}{D}_{i,j}\left(U\right)=\frac{1}{4}{(F_i\left(U\right)-F_j\left(U\right))}^T{(S_i\left(U\right)+S_j\left(U\right))}^{-1}(F_i\left(U\right)-F_j\left(U\right))+\hfill \\ +\frac{1}{2}ln|\frac{1}{2}(S_i\left(U\right)+S_j\left(U\right))|-\frac{1}{4}ln|S_i\left(U\right)\left|\right|S_j\left(U\right)|.\hfill \end{array}$$

(50)

6.2. Non-Gaussian Models

If $p\left(Z\right|\theta ,U)$ is non-Gaussian distribution, then it is possible, on the basis of Equation (10) in [24], to construct an information lower bound of the form

$$I_l\left(U\right)=-\sum _{i=1}^r{p}_{0,i}ln\left(\sum _{j=1}^r{p}_{0,j}{e}^{-C_{\alpha }(p_i\left|\right|p_j)}\right),$$

(51)

where

$$C_{\alpha }(p_i\left|\right|p_j)=-ln\int p{\left(Z\right|i,U)}^{\alpha }p{\left(Z\right|j,U)}^{1-\alpha }dZ$$

(52)

is the Chernoff $\alpha$-divergence and $\alpha \in [0,1]$. Unfortunately, calculation of (52) is difficult if $n_Y$ is large.

6.3. Infinite Set of Parameters

Let us consider following model:

$$Y=F\left(\theta \right)U+Z,$$

(53)

where $p\left(Z\right|\theta )=N(Z,0,S\left(\theta \right))$ and $\theta \in R^p$. If we assume that prior density of $\theta$ is Gaussian; that is,

$$p_0\left(\theta \right)=N(\theta ,m_{\theta },S_{\theta }),S_{\theta }>0,$$

(54)

then

$$p\left(Y\right)=\int p_0\left(\theta \right)N(Y,F\left(\theta \right)U,S\left(\theta \right)),$$

(55)

and $p\left(Y\right)$ can be approximated by a finite Gaussian mixture

$$p\left(Y\right|U)=\int p_0\left(\theta \right)N(Y,F\left(\theta \right)U,S\left(\theta \right))d\theta \approx \sum _{j=1}^{N_a}{p}_{0,j}N(Y,F\left({\theta }_j\right)U,S\left({\theta }_j\right)),$$

(56)

where $p_{0,j}⩾0$ and $\sum _{j=1}^{N_a}{p}_{0,j}=1$. It’s possible to calculate weights and nodes in (56) by using multidimensional quadrature. If $n_{\theta }$ is large, then an appropriate sparse grid should be used. To illustrate the method, we will show only a simple, second-order quadrature with $2n_{\theta }$ points.

Lemma 3.

The approximate value of the integral $J\left(f\right)=\int N(\theta ,m_{\theta },S_{\theta })f\left(\theta \right)d\theta$ is given by

$$J\left(f\right)\approx \frac{1}{2n_{\theta }}\sum _{j=1}^{2n_{\theta }}f\left({\theta }_j\right),$$

(57)

where

$${\theta }_{2i-1}=m_{\theta }-S_{\theta }^{0.5}{e}_i,{\theta }_{2i}=m_{\theta }+S_{\theta }^{0.5}{e}_i,i=1,...,n_{\theta }$$

(58)

and $e_i$ is $i^{th}$ basis vector. If $f\left(\theta \right)=\frac{1}{2}{\theta }^{T}A\theta +b^T\theta +c$, then equality holds in (57).

Proof. 

Direct calculation. □

Application of Lemma 3 to (56) gives $p_{0,j}=p_0={\left(2n_{\theta }\right)}^{-1}$, $N_a=2n_{\theta }$. Now, since (56) is a Gaussian mixture, the results of Section 2 can be utilized and the information lower bound takes the form

$$I_l\left(U\right)=-p_0\sum _{i=1}^{r}ln\left(\sum _{j=1}^r{e}^{-D_{i,j}\left(U\right)}\right)-ln{p}_0,$$

(59)

where $D_{i,j}$ and ${\theta }_j$ are given by (15), (16), and (58), respectively, and $F_j=F\left({\theta }_j\right)$, $S_j=S\left({\theta }_j\right)$. The approximate solution of (12) can be found by maximization of (59) with the constraints (3).

7. Discussion and Conclusions

An effective Bayesian design method for linear dynamical model discrimination has been developed. The discrimination task is formulated as an estimation problem, where the estimated parameter $\theta$ indexes particular models. To overcome the computational complexity, instead of $I(Y;\theta |U)$, its lower bound ($I_l\left(U\right)$), proposed by Kolchinsky & Tracey [24], was used as a utility function. This bound is especially useful, as it is differentiable, tight, and reaches the maximum available information $H\left(\theta \right)$. It has been proved, on the basis of the results of Feder & Merhav [20], that the error probability of the MAP estimator (see Lemma 2) is upper bounded by $\frac{1}{2}(H\left(\theta \right)-I_l\left(U\right)){log}_{2}e$. The maximization of $I_l\left(U\right)$ has been considered under the signal energy constraint, but other kinds of constraints can also be easily implemented. It was shown that the maximization of $I_l\left(U\right)$, in the case of two parameters, is equivalent to maximization of a quadratic form on the sphere (see also [3]). Next, the small energy limit was analyzed. It was proved that the solution is given by an eigenvector corresponding to maximal eigenvalue of some Hermitian matrix. This result can serve as starting point for numerical maximization of $I_l\left(U\right)$. If the energy of the signal tends to infinity, then almost all (in the sense of the Lebesgue measure) signals generate maximum information, provided that the impulse responses of the models are pairwise different. Under these conditions, it was proved that $P_e$ of the MAP estimator tends to zero.

An example of discrimination of three stochastic models with different structures was given. It is easy to observe, from Figure 1, that, in contrast to the step signal, the optimal signal clearly distinguished the systems, although the energy of both signals was the same. The $P_e$ of the MAP estimator was calculated by Monte Carlo simulation. It was observed that the square signal gave an error probability several thousand times greater than the optimal signal with the same energy. Hence, we conclude that the error probability and the accuracy of MAP estimator depends very strongly on the excitation signal. Although Theorem 2 implies that ${lim}_{\varrho \to \infty }{P}_e\left(\varrho U\right)=0$ for almost all U, there exist signals such that $P_e\left(\varrho U\right)$ is locally increasing. This is the case of a high-frequency square signal, as illustrated in Figure 2.

It was shown, in Section 5 (see Figure 3), that $P_e$ of the MAP estimator corresponding to Bayesian design was a few times smaller than $P_e$ generated by D-optimal design, at least in the analyzed example. This result suggests that Bayesian design can be applied to non-linear problems and that it is superior to classical D-optimal design.

Some extensions of the results to the infinite set of parameters and beyond linear and Gaussian assumptions were briefly discussed in Section 6. Extension to non-linear models seems to be easy, but the non-Gaussian case is difficult and a deeper analysis is required. The case of an infinite set of parameters was discussed in Section 6.3. It was shown that the measurement density can be approximated by a finite Gaussian mixture, after which the results of Section 2 could be directly applied. The general conclusion of this analysis is that the information bounds can be easily constructed, as long as the measurement density is a mixture of Gaussian distributions.

An analytical gradient formula must be provided for the effective numerical maximization of $I_l\left(U\right)$. The matrix inversions and the determinants in (15) and (16) should be calculated by SVD. To reduce the computational complexity and required memory resources, the symmetries appearing in (15) and (16), and the fact that $D_{i,i}=0$, should be utilized. The determinants in (15) and the matrices $F_i,S_i$ can be calculated off-line, but the matrices $Q_{i,j}$ may require too much memory if N and r are large. Therefore, $Q_{i,j}$ was calculated on-line.

Applications of the presented methods in dual control problems [5] are expected as part of our future work. An additional area in applications is the issue of automated testing of dynamical systems.