An Approximate Bayesian Approach to Optimal Input Signal Design for System Identification

1. Introduction

The design of informative input signals is the cornerstone of modern system identification. Without properly chosen excitation, even advanced estimation algorithms may fail to provide accurate parameter estimates, leading to unreliable prediction and control. Classical references such as [1,2,3] and reviews by [4,5,6], emphasize that identification is not only a matter of statistical estimation but also of experimental design, where the input signal determines the achievable information content. Optimal experimental design (OED) methods therefore play a crucial role in practical applications. Traditionally, OED has relied on the Fisher Information Matrix (FIM), with criteria such as D or A optimality widely used due to their computational efficiency and asymptotic guarantees [4,7]. However, FIM-based approaches are inherently local, relying on linearization and asymptotic normality. They may thus be fragile in scenarios with large model uncertainty or strongly non-linear stochastic dynamics.

A natural alternative is the Bayesian approach, which evaluates an experiment through its expected information gain, typically quantified by the mutual information between model parameters and observations [5,6,8,9,10,11]. Bayesian design has several advantages: it is globally valid over the parameter space, it naturally incorporates prior information, and it is applicable to non-linear, stochastic models. Its main drawback is computational intractability, since mutual information requires high-dimensional integration over parameters and observations. Since the general case is challenging and difficult to solve, in this article we focus on models that can be represented in the form

$$\mathit{Y}=\mathit{F}(\mathit{\theta },\mathit{U})+\mathit{Z},$$

(1)

where $\mathit{\theta }\in \{{\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_r\}$ is a parameter with prior distribution $P(\mathit{\theta }={\mathit{\theta }}_j)=p_{0,j}$. The noise $\mathit{Z}$ is conditionally normal, that is, $p\left(\mathit{Z}\right|\mathit{\theta })=\mathcal{N}(\mathit{Z},0,\mathbf{S}(\mathit{\theta },\mathit{U}))$, and $\mathit{U}$ is a design variable. For this class of models, the density of the observations $\mathit{Y}$ is a finite Gaussian mixture of the form $p\left(\mathit{Y}\right)={\displaystyle \sum _{j=1}^r}{p}_{0,j}\mathcal{N}(\mathit{Y},\mathit{F}({\mathit{\theta }}_j,\mathit{U}),\mathbf{S}(\mathit{\theta },\mathit{U}))$. Within these Gaussian mixtures, the mutual information between $\mathit{Y}$ and $\mathit{\theta }$ can be estimated from below by the effective and tractable pairwise-distance-based lower bound given by Kolchinsky and Tracey [12,13]. We maximize this bound to achieve an approximate optimal design parameter $\mathit{U}$ and then generalize the method to a parameter space of continuum cardinality. In particular, we show how to treat the Gaussian prior $p_0\left(\mathit{\theta }\right)=\mathcal{N}(\mathit{\theta },{\mathit{m}}_{\theta },{\mathbf{S}}_{\theta })$ and prior distributions with compact support.

In this work, we focus on quasi-linear systems, namely stochastic dynamical systems that are linear in the state variables but non-linear in the control variables. Such systems occur ubiquitously in science and engineering. In quantum mechanics, Hamiltonians and Lindblad dissipators depend on external control fields such as laser intensities, magnetic fields, or gate voltages, leading to a non-linear dependence on the control [14,15]. In chemical processes, flow rates directly determine reaction speeds in continuous stirred tank reactors [16,17]. In thermal plants, convection coefficients scale non-linearly with flow, giving rise to quasi-linear heat transfer dynamics [16]. This broad applicability makes quasi-linear systems a natural and important class for advanced input design. However, there is a notable lack of Bayesian design methods and software tools tailored to this class of systems. Therefore, in this paper, we address this gap by developing such a method. Specifically, we show that a finite sequence of observations generated by a quasi-linear system can always be expressed in the form of the model (1), and we provide an effective algorithm for calculating the lower bound on mutual information.

The study reported here constitutes a substantial and far-reaching extension of the initial results presented in [18], as well as related research reported in [19,20,21]. The article’s main contributions can be summarized as follows. We first introduce an Information-Theoretic Lower Bound (ITB) on the estimation error of any estimator, [22,23], and briefly discuss its relation to the Bayesian Cramér–Rao Bound (BCRB), [24,25,26] [23]. We conclude that maximizing mutual information is superior to maximizing Bayesian or classical Fisher information, which is consistent with the arguments presented in [5,8,10,11]. Building on this result, we introduce a novel Bayesian design method for model (1). To address the intractability and computational complexity of direct mutual information evaluation, we discretize the parameter space and maximize the Kolchinsky–Tracey lower bound [12,13], and subsequently extend this approach to a parameter space of continuum cardinality. We then focus on the application to linear and quasi-linear system identification. Since the information-theoretic bound requires inversion of large covariance matrices of the observations, whose dimensions grow linearly with the number of data points N (with $N\sim {10}^3$–${10}^6$ in applications), direct inversion is computationally problematic. To overcome this challenge, we develop an algorithm that reduces the dimension of the matrices needed for inversion by a factor of N, thus making the approach feasible for long experiments. The proposed Bayesian method is compared with the average D-optimal design [1,4,27], a semi-Bayesian method, and its advantages are demonstrated. The effectiveness of the method is further illustrated by four examples. The first two, intentionally elementary, highlights the effectiveness of our approach. The third and fourth examples are drawn from atomic sensor models, a domain where optimal input design is particularly critical. We analyze a controlled harmonic oscillator with stochastic disturbances as a paradigmatic atomic sensor model [28], and a complex magnetometer model with a non-linear dependence of the system matrices on the input [29]. In the latter case, we provide a simplified model, derive the optimal input, and demonstrate that it significantly outperforms the harmonic signal, which might otherwise be presumed optimal. Moreover, we show that the estimation error of the MAP estimator achieves the theoretical lower bound.

The article is organized as follows. Section 2 formulates the problem. Section 3 develops the approximate Bayesian solution for finite and infinite parameter spaces. Section 4 applies the method to quasi-linear systems. Section 5 compares the approach with classical design methods. Section 6 presents examples. Section 7 provides discussion and conclusions.

2. Formulation of the Problem

Let us consider a family of models

$$\mathit{Y}=\mathit{F}(\mathit{\theta },\mathit{U})+\mathit{Z},$$

(2)

where $\mathit{Y},\mathit{Z}\in {\mathbb{R}}^{n_Y}$, $\mathit{U}\in {\mathbb{R}}^{n_U}$, $\mathit{\theta }\in \Theta \subset {\mathbb{R}}^{n_{\theta }}$. The set $\Theta$ will be called parameter space. Parameter $\mathit{\theta }$ is unknown. Prior distribution of $\mathit{\theta }$ is denoted by $p_0$. Random variable $\mathit{Z}$ is conditionally normal i.e. $p\left(\mathit{Z}\right|\mathit{\theta })=\mathcal{N}(\mathit{Z},0,\mathbf{S}(\mathit{\theta },\mathit{U}))$, where $\mathbf{S}(\mathit{\theta },\mathit{U})\in {\mathbf{S}}^{+}\left(n_Y\right)$, for all $\mathit{\theta }\in \Theta$, $\mathit{U}\in {\mathbb{R}}^{n_U}$. Functions $\mathit{F}$ and $\mathbf{S}$ are smooth. The variable $\mathit{U}$ is called the design parameter or, in the context of dynamical systems, the input signal. The set of admissible signals is given by

$${\mathcal{U}}_{ad}=\{\mathit{U}\in {\mathbb{R}}^{n_U};|\mathit{U}-\tilde{\mathit{U}}|⩽\varrho \},$$

(3)

where $\tilde{\mathit{U}}$ is the given vector and $\varrho$ is maximal norm of the signal. We will also consider an alternative, and useful in some applications, definition of ${\mathcal{U}}_{ad}$:

$${\mathcal{U}}_{ad}=\{\mathit{U}\in {\mathbb{R}}^{n_U};{\mathit{U}}_{min}⩽\mathit{U}⩽{\mathit{U}}_{max}\},$$

(4)

where ${\mathit{U}}_{min},{\mathit{U}}_{max}\in {\mathbb{R}}^{n_U}$ are fixed vectors. Under these assumptions and after applying the Bayes rule we get the likelihood, evidence and posterior distribution of $\mathit{\theta }$:

$$p\left(\mathit{Y}\right|\mathit{\theta },\mathit{U})=\mathcal{N}(\mathit{Y},\mathit{F}(\mathit{\theta },\mathit{U}),\mathbf{S}(\mathit{\theta },\mathit{U})),$$

(5)

$$p\left(\mathit{Y}\right|\mathit{U})=\int p_0\left(\mathit{\theta }\right)\mathcal{N}(\mathit{Y},\mathit{F}(\mathit{\theta },\mathit{U}),\mathbf{S}(\mathit{\theta },\mathit{U}))d\mathit{\theta },$$

(6)

$$p\left(\mathit{\theta }\right|\mathit{Y},\mathit{U})={\displaystyle \frac{p_0\left(\mathit{\theta }\right)p\left(\mathit{Y}\right|\mathit{\theta },\mathit{U})}{p\left(\mathit{Y}\right|\mathit{U})}}.$$

(7)

The Minimum Mean Squared Error (MMSE) estimator of $\mathit{\theta }$ is then given by

$$\widehat{\mathit{\theta }}(\mathit{Y},\mathit{U})=\int \mathit{\theta }p\left(\mathit{\theta }\right|\mathit{Y},\mathit{U})d\mathit{\theta }.$$

(8)

To avoid the difficulties involved in calculating the integral (8), instead of the MSE, the Maximum a Posteriori (MAP) estimator is typically used. Taking the negative logarithm of both sides of (7) and omitting the terms independent on $\mathit{\theta }$, we get the following:

$$\mathcal{L}(\mathit{\theta },\mathit{Y},\mathit{U})=\frac{1}{2}{|\mathit{Y}-\mathit{F}(\mathit{\theta },\mathit{U})|}_{{\mathbf{S}}^{-1}(\mathit{\theta },\mathit{U})}^2+\frac{1}{2}ln\left|\mathit{S}(\mathit{\theta },\mathit{U})\right|-ln{p}_0\left(\mathit{\theta }\right).$$

(9)

Thus, the MAP estimator of $\mathit{\theta }$ is given by

$$\widehat{\mathit{\theta }}(\mathit{Y},\mathit{U})=\arg \underset{\mathit{\theta }\in \Theta }{min}\mathcal{L}(\mathit{\theta },\mathit{Y},\mathit{U}).$$

(10)

Estimators (8) or (10) may be biased, so the Cramér-Rao Bound (CRB) cannot be applied directly to them. However, a Bayesian version of the CRB exists and can be used to estimate the error of any, also biased, estimator [23,24,25,26]. Let us introduce Bayesian information (BI):

$${\mathbf{J}}_B={\mathbb{E}}_{p(\mathit{\theta },\mathit{Y}|\mathit{U})}\left[{\nabla }_{\mathit{\theta }}\mathcal{L}(\mathit{\theta },\mathit{Y},\mathit{U})\left({\nabla }_{\mathit{\theta }}\mathcal{L}{(\mathit{\theta },\mathit{Y},\mathit{U})}^{\mathsf{T}}\right)\right]={\mathbf{J}}_D+{\mathbf{J}}_P,$$

(11)

where

$${\mathbf{J}}_P={\mathbb{E}}_{p_0\left(\mathit{\theta }\right)}\left[{\nabla }_{\mathit{\theta }}ln{p}_0\left(\mathit{\theta }\right){({\nabla }_{\mathit{\theta }}ln{p}_0\left(\mathit{\theta }\right))}^{\mathsf{T}}\right],$$

(12)

is the fraction of BI associated to a prior and

$${\mathbf{J}}_D\left(\mathit{U}\right)={\mathbb{E}}_{p(\mathit{\theta },\mathit{Y}|\mathit{U})}\left[{\nabla }_{\mathit{\theta }}lnp\left(\mathit{Y}\right|\mathit{\theta },\mathit{U}){({\nabla }_{\mathit{\theta }}lnp\left(\mathit{Y}\right|\mathit{\theta },\mathit{U}))}^{\mathsf{T}}\right],$$

(13)

is part of BI provided by observations $\mathit{Y}$. The matrix ${\mathbf{J}}_D$ is the Bayesian equivalent of the Fisher information matrix. The Fisher information matrix can be recovered from (13) assuming $p_0\left(\mathit{\theta }\right)=\delta (\mathit{\theta }-{\mathit{\theta }}^0)$, where ${\mathit{\theta }}^0$ is the true value of the parameter. Let $\widehat{\mathit{\theta }}(\mathit{Y},\mathit{U})$ be any estimator of $\mathit{\theta }$. Assuming a sufficiently regular prior, it can be proven that

$$\mathbb{E}|\mathit{\theta }-\widehat{\mathit{\theta }}{(\mathit{Y},\mathit{U})|}^2⩾{\displaystyle \frac{n_{\theta }}{|{\mathbf{J}}_P+{\mathbf{J}}_D{\left(\mathbf{U}\right)|}^{1/n_{\theta }}}},$$

(14)

which is usually known as the Van Trees inequality or Bayesian Cramér-Rao Bound (BCRB) [23,24,30][inequality (2.9), p. 17]. Formulas (11)-(13) are well defined only under rather restrictive assumptions. In particular, the joint density $p(\mathit{Y},\mathit{\theta })$ must be differentiable with respect to $\mathit{\theta }$, and it must satisfy the regularity condition $\int {\nabla }_{\mathit{\theta }}p(\mathit{Y},\mathit{\theta })d\mathit{Y}=0$. Moreover, the prior and likelihood distributions must guarantee the existence of the expectations in (12) and (13). This excludes the uniform and many other useful prior distributions and considerably limits the applicability of inequality (14) (see [23] for details). Beyond inequality (14), there exists a large class of Bayesian bounds reported in [26]. Many of these bounds can serve as a utility function.

Probably one of the best design criterion is the Ziv-Zakai lower bound [31]. However, to compute this estimate, an additional, and rather complex, optimization sub-problem must be solved as shown in [31]. Therefore, due to the high computational complexity of the multivariate Ziv-Zakai bound, we will not consider it here. Given the application-oriented focus of this article and following the arguments presented in [5,6], we conclude that the entropy-based lower bound [22][p. 255], [23][Sect. 2.2, p.16-17], is a reasonable optimality criterion and provides slightly tighter estimates than the BCRB (cf. [31][Sect. V.D.]). To proceed, let us define the entropies of $\mathit{Y}$ and $\mathit{\theta }$ and the corresponding conditional entropies:

$$\begin{array}{cc}\hfill & H_Y\left(\mathit{U}\right)=\mathbb{E}(-lnp\left(\mathit{Y}\right|\mathit{U})),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & H_{\theta }=\mathbb{E}(-ln{p}_0\left(\mathit{\theta }\right)),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & H_{Y|\theta }\left(\mathit{U}\right)=\mathbb{E}(-lnp\left(\mathit{Y}\right|\mathit{\theta },\mathit{U}))=\frac{1}{2}\int p_0\left(\mathit{\theta }\right)ln\left({\left(2\pi e\right)}^{n_Y}\left|\mathbf{S}(\mathit{\theta },\mathit{U})\right|\right)d\mathit{\theta },\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & H_{\theta |Y}\left(\mathit{U}\right)=\mathbb{E}(-lnp\left(\mathit{\theta }\right|\mathit{Y},\mathit{U})).\hfill \end{array}$$

(18)

The mutual information (MI) between $\mathit{\theta }$ and $\mathit{Y}$ is defined as

$$I_{\mathit{\theta };Y}\left(\mathit{U}\right)=H_{\theta }-H_{\theta |Y}\left(\mathit{U}\right)=H_Y\left(\mathit{U}\right)-H_{Y|\theta }\left(\mathit{U}\right).$$

(19)

The following theorem establishes the ultimate limit of the estimation error expressed in terms of mutual information, demonstrating that the Bayesian Cramér–Rao Bound (BCRB) does not constitute a fundamental limit.

Theorem 1.

Let $\widehat{\mathit{\theta }}(\mathit{Y},\mathit{U})$ be any estimator of $\mathit{\theta }$. Then the following inequalities hold:

$$\mathbb{E}|\mathit{\theta }-\widehat{\mathit{\theta }}{(\mathit{Y},\mathit{U})|}^2⩾n_{\theta }{\left(2\pi e\right)}^{-1}{e}^{2n_{\theta }^{-1}(H_{\theta }-I_{\theta ;Y}\left(\mathit{U}\right))}⩾{\displaystyle \frac{n_{\theta }}{|{\mathbf{J}}_P+{\mathbf{J}}_D{\left(\mathit{U}\right)|}^{1/n_{\theta }}}}.$$

(20)

The proof is given in Appendix A. The first of the inequalities (20) will be called the Information-Theoretic Lower Bound (ITB). The last part of (20) is known as the Efroimovich inequality [23,25][inequality (2.7), Ch. 2.2, p. 16].

To determine the optimal signal, one may maximize the MI, the determinant of ${\mathbf{J}}_B$ or the determinant of FIM. The latter corresponds to the classical design methods outlined in Sect. Section 5. However, the right-hand side of (14) generally underestimates the estimation error, and large values of ${\mathbf{J}}_B$ do not necessarily guarantee a small error. We illustrate this problem in Appendix B. In contrast, the ITB (20), which plays a central role in our subsequent analysis, shows that maximizing $I_{\theta ;Y}\left(\mathit{U}\right)$ is essential for reducing the estimation error and provides a more fundamental criterion than maximizing either the Bayesian or classical Fisher information. In particular, in the context of optimal experimental design and input signal design for system identification, maximizing the mutual information between the parameters and the observations constitutes the most principled optimality criterion, as it directly quantifies the amount of knowledge gained about the parameters [5]. Accordingly, we define the optimal signal as the solution of the following optimization problem:

$${\mathit{U}}^{*}=\arg \underset{\mathit{U}\in {\mathcal{U}}_{ad}}{max}{I}_{\theta ;Y}\left(\mathit{U}\right).$$

(21)

As $I_{\theta ;Y}$ is smooth and ${\mathcal{U}}_{ad}$ is compact, then (21) is well defined. After solving the task, the MMSE, MAP, or any other estimator can be used to determine $\mathit{\theta }$.

Computing mutual information (MI) or its lower-bound remains a significant challenge. The analyses presented in the literature [5,6,10,11], show that this can be done in three main ways: 1) using Monte Carlo or nested Monte Carlo (MC) simulations, [5][Sect. 3.1.], [6]; 2) applying variational lower-bound (VLB) estimates of the MI, [5][Sect. 3.3.1], [6]; 3) utilizing existing, easily computable estimates of conditional entropy or the MI. Since we aim to numerically maximize MI, which also depends on the design parameter U, the procedure for calculating MI will be called by the optimization solver millions of times and must therefore be sufficiently fast. Consequently, although MC methods provide good estimates of MI, they are of limited use here. The VLB methods require simultaneous optimization of the variational distribution with respect to its parameters and the signal U, [5][Sect. 3.3.1 and 4.3.4], [6]. In addition, stochastic simulations are also used to compute the expected values in the VLB. As the goal of this article is to develop a simple design method that does not require hours of computation, we focus on the third option and use existing, easily computable lower-bounds of entropy or MI provided in [12] or [32].

3. Approximate Solutions

The optimization problem (21) becomes considerably more tractable when the parameter space $\Theta$ is finite. Consequently, we first derive an approximate solution for a finite set of parameters and subsequently extend this result to obtain an approximate solution for the case where $\Theta$ is an uncountable subset of ${\mathbb{R}}^{n_{\theta }}$.

3.1. Finite Parameter Space

Let $\Theta =\{{\mathit{\theta }}_1,...,{\mathit{\theta }}_r\},{\mathit{\theta }}_i\in {\mathbb{R}}^{n_{\theta }}$, ${\mathit{\theta }}_i\ne {\mathit{\theta }}_j$ and assume that $p_0$ is a discrete distribution of the form

$$p_0\left(\mathit{\theta }\right)=\sum _{j=1}^r{p}_{0,j}\delta (\mathit{\theta }-{\mathit{\theta }}_j),$$

(22)

where $p_{0,j}=P(\theta ={\theta }_j)$. Then, on the basis of (6) and (22), the density of $\mathit{Y}$ becomes a Gaussian mixture:

$$p\left(\mathit{Y}\right|\mathit{U})=\sum _{j=1}^r{p}_{0,j}\mathcal{N}(\mathit{Y},\mathit{F}({\mathit{\theta }}_j,\mathit{U}),\mathbf{S}({\mathit{\theta }}_j,\mathit{U})).$$

(23)

The application of the Bayes rule gives the posterior:

$$p\left({\mathit{\theta }}_j\right|\mathit{Y},\mathit{U})={\displaystyle \frac{p_{0,j}\mathcal{N}(\mathit{Y},\mathit{F}({\mathit{\theta }}_j,\mathit{U}),\mathbf{S}({\mathit{\theta }}_j,\mathit{U}))}{p\left(\mathit{Y}\right|\mathit{U})}}.$$

(24)

The discrete counterpart of formulas (15)-(19) takes the form:

$$\begin{array}{cc}\hfill & H_Y\left(\mathit{U}\right)=-\int p\left(\mathit{Y}\right|\mathit{U})lnp\left(\mathit{Y}\right|\mathit{U})d\mathit{Y},\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & H_{\theta }=-\sum _{j=1}^r{p}_{0,j}ln{p}_{0,j},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & H_{Y|\theta }\left(\mathit{U}\right)=\frac{1}{2}\sum _{j=1}^r{p}_{0,j}ln\left({\left(2\pi e\right)}^{n_Y}\left|\mathbf{S}({\mathit{\theta }}_j,\mathit{U})\right|\right),\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & H_{\theta |Y}\left(\mathit{U}\right)=\int p\left(\mathit{Y}\right|\mathit{U})\left(-\sum _{j=1}^{r}p\left({\mathit{\theta }}_j\right|\mathit{Y},\mathit{U})lnp\left({\mathit{\theta }}_j\right|\mathit{Y},\mathit{U})\right)d\mathit{Y},\hfill \end{array}$$

(28)

$$I_{\mathit{\theta };\mathit{Y}}\left(\mathit{U}\right)=H_{\theta }-H_{\theta |Y}\left(\mathit{U}\right)=H_Y\left(\mathit{U}\right)-H_{Y|\theta }\left(\mathit{U}\right).$$

(29)

The direct computation of mutual information (29) remains difficult and, in many cases, intractable. Hence, our central idea is to overcome this difficulty by replacing mutual information (29) with a computationally tractable and non-trivial lower bound. In particular, we observe that (23) is a finite Gaussian mixture. For such mixtures, one of the most effective lower bounds on $I_{\mathit{\theta }|\mathit{Y}}$ is the inequality introduced in [12].

Lemma 1.(Information bounds, [12]). For the Gaussian mixture (23) with $p_{0,j}=P(\theta ={\theta }_j)$, the following inequalities holds:

$$I_l\left(\mathit{U}\right)⩽I_{\theta ;Y}\left(\mathit{U}\right)⩽H_{\theta },$$

(30)

where

$$\begin{array}{cc}\hfill & I_l\left(\mathit{U}\right)=-\sum _{i=1}^r{p}_{0,i}ln\left(\sum _{j=1}^r{p}_{0,j}{e}^{-d_{i,j}\left(\mathit{U}\right)}\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & d_{i,j}\left(\mathit{U}\right)=\frac{1}{8}{\Delta }_{i,j}^{\mathsf{T}}{\left(\frac{1}{2}({\mathbf{S}}_i+{\mathbf{S}}_j)\right)}^{-1}{\Delta }_{i,j}+\frac{1}{2}ln\left|\frac{1}{2}({\mathbf{S}}_i+{\mathbf{S}}_j)\right|-\frac{1}{4}ln\left(|{\mathbf{S}}_i\left|\right|{\mathbf{S}}_j|\right),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & {\Delta }_{i,j}=\mathit{F}({\mathit{\theta }}_i,\mathit{U})-\mathit{F}({\mathit{\theta }}_j,\mathit{U}),{\mathbf{S}}_i=\mathbf{S}({\mathit{\theta }}_i,\mathit{U}),{\mathbf{S}}_j=\mathbf{S}({\mathit{\theta }}_j,\mathit{U}).\hfill \end{array}$$

(33)

Detailed proof is given in [13][Sect. IIIb, inequality (11), and Sect. IV, formula (15) with $\alpha =0.5$] and also in [12]. Now, the approximate solution of (21) is given by

$${\mathit{U}}^{*}=\arg \underset{\mathit{U}\in {\mathcal{U}}_{ad}}{max}{I}_l\left(\mathit{U}\right).$$

(34)

Since ${\mathcal{U}}_{ad}$ is compact and $I_l$ is smooth and bounded, the solution of (34) exists. We also note that in the case of two alternatives, that is, when $r=2$ in (31), we get

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

(35)

Accordingly, the optimal signal in this case arises as the solution of a somewhat simplified optimization problem:

$$\underset{\mathit{U}\in {\mathcal{U}}_{ad}}{max}{d}_{1,2}\left(\mathit{U}\right).$$

(36)

If the function $\mathit{F}$ in (2) is affine with respect to $\mathit{U}$ and the covariance $\mathbf{S}$ does not depend on $\mathit{U}$, then it follows from Lemma 1, that $d_{1,2}$ is a positive (semi-) definite quadratic form with respect to $\mathit{U}$. For constraints (4), we thus obtain a convex quadratic programming problem. In the case of constraints (3), one needs to find the minimum of $d_{1,2}$ on a closed ball in ${\mathbb{R}}^N$. This is also a convex problem, and it can be reduced to finding zeros of a scalar function [33][Thm. 4.1, p. 70 and Sect. 4.3]. Furthermore, if $\mathit{F}({\mathit{\theta }}_i,\mathit{U})={\mathit{F}}_i\mathit{U},i=1,2$, and the constraints are defined by (3), then the solution of (36) is the eigenvector of the matrix $\mathbf{Q}={({\mathit{F}}_1-{\mathit{F}}_2)}^T{({\mathbf{S}}_1+{\mathbf{S}}_2)}^{-1}({\mathit{F}}_1-{\mathit{F}}_2)$, corresponding to its largest eigenvalue (see [18] [Sect. 2.1], for details).

3.2. Infinite Parameter Space

Let us assume that $\Theta ={\mathbb{R}}^{n_{\theta }}$ and consider the Gaussian prior

$$p_0\left(\mathit{\theta }\right)=\mathcal{N}(\mathit{\theta },{\mathit{m}}_{\theta },{\mathbf{S}}_{\theta }),{\mathbf{S}}_{\theta }>0.$$

(37)

Then the integral in (6) can be approximated with a finite Gaussian mixture

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

(38)

where $p_{0,j}⩾0$ and ${\displaystyle \sum _{j=1}^{N_a}}{p}_{0,j}=1$. The weights $p_{0,j}$ and the nodes ${\mathit{\theta }}_j$ in (38) can be calculated by using the multidimensional Gauss-Hermite quadrature rule or any other suitable method.

The Gauss-Hermite quadrature of order p is exact for polynomials of degree at most $2p-1$. The approximation error of the integral (38) using the Gauss-Hermite quadrature of order p depends on the $2p$-th derivatives of the functions F and S. In the single-parameter case, with Gaussian prior $\mathcal{N}(\theta ,m_{\theta },{\sigma }_{\theta })$, the error estimate is given by the formula:

$$e⩽{\sigma }_{\theta }{\displaystyle \frac{4^{p}C}{p!}}\underset{\theta ,\mathit{U},\mathit{Y}}{sup}{\displaystyle \frac{d^{2p}}{d{\theta }^{2p}}}\mathcal{N}(\mathit{Y},\mathit{F}(\theta ,\mathit{U}),\mathbf{S}(\theta ,\mathit{U})),$$

(39)

where C is constant. The error tends to zero as $p\to \infty$ or ${\sigma }_{\theta }\to 0$. Therefore, the Gauss-Hermite approximation of the integral (38) is especially useful when the prior is narrow (small ${\sigma }_{\theta }$) or when the integrand, in the neighbourhood of the point $m_{\theta }$, can be well approximated by low-degree polynomials. To illustrate the method, we will show only a very simple second-order Gaussian quadrature rule with $2n_{\theta }$ points.

Lemma 2.

Approximate value of the integral $J\left(f\right)=\int \mathcal{N}(\mathit{\theta },{\mathit{m}}_{\theta },{\mathbf{S}}_{\theta })f\left(\mathit{\theta }\right)d\mathit{\theta }$ is given by

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

(40)

where

$${\mathit{\theta }}_{2i-1}={\mathit{m}}_{\theta }-{\mathbf{S}}_{\theta }^{0.5}\sqrt{n_{\theta }}{\mathit{e}}_i,{\mathit{\theta }}_{2i}={\mathit{m}}_{\theta }+{\mathbf{S}}_{\theta }^{0.5}\sqrt{n_{\theta }}{\mathit{e}}_i,i=1,...,n_{\theta }$$

(41)

and ${\mathit{e}}_i$ is $i^{th}$ basis vector in ${\mathbb{R}}^{n_{\theta }}$. If $f\left(\mathit{\theta }\right)=\frac{1}{2}{\mathit{\theta }}^{\mathsf{T}}\mathbf{A}\mathit{\theta }+{\mathit{b}}^{\mathsf{T}}\mathit{\theta }+c$, then the equality holds in (40).

Proof. 

Direct calculation. □

An analogous method can be used for prior distributions defined on compact subsets of ${\mathbb{R}}^{n_{\theta }}$ (e.g. an n-dimensional hypercube), but the formulas for nodes and weights in (38) will then change. For example, if $\theta$ is a scalar parameter and the prior distribution is uniform, that is, $p_0\left(\theta \right)=U[a,b]$, then, using a second-order Gauss-Legendre quadrature, the approximate value of the integral $\underset{a}{\overset{b}{\int }}{p}_0\left(\theta \right)f\left(\theta \right)d\theta$ is computed using the formula:

$$\underset{a}{\overset{b}{\int }}{p}_0\left(\theta \right)f\left(\theta \right)d\theta \approx p_{0,1}f\left({\theta }_1\right)+p_{0,2}f\left({\theta }_2\right),$$

(42)

where $p_{0,1}=p_{0,2}=0.5$ and

$${\theta }_1=\frac{1}{2}\left(a+b-{\displaystyle \frac{b-a}{\sqrt{3}}}\right),{\theta }_2=\frac{1}{2}\left(a+b+{\displaystyle \frac{b-a}{\sqrt{3}}}\right).$$

(43)

The formula (42) is exact for polynomials of degree 3. The error estimate is analogous to (39) and tends to zero when $b-a\to 0$. More general multidimensional formulas, integration methods and error estimates are given in [34,35,36].

The application of lemma 2 to (38) gives $N_a=2n_{\theta }$, $p_{0,j}={\left(2n_{\theta }\right)}^{-1}$. Now, since (38) is approximated by a Gaussian mixture, the results of Section 3.1 can be used. Based on Eqs. (31-33), (38) and Lemma 1, the information lower bound takes the form:

$$I_l\left(\mathit{U}\right)=-{\displaystyle \frac{1}{2n_{\theta }}}\sum _{i=1}^{2n_{\theta }}ln\left({\displaystyle \frac{1}{2n_{\theta }}}\sum _{j=1}^{2n_{\theta }}{e}^{-d_{i,j}\left(\mathit{U}\right)}\right),$$

(44)

where $d_{i,j}$ and ${\mathit{\theta }}_j$ are given by Eqs. (31-33) and (41) or (43). The approximate solution of (21) can be found by maximizing (44) with constraints (3) or (4).

4. Bayesian Input Signal Design in Quasi-Linear Control Systems

Consider the family of quasi-linear systems

$${\mathit{x}}_{k+1}=\mathbf{A}(\mathit{\theta },{\mathit{u}}_k){\mathit{x}}_k+\mathbf{B}(\mathit{\theta },{\mathit{u}}_k)+\mathbf{G}(\mathit{\theta },{\mathit{u}}_k){\mathit{w}}_k,$$

(45)

$${\mathit{y}}_k=\mathbf{C}{\mathit{x}}_k+{\mathit{v}}_k,$$

(46)

where $k=0,1,...,N$, $N⩾1$, ${\mathit{x}}_k\in {\mathbb{R}}^n,{\mathit{y}}_k\in {\mathbb{R}}^{n_y},{\mathit{w}}_k\in {\mathbb{R}}^{n_w},{\mathit{v}}_k\in {\mathbb{R}}^{n_y}$, ${\mathit{w}}_k\sim \mathcal{N}(0,\mathbb{I})$, ${\mathit{v}}_k\sim \mathcal{N}(0,{\mathbf{S}}_v),{\mathbf{S}}_v>0$. Variables ${\mathit{x}}_0,{\mathit{w}}_0,...,{\mathit{w}}_{N-1},{\mathit{v}}_0,...,{\mathit{v}}_N$ are mutually independent. The initial state ${\mathit{x}}_0$ is conditionally normal i.e. $p\left({\mathit{x}}_0\right|\mathit{\theta })=\mathcal{N}({\mathit{x}}_0,{\mathit{m}}_0^{-}\left(\mathit{\theta }\right),{\mathbf{S}}_0^{-}\left(\mathit{\theta }\right))$, where ${\mathit{m}}_0^{-}$, ${\mathbf{S}}_0^{-}$ are smooth and ${\mathbf{S}}_0^{-}\left(\mathit{\theta }\right)>0$, for all $\mathit{\theta }\in \Theta$. The joint prior distribution of the initial state ${\mathit{x}}_0$ and the parameter $\mathit{\theta }$ is given by $p_0({\mathit{x}}_0,\mathit{\theta })=p_0\left(\mathit{\theta }\right)\mathcal{N}({\mathit{x}}_0,{\mathit{m}}_0^{-}\left(\theta \right),{\mathbf{S}}_0^{-}\left(\mathit{\theta }\right))$. Let’s define ${\mathbf{A}}_k=\mathbf{A}(\mathit{\theta },{\mathit{u}}_k)$, ${\mathbf{B}}_k=\mathbf{B}(\mathit{\theta },{\mathit{u}}_k)$, ${\mathbf{G}}_k=\mathbf{G}(\mathit{\theta },{\mathit{u}}_k)$. Then the solution of (45) has the form:

$$\begin{array}{cc}\hfill & {\mathit{x}}_0=\mathbb{I}{\mathit{x}}_0,\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill & {\mathit{x}}_1={\mathbf{A}}_0{\mathit{x}}_0+{\mathbf{B}}_0+{\mathbf{G}}_0{\mathit{w}}_0,\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill & {\mathit{x}}_2={\mathbf{A}}_1{\mathit{x}}_1+{\mathbf{B}}_1+{\mathbf{G}}_1{\mathit{w}}_1={\mathbf{A}}_1{\mathbf{A}}_0{\mathit{x}}_0+{\mathbf{A}}_1{\mathbf{B}}_0+{\mathbf{B}}_1+{\mathbf{A}}_1{\mathbf{G}}_0{\mathit{w}}_0+{\mathbf{G}}_1{\mathit{w}}_1,\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill & ⋮\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill & {\mathit{x}}_N=\Phi (N,0){\mathit{x}}_0+\sum _{j=0}^{N-1}\Phi (N,j+1){\mathbf{B}}_j+\sum _{j=0}^{N-1}\Phi (N,j+1){\mathbf{G}}_j{\mathit{w}}_j,\hfill \end{array}$$

(51)

where $\Phi (n,n)=\mathbb{I}$ and

$$\Phi (n,j)=\prod _{i=1}^{n-j}{\mathbf{A}}_{n-i},j<n.$$

(52)

Now, if we denote $\mathit{X}=\mathrm{col}({\mathit{x}}_0,...,{\mathit{x}}_N)$, $\mathit{Y}=\mathrm{col}({\mathit{y}}_0,...,{\mathit{y}}_N)$, $\mathit{U}=\mathrm{col}({\mathit{u}}_0,...,{\mathit{u}}_{N-1})$, $\mathit{W}=\mathrm{col}({\mathit{w}}_0,...,{\mathit{w}}_{N-1})$, $\mathit{V}=\mathrm{col}({\mathit{v}}_0,...,{\mathit{v}}_N)$, we can rewrite Eqs. (47-51) and (46), in matrix-vector form

$$\mathit{X}=\mathcal{A}(\mathit{\theta },\mathit{U}){\mathit{x}}_0+\mathcal{B}(\mathit{\theta },\mathit{U})+\mathcal{G}(\mathit{\theta },\mathit{U})\mathit{W},$$

(53)

$$\mathit{Y}=\mathcal{C}\mathit{X}+\mathit{V},$$

(54)

where the matrices $\mathcal{A}$, $\mathcal{B}$, $\mathcal{G}$, $\mathcal{C}={\mathbb{I}}_{N+1}\otimes \mathbf{C}$ follow directly from Eqs. (47-51), (46), and $\mathit{W}\sim \mathcal{N}(0,{\mathbb{I}}_{N{n}_w})$, $\mathit{V}\sim \mathcal{N}(0,{\mathbb{I}}_{N+1}\otimes {\mathbf{S}}_v)$. Substituting (53) into (54) and taking into account that $p\left({\mathit{x}}_0\right|\mathit{\theta })=\mathcal{N}({\mathit{x}}_0,{\mathit{m}}_0^{-}\left(\mathit{\theta }\right),{\mathbf{S}}_0^{-}\left(\mathit{\theta }\right))$, we get

$$\begin{array}{cc}\hfill & \mathit{Y}=\mathcal{C}\mathcal{A}(\mathit{\theta },\mathit{U}){\mathit{m}}_0^{-}\left(\mathit{\theta }\right)+\mathcal{C}\mathcal{B}(\mathit{\theta },\mathit{U})+\mathit{Z},\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill & \mathit{Z}=\mathcal{C}\mathcal{A}(\mathit{\theta },\mathit{U})({\mathit{x}}_0-{\mathit{m}}_0^{-}\left(\mathit{\theta }\right))+\mathcal{C}\mathcal{G}(\mathit{\theta },\mathit{U})\mathit{W}+\mathit{V}.\hfill \end{array}$$

(56)

The conditional density of variable $\mathit{Z}$ has the form $p\left(\mathit{Z}\right|\mathit{\theta })=\mathcal{N}(\mathit{Z},0,\mathbf{S}(\mathit{\theta },\mathit{U}))$, where the covariance matrix $\mathbf{S}$ is given by

$$\mathbf{S}(\mathit{\theta },\mathit{U})=\mathcal{C}(\mathcal{A}(\mathit{\theta },\mathit{U}){\mathbf{S}}_0^{-}\left(\mathit{\theta }\right)\mathcal{A}{(\mathit{\theta },\mathit{U})}^{\mathsf{T}}+\mathcal{G}(\mathit{\theta },\mathit{U})\mathcal{G}{(\mathit{\theta },\mathit{U})}^{\mathsf{T}}){\mathcal{C}}^{\mathsf{T}}+{\mathbb{I}}_{N+1}\otimes {\mathbf{S}}_v.$$

(57)

Finally, if we define

$$\mathit{F}(\mathit{\theta },\mathit{U})=\mathcal{C}\mathcal{A}(\mathit{\theta },\mathit{U}){\mathit{m}}_0\left(\mathit{\theta }\right)+\mathcal{C}\mathcal{B}(\mathit{\theta },\mathit{U}),$$

(58)

we can rewrite (55) in the form $\mathit{Y}=\mathit{F}(\mathit{\theta },\mathit{U})+\mathit{Z}$, which is exactly the model (2). To find the optimal input signal, we maximize one of the criteria (31) or (44) with constraints (3) or (4).

With a large number of data (large N), calculating the inverse and determinant of a very large matrix $\mathbf{S}(\mathit{\theta },\mathit{U})$ in (9) and calculating the quantities $d_{i,j}\left(\mathit{U}\right)$ in Eqs. (31-33) is numerically ill-conditioned and requires special treatment. The algorithms below reduce the size of the matrices necessary to invert by a factor of $N+1$.

Lemma 3.

The following identities hold:

$$p\left(\mathit{Y}\right|\mathit{\theta },\mathit{U})=\prod _{k=0}^N\mathcal{N}({\mathit{y}}_k,\mathbf{C}{\mathit{m}}_k^{-}\left(\mathit{\theta }\right),{\Sigma }_k\left(\mathit{\theta }\right)),$$

(59)

$${|\mathit{Y}-\mathit{F}(\mathit{\theta },\mathit{U})|}_{\mathbf{S}{(\mathit{\theta },\mathit{U})}^{-1}}^2=\sum _{k=0}^N{|{\mathit{y}}_k-\mathbf{C}{\mathit{m}}_k^{-}\left(\mathit{\theta }\right)|}_{{\Sigma }_k^{-1}\left(\mathit{\theta }\right)}^2,$$

(60)

$$\left|\mathbf{S}(\mathit{\theta },\mathit{U})\right|=\prod _{k=0}^N\left|{\Sigma }_k\left(\mathit{\theta }\right)\right|,$$

(61)

$$\mathcal{L}(\mathit{\theta },\mathit{Y},\mathit{U})=\frac{1}{2}\sum _{k=0}^N\left(|{\mathit{y}}_{\mathit{k}}-\mathbf{C}{\mathit{m}}_k^{-}{\left(\mathit{\theta }\right)|}_{{\Sigma }_k^{-1}\left(\mathit{\theta }\right)}^2+ln\left|{\Sigma }_k\left(\theta \right)\right|\right)-ln{p}_0\left(\mathit{\theta }\right),$$

(62)

where $\mathcal{L}$ is given by (9) and ${\mathit{m}}_k^{-}$, ${\Sigma }_k$, are calculated recursively by the Kalman filter

$$\begin{array}{cc}\hfill & {\Sigma }_k\left(\mathit{\theta }\right)={\mathbf{S}}_v+\mathbf{C}{\mathbf{S}}_k^{-}\left(\mathit{\theta }\right){\mathbf{C}}^{\mathsf{T}},\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill & {\mathbf{L}}_k\left(\mathit{\theta }\right)={\mathbf{S}}_k^{-}\left(\mathit{\theta }\right){\mathbf{C}}^{\mathsf{T}}{\Sigma }_k^{-1}\left(\mathit{\theta }\right),\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill & {\mathit{m}}_k\left(\mathit{\theta }\right)={\mathit{m}}_k^{-}\left(\mathit{\theta }\right)+{\mathbf{L}}_k\left(\mathit{\theta }\right)({\mathit{y}}_k-\mathbf{C}{\mathit{m}}_k^{-}\left(\mathit{\theta }\right)),\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill & {\mathbf{S}}_k\left(\mathit{\theta }\right)={\mathbf{S}}_k^{-}\left(\mathit{\theta }\right)-{\mathbf{L}}_k\left(\mathit{\theta }\right){\Sigma }_k\left(\mathit{\theta }\right){\mathbf{L}}_k{\left(\mathit{\theta }\right)}^{\mathsf{T}},\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill & {\mathit{m}}_{k+1}^{-}\left(\mathit{\theta }\right)={\mathbf{A}}_k{\mathit{m}}_k\left(\mathit{\theta }\right)+{\mathbf{B}}_k,\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill & {\mathbf{S}}_{k+1}^{-}\left(\mathit{\theta }\right)={\mathbf{A}}_k{\mathbf{S}}_k\left(\mathit{\theta }\right){\mathbf{A}}_k^{\mathsf{T}}+{\mathbf{G}}_k{\mathbf{G}}_k^{\mathsf{T}},k=0,1...,N,\hfill \end{array}$$

(68)

with initial conditions ${\mathit{m}}_0^{-}\left(\mathit{\theta }\right)$, ${\mathbf{S}}_0^{-}\left(\mathit{\theta }\right)$.

The proof is given in Appendix A. The Eqs. (63-68), are, in fact, a family of discrete-time Kalman filters indexed by $\mathit{\theta }$. The first four formulas describe the correction step. The prediction step is given by the last two equations. The matrix ${\mathbf{L}}_k$ is the Kalman gain and ${\Sigma }_k$ is the covariance matrix of the output prediction error ${\mathit{ϵ}}_k={\mathit{y}}_k-\mathbf{C}{\mathit{m}}_k^{-}$.

Lemma 4.

Let us define

$${\tilde{\mathbf{A}}}_k=\left[\begin{array}{cc}\mathbf{A}({\mathit{\theta }}_i,{\mathit{u}}_k)& 0\\ 0& \mathbf{A}({\mathit{\theta }}_j,{\mathit{u}}_k)\end{array}\right],{\tilde{\mathbf{B}}}_k=\left[\begin{array}{c}\mathbf{B}({\mathit{\theta }}_i,{\mathit{u}}_k)\\ \mathbf{B}({\mathit{\theta }}_j,{\mathit{u}}_k)\end{array}\right],$$

(69)

$${\tilde{\mathbf{G}}}_k=\left[\begin{array}{cc}\mathbf{G}({\mathit{\theta }}_i,{\mathit{u}}_k)& 0\\ 0& \mathbf{G}({\mathit{\theta }}_j,{\mathit{u}}_k)\end{array}\right],\tilde{\mathbf{C}}={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}\mathbf{C}& -\mathbf{C}\end{array}\right]$$

(70)

and let

$$\begin{array}{cc}\hfill & {\tilde{\Sigma }}_k={\mathbf{S}}_v+\tilde{\mathbf{C}}{\tilde{\mathbf{S}}}_k^{-}{\tilde{\mathbf{C}}}^{\mathsf{T}},\hfill \end{array}$$

(71)

$$\begin{array}{cc}\hfill & {\tilde{\mathbf{L}}}_k={\tilde{\mathbf{S}}}_k^{-}{\tilde{\mathbf{C}}}^{\mathsf{T}}{\tilde{\Sigma }}_k^{-1},\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill & {\tilde{\mathbf{S}}}_k={\tilde{\mathbf{S}}}_k^{-}-{\tilde{\mathbf{L}}}_k{\tilde{\Sigma }}_k{\tilde{\mathbf{L}}}_k^{\mathsf{T}},\hfill \end{array}$$

(73)

$$\begin{array}{cc}\hfill & {\tilde{\mathit{m}}}_{k+1}^{-}={\tilde{\mathbf{A}}}_k(\mathbb{I}-{\tilde{\mathbf{L}}}_k\tilde{\mathbf{C}}){\tilde{\mathit{m}}}_k^{-}+{\tilde{\mathbf{B}}}_k,\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill & {\tilde{\mathbf{S}}}_{k+1}^{-}={\tilde{\mathbf{A}}}_k{\tilde{\mathbf{S}}}_k{\tilde{\mathbf{A}}}_k^{\mathsf{T}}+{\tilde{\mathbf{G}}}_k{\tilde{\mathbf{G}}}_k^{\mathsf{T}},k=0,1,...,N,\hfill \end{array}$$

(75)

with initial conditions

$${\tilde{\mathit{m}}}_0^{-}=\left[\begin{array}{c}{\mathit{m}}_0^{-}\left({\mathit{\theta }}_i\right)\\ {\mathit{m}}_0^{-}\left({\mathit{\theta }}_j\right)\end{array}\right],{\tilde{\mathbf{S}}}_0^{-}=\left[\begin{array}{cc}{\mathbf{S}}_0^{-}\left({\mathit{\theta }}_i\right)& 0\\ 0& {\mathbf{S}}_0^{-}\left({\mathit{\theta }}_j\right)\end{array}\right].$$

(76)

Then the quantity $d_{i,j}\left(\mathit{U}\right)$ in formula (31) is given by:

$$d_{i,j}\left(\mathit{U}\right)=\frac{1}{4}\sum _{k=0}^N|\tilde{\mathbf{C}}{\tilde{\mathit{m}}}_k^{-}{|}_{{\tilde{\Sigma }}_k^{-1}}^2+\frac{1}{2}\sum _{k=0}^{N}ln\left|{\tilde{\Sigma }}_k\right|-\frac{1}{4}ln\left(|{\mathbf{S}}_i\left|\right|{\mathbf{S}}_j|\right),$$

(77)

where $|{\mathbf{S}}_i|=|\mathbf{S}\left({\mathit{\theta }}_{\mathit{i}},\mathit{U}\right)|$, $|{\mathbf{S}}_j|=|\mathbf{S}\left({\mathit{\theta }}_{\mathit{j}},\mathit{U}\right)|$, are calculated accordingly to Lemma 3, Eq. (61).

The proof is given in Appendix A. Let us observe that, instead of calculating the inverse and determinant of the large matrices ${\mathbf{S}}_i$, ${\mathbf{S}}_j$, $\frac{1}{2}({\mathbf{S}}_i+{\mathbf{S}}_j)$, of dimension $(N+1)n_y$, we only need to calculate the determinants and inverses of the much smaller matrices ${\Sigma }_k$, ${\tilde{\Sigma }}_k$, whose dimension is $n_y$, which is usually a small number.

5. Comparison with Classical Methods of Input Signal Design

Classical methods for input signal design in system identification are primarily concerned with LTI state space or transfer function models (such as ARMAX) and are usually based on maximizing some functions of error covariance or Fisher information matrix. For the prediction error method (PEM) estimator, the asymptotic form of the error covariance matrix (or Fisher information) is well known, both in the time and frequency domains. In the time domain, the solution corresponds to a specific input signal, whereas in the frequency domain, the solution yields the optimal power spectral density of the input signal. Below, we provide a brief overview of these methods, following the methodology presented in [1,2][Ch.9. Sect. 9.3, 9.4] and [4][Sect. 6.1].

Consider the LTI, SISO system

$$\begin{array}{cc}\hfill & {\mathit{x}}_{k+1}=\mathbf{A}\left(\mathit{\theta }\right){\mathit{x}}_k+\mathbf{B}\left(\mathit{\theta }\right)u_k+\mathbf{G}\left(\mathit{\theta }\right){\mathit{w}}_k,\hfill \end{array}$$

(78)

$$\begin{array}{cc}\hfill & y_k=\mathbf{C}{\mathit{x}}_k+v_k,\hfill \end{array}$$

(79)

under the assumptions stated in Section 4. System (78), (79) is equivalent to the transfer function model

$$y_k=G(\mathit{\theta },z)u_k+H(\mathit{\theta },z)e_k,$$

(80)

where $e_k\sim \mathcal{N}(0,{\sigma }_e^2)$ is a sequence of mutually independent Gaussian variables. The filters G and H are determined by the formulas

$$G(\mathit{\theta },z)=\mathbf{C}{(z\mathbb{I}-\mathbf{A}\left(\mathit{\theta }\right))}^{-1}\mathbf{B}\left(\mathit{\theta }\right),H(\mathit{\theta },z)=1+\mathbf{C}{(z\mathbb{I}-\mathbf{A}\left(\mathit{\theta }\right))}^{-1}\mathbf{K}\left(\mathit{\theta }\right),$$

(81)

where the Kalman gain $\mathbf{K}\left(\mathit{\theta }\right)$ is given by

$$\mathbf{K}\left(\mathit{\theta }\right)=\mathbf{A}\left(\mathit{\theta }\right)\mathbf{S}\left(\mathit{\theta }\right){\mathbf{C}}^{\mathsf{T}}{(\mathbf{C}\mathbf{S}\left(\mathit{\theta }\right){\mathbf{C}}^{\mathsf{T}}+{\sigma }_v^2)}^{-1},$$

(82)

with a non-negative matrix $\mathbf{S}$ being a solution of the Riccati equation (cf. [2])

$$\mathbf{S}=\mathbf{A}\mathbf{S}{\mathbf{A}}^{\mathsf{T}}+\mathbf{G}{\mathbf{G}}^{\mathsf{T}}-\mathbf{A}\mathbf{S}{\mathbf{C}}^{\mathsf{T}}{(\mathbf{C}\mathbf{S}{\mathbf{C}}^{\mathsf{T}}+{\sigma }_v^2)}^{-1}\mathbf{C}\mathbf{S}{\mathbf{A}}^{\mathsf{T}}.$$

(83)

The prediction errors are given by the recurrence

$${ϵ}_k(\mathit{\theta },\mathit{Y},\mathit{U})=H^{-1}(\mathit{\theta },z)\left(y_k-G(\mathit{\theta },z)u_k\right).$$

(84)

The cost function used in the Prediction Error Method (PEM) is expressed as:

$$V(\mathit{\theta },\mathit{Y},\mathit{U})={\displaystyle \frac{1}{2N{\sigma }_e^2}}\sum _{k=1}^N{ϵ}_k^2(\mathit{\theta },\mathit{U}).$$

(85)

Minimization of (85) with reference to $\mathit{\theta }$ yields the PEM estimator

$$\widehat{\mathit{\theta }}(\mathit{Y},\mathit{U})=\arg \underset{\mathit{\theta }\in \Theta }{min}V(\mathit{\theta },\mathit{Y},\mathit{U}).$$

(86)

The above estimator, under rather weak identifiability conditions, is consistent, asymptotically normal, and efficient, i.e. it achieves the Cramér-Rao lower bound. Following the reasoning presented in [1,2][Ch. 9, Sect. 9.3 and 9.4] or [4][Sect. 6.1], we divide the parameter vector into two groups related to the parameters appearing in G and H, that is, $\mathit{\theta }=\mathrm{col}({\mathit{\theta }}_{\mathit{H}},{\mathit{\theta }}_{\mathbf{G}})$. The sensitivity of ${ϵ}_k$ to changes in ${\mathit{\theta }}_G$ is calculated recursively according to the following equations:

$${\mathit{\psi }}_k(\mathit{\theta },\mathit{U})=H^{-1}(\mathit{\theta },z)\left({\nabla }_{{\mathit{\theta }}_G}G(\mathit{\theta },z)\right)u_k=F_z(\mathit{\theta },z)u_k,$$

(87)

where ${\nabla }_{{\mathit{\theta }}_G}$ means differentiation only with respect to the parameters that occur in G. The information matrix, which is also the inverse of the error covariance ${\mathit{P}}_{{\mathit{\theta }}_G}$, is given by:

$$\mathit{M}(\mathit{\theta },\mathit{U})={\mathit{P}}_{{\mathit{\theta }}_G}^{-1}(\mathit{\theta },\mathit{U})={\mathit{R}}_e\left(\mathit{\theta }\right)+{\displaystyle \frac{1}{N{\sigma }_e^2}}\sum _{k=1}^N{\mathit{\psi }}_k(\mathit{\theta },\mathit{U}){\mathit{\psi }}_k{(\mathit{\theta },\mathit{U})}^{\mathsf{T}},$$

(88)

where ${\mathit{R}}_e$ does not depend on $\mathit{U}$. Using the D-optimal criterion, the optimal signal is given by maximization of $det\left(\mathit{M}({\mathit{\theta }}^0,\mathit{U})\right)$, where ${\mathit{\theta }}^0$ is the true value of the parameter. Since ${\mathit{\theta }}^0$ is unknown, one can use the prior distribution and maximize the average D-optimal criterion:

$$Q\left(\mathit{U}\right)={\mathbb{E}}_{p_0\left(\mathit{\theta }\right)}det\left(\mathit{M}(\mathit{\theta },\mathit{U})\right),$$

(89)

with constraints (3) or (4). The asymptotic error covariance can also be expressed in terms of the power spectral density of the input signal $u_k$. Let ${\Phi }_u$ denote the spectral density of $u_k$. As it was shown in [2][p. 291] or [4] [Sect. 6.2], we have:

$$\mathit{M}({\Phi }_u,{\mathit{\theta }}^0)={\mathit{P}}_{{\mathit{\theta }}_G}^{-1}({\Phi }_u,{\mathit{\theta }}^0)={\displaystyle \frac{N}{2\pi {\sigma }_e^2}}\underset{-\pi }{\overset{\pi }{\int }}{F}_z(e^{i\omega },{\mathit{\theta }}^0)F_z{(e^{-i\omega },{\mathit{\theta }}^0)}^{\mathsf{T}}{\Phi }_u\left(\omega \right)d\omega +{\mathit{R}}_e\left({\mathit{\theta }}^0\right),$$

(90)

where $F_z$ is defined by (87) and the term ${\mathit{R}}_e$ in (90) does not depend on ${\Phi }_u$. Similarly as before, the parameter-averaged determinant of the matrix $\mathit{M}$ is maximized with respect to ${\Phi }_u$, subject to the signal power and frequency constraints. Typically, the spectrum ${\Phi }_u$ is parametrized by a finite number of coefficients $c_k$, so that the resulting optimization problem is convex; see [37] for details. After performing spectral factorization of ${\Phi }_u$, a filter is obtained, whose input is white noise and whose output yields the optimal signal $u_k$. This has been implemented in the MOOSE-2 solver [38]. Unfortunately, MOOSE-2 does not allow for averaging over prior and involves true unknown value of the parameter.

Numerous variants of the aforementioned methods can be found in the literature. For example, instead of the D-optimality criterion, one may also consider maximizing $\mathrm{tr}\left(\mathit{M}\right)$ or ${\lambda }_{min}\left(\mathit{M}\right)$. However, the vast majority of methods are based on the principles stated above (see, e.g., [4]), that is, maximization of some functions of the Fisher information matrix. Finally, we note that the above methods employ the classical optimality criterion, averaged only over the prior distribution. Consequently, they are not fully Bayesian and, following the terminology of [11], should rather be referred to as pseudo-Bayesian methods.

6. Examples of Input Signal Design

In the following, we present four examples of optimal input signal design using both Bayesian and classical methods. Examples 1-3 are classical in nature and concern time-invariant linear stochastic systems. Examples 1 and 2 are elementary, while Example 3, taken from [28], addresses the design of a control signal for a paradigmatic model of the atomic sensor. The sensor is modelled as a harmonic oscillator with the natural frequency being the parameter of interest. In examples 1-3, the Bayesian approach is compared with classical methods.

Example 4, adapted from [29], is more advanced and considers the design of the pump laser control signal in an optically pumped magnetometer. The magnetometer is modelled as a quasi-linear stochastic system, where the matrices $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{G}$ depend nonlinearly on the control signal u. For this system, classical methods cannot be applied. Therefore, estimation errors are compared with the Information-Theoretic Lower Bound (ITB) provided in Theorem 1 and with the errors obtained by using an appropriately selected harmonic input signal.

In all examples, the errors were computed using the Monte Carlo method. The parameter $\mathit{\theta }$ and the initial conditions ${\mathit{x}}_0$ were sampled from the prior distribution $p_0({\mathit{x}}_0,\mathit{\theta })$. Observations $y_0,\dots ,y_N$ corresponding to the sampled parameters and initial conditions were then generated, and the error of the MAP estimator was calculated. This error was subsequently averaged over many repetitions of the procedure.

6.1. Elementary Example

We begin with a very simple first-order system

$$\begin{array}{cc}\hfill & x_{k+1}={\theta }_1{x}_k+{\theta }_2{u}_k+g{w}_k,\hfill \end{array}$$

(91)

$$\begin{array}{cc}\hfill & y_k=x_k+{\sigma }_v{v}_k,k=0,1,...,N,\hfill \end{array}$$

(92)

with ${\sigma }_v=0.1$, $g=0.01$. The parameter vector $\mathit{\theta }={[{\theta }_1,{\theta }_2]}^{\mathsf{T}}$, has a prior distribution $p_0\left(\mathit{\theta }\right)=\mathcal{N}(\mathit{\theta },{\mathit{m}}_{\theta },{\mathbf{S}}_{\theta })$, where ${\mathit{m}}_{\theta }={\left[\begin{array}{cc}0.8& 0.2\end{array}\right]}^{\mathsf{T}}$, ${\mathbf{S}}_{\theta }={10}^{-2}\mathbb{I}$. As assumed in Section 4, the initial condition $x_0$ is conditionally Gaussian, that is, $p\left(x_0\right|\mathit{\theta })=\mathcal{N}(x_0,m_0^{-}\left(\mathit{\theta }\right),s_0^{-}\left(\mathit{\theta }\right))$, with $m_0^{-}\left(\mathit{\theta }\right)=0$, $s_0^{-}\left(\mathit{\theta }\right)=0.01$. The length of the signal $N=100$ and the set of admissible signals is given by (3) with $\tilde{\mathit{U}}=0$, that is, the norm of the signal cannot be greater than $\varrho$. To minimize the averaged D-optimal criterion (89), we need to calculate the sensitivity of the prediction error. The sensitivity equations (87) now take the form

$$\begin{array}{cc}\hfill & (1+(K\left({\theta }_1\right)-{\theta }_1)z^{-1})(1-{\theta }_1{z}^{-1}){\psi }_{1,k}={\theta }_2{z}^{-2}{u}_k,\hfill \end{array}$$

(93)

$$\begin{array}{cc}\hfill & (1+(K\left({\theta }_1\right)-{\theta }_1)z^{-1}){\psi }_{2,k}=z^{-1}{u}_k,\hfill \end{array}$$

(94)

where the Kalman gain $K\left({\theta }_1\right)$ is given by (82), (83) with $\mathbf{A}={\theta }_1$, $\mathbf{G}=g$, $\mathbf{C}=1$.

The optimal input signals were designed by maximizing the Bayesian criterion (44), the averaged D-optimal criterion (89), and the spectral criterion (90), subject to the constraint (3) with $\tilde{\mathit{U}}=0$.

Maximization of the spectral criterion (90) was performed using the MOOSE-2 solver [38], evaluated at $\mathit{\theta }={\mathit{m}}_{\theta }$ with default parameters, that is, the input spectrum was FIR-type with 20 lags and the spectrum power constraint was set to 1 (prob.spectrum.signal.power.ub=1). There were no additional constraints on the shape of the spectrum.

The optimal signals and the corresponding estimation errors of ${\theta }_1$ and ${\theta }_2$ are shown in Figure 1 and Figure 2. In Fig. Figure 2 we also calculate the estimation errors for the constant (step) signal, which is certainly not optimal. The constant signal and the MOOSE signal were always assigned a norm equal to $\varrho$.

6.2. Example with a Non-Gaussian Prior Distribution

Consider the following system:

$$d\mathit{x}=\left({\mathbf{A}}_c\left(\theta \right)+{\mathbf{B}}_c\left(\theta \right)u\right)dt+{\mathbf{G}}_c\left(\theta \right)d\mathit{w},$$

(95)

$$y_k=\mathbf{C}{\mathit{x}}_k+s_v{v}_k,$$

(96)

where

$${\mathbf{A}}_c\left(\theta \right)=\left[\begin{array}{cc}0& 1\\ 0& -\theta \end{array}\right],{\mathbf{B}}_c=\left[\begin{array}{c}0\\ \theta \end{array}\right],{\mathbf{G}}_c\left(\theta \right)=\left[\begin{array}{c}0\\ \sqrt{d_c\theta }\end{array}\right],\mathbf{C}=\left[\begin{array}{cc}1& 0\end{array}\right],$$

(97)

${\mathit{x}}_k=\mathit{x}\left(t_k\right)$, $t_k=k\Delta$, $\Delta =0.05·{10}^{-3}$, $d_c=0.01$, $s_v=0.1$. This system can be considered as controlled Brownian motion or DC motor with stochastic disturbances. The parameter $\theta$ is the unknown damping rate of the system. Assuming that $u\left(t\right)=u_k,t\in [t_k,t_{k+1}]$, the discrete-time system corresponding to (95) has the form

$${\mathit{x}}_{k+1}=\mathbf{A}\left(\theta \right){\mathit{x}}_k+\mathbf{B}\left(\theta \right)u_k+\mathbf{G}\left(\theta \right)w_k,$$

(98)

where, according to the procedure given in Appendix C:

$$\begin{array}{cc}\hfill & \mathbf{A}\left(\theta \right)=\left[\begin{array}{cc}1& {\displaystyle \frac{1-e^{-\theta \Delta }}{\theta }}\\ 0& e^{-\theta \Delta }\end{array}\right],\mathbf{B}\left(\theta \right)=\left[\begin{array}{c}\Delta -{\displaystyle \frac{1-e^{-\theta \Delta }}{\theta }}\\ 1-e^{-\theta \Delta }\end{array}\right],\hfill \\ \hfill & \mathbf{G}\left(\theta \right)=\sqrt{d_c\theta }\left[\begin{array}{cc}\sqrt{D_{1,1}\left(\theta \right)}& 0\\ {\displaystyle \frac{D_{1,2}\left(\theta \right)}{\sqrt{D_{1,1}\left(\theta \right)}}}& \sqrt{{\displaystyle \frac{D_{1,1}\left(\theta \right)D_{2,2}\left(\theta \right)-D_{1,2}{\left(\theta \right)}^2}{D_{1,1}\left(\theta \right)}}}\end{array}\right],\hfill \end{array}$$

(99)

where

$$\begin{array}{cc}\hfill & D_{1,1}\left(\theta \right)={\displaystyle \frac{4e^{-\theta \Delta }-e^{-2\theta \Delta }+2\theta \Delta -3}{2{\theta }^3}},\hfill \end{array}$$

(100)

$$\begin{array}{cc}\hfill & D_{1,2}\left(\theta \right)={\displaystyle \frac{1-2e^{-\theta \Delta }+e^{-2\theta \Delta }}{2{\theta }^2}},D_{2,2}\left(\theta \right)={\displaystyle \frac{1-e^{-2\theta \Delta }}{2\theta }}.\hfill \end{array}$$

(101)

The initial condition is Gaussian with ${\mathit{m}}_0^{-}=0$, ${\mathbf{S}}_0^{-}=\mathrm{diag}[0.001,0.005]$. Unlike in the previous example, here we assume that the prior distribution of $\theta$ is uniform, that is, $p_0\left(\theta \right)=U[a,b]$ with $a=0.05$, $b=2$. Following the Gauss-Legendre formula (43), we get ${\theta }_1=0.5(a+b-(b-a)/\sqrt{3})\approx 0.462$, ${\theta }_2=0.5(a+b+(b-a)/\sqrt{3})\approx 1.588$, $p_{0,1}=p_{0,2}=0.5$. Thus, $r=2$ in (31) and, accordingly to (35), the Bayesian optimal signal is a solution of the simplified and convex optimization problem (36) with $d_{1,2}$ defined by Lemmas 3 and 4. Moreover, since the matrices in (99), do not depend on $u_k$, the last two terms in (77) can be omitted. The set of admissible signals is given by (3) with $\tilde{\mathit{U}}=0$, that is, the signal norm cannot be greater than $\varrho$.

In order to employ the classical methods described in Sect. Section 5, it is necessary to first evaluate the sensitivity of the prediction error. The transfer functions G and H in (80) have the form

$$G(\theta ,z)={\displaystyle \frac{B(\theta ,z)}{A(\theta ,z)}}{z}^{-1},H(\theta ,z)={\displaystyle \frac{C(\theta ,z)}{A(\theta ,z)}},$$

(102)

where

$$\begin{array}{cc}\hfill & A(\theta ,z)=1-(1+e^{\theta \Delta }))z^{-1}+e^{-\theta \Delta }{z}^{-2},\hfill \end{array}$$

(103)

$$\begin{array}{cc}\hfill & B(\theta ,z)=\Delta -{\displaystyle \frac{1-e^{-\theta \Delta }}{\theta }}+\left({\displaystyle \frac{1}{\theta }}-\left({\displaystyle \frac{1}{\theta }}+\Delta \right)e^{-\theta \Delta }\right)z^{-1},\hfill \end{array}$$

(104)

$$\begin{array}{cc}\hfill & C(\theta ,z)=1+\left({\mathbf{K}}_1\left(\theta \right)-\left(1-e^{-\theta \Delta }\right)\right)z^{-1}+\left({\displaystyle \frac{1-e^{-\theta \Delta }}{\theta }}{\mathbf{K}}_2\left(\theta \right)+(1-{\mathbf{K}}_1\left(\theta \right))e^{-\theta \Delta }\right)z^{-2},\hfill \end{array}$$

(105)

and the Kalman gain $\mathbf{K}$ is given by (82). Since we only have one parameter, the sensitivity ${\psi }_k$ is a number, and the sensitivity equation (87) now takes the form

$$A(\theta ,z)C(\theta ,z){\psi }_k(\theta ,\mathit{U})=\left({\displaystyle \frac{\partial B(\theta ,z)}{\partial \theta }}A(\theta ,z)-B(\theta ,z){\displaystyle \frac{\partial A(\theta ,z)}{\partial \theta }}\right)z^{-1}{u}_k.$$

(106)

The D-optimal signal is then obtained by maximization of the averaged D-optimal criterion

$$Q\left(\mathit{U}\right)={\mathbb{E}}_{p_0\left(\theta \right)}\left[{\displaystyle \frac{1}{N}}\sum _{k=1}^N{\psi }_k^2(\theta ,\mathit{U})\right],$$

(107)

with constraints (3).

The optimal input signals were designed by maximizing the Bayesian criterion (44) and the averaged D-optimal criterion (107), subject to the constraint (3) with $\tilde{\mathit{U}}=0$. The results are presented in Figure 3 and Figure 4. Figure 4 also shows the estimation error for the step signal (constant) and the PRBS signal. The constant and PRBS signals were always assigned a norm equal to $\varrho$.

6.3. Optimal Input Design for Atomic Sensor Model

In [28], a simplified paradigmatic model of an atomic sensor (an atomic magnetometer, [39,40]) is introduced, in which the dynamics is governed by oscillations of the collective spin of an atomic ensemble subjected to an external magnetic field. The system is driven by circularly polarized light from a pump laser, whose frequency acts as the input signal. A linearly polarized probe laser illuminates the atoms, and upon transmission through the medium, its polarization undergoes a Faraday rotation. The $J_z$ component of the collective spin is inferred from the measurement of the probe laser’s polarization angle. The model presented in [28] describes the dynamics of the spin components $\mathit{J}={[J_y,J_z]}^{\mathsf{T}}$ and has the form

$$d\mathbf{J}=\left[\begin{array}{cc}-{\displaystyle \frac{1}{T_2}}& {\omega }_L\\ -{\omega }_L& -{\displaystyle \frac{1}{T_2}}\end{array}\right]Jdt+\left[\begin{array}{c}0\\ 1\end{array}\right]\mathcal{E}\left(t\right)dt+d{\mathit{w}}^{\left(J\right)},$$

(108)

where ${\omega }_L$ is the Larmor frequency, $T_2=0.87$ ms, is the relaxation time, $\mathcal{E}$ is the pumping laser frequency, and ${\mathit{w}}^{\left(J\right)}$ is the Wiener process with known covariance $q\mathbb{I}$. The observation has the form $I_k=g_D{J}_z(k\Delta )+{\xi }_k$, $k=0,1,...$, where $I_k$ is the photocurrent, $\Delta =5$$\mu \mathrm{s}$, is the sampling time, ${\xi }_k\sim \mathcal{N}(0,{\sigma }_{\xi }^2)$ and $g_D$, ${\sigma }_{\xi }$ are known parameters. The Larmor frequency and the external magnetic field B are related to each other by the formula ${\omega }_L={\gamma }_{e}B$, where ${\gamma }_e$ is the gyromagnetic ratio. Hence, by measuring ${\omega }_L$, one can determine the field B. Taking $T_2$ as the time unit and after rescaling the time, state variables, observations, and the input signal $\mathcal{E}$, we get the following, equivalent to (108), stochastic system:

$$d\mathit{x}=\left({\mathbf{A}}_c\left(\theta \right)+{\mathbf{B}}_{c}u\right)dt+{\mathbf{G}}_{c}d\mathit{w},$$

(109)

$$y_k=\mathbf{C}{\mathit{x}}_k+s_v{v}_k,$$

(110)

where

$${\mathbf{A}}_c\left(\theta \right)=\left[\begin{array}{cc}-1& \theta \\ -\theta & -1\end{array}\right],{\mathbf{B}}_c=\left[\begin{array}{c}0\\ b_c\end{array}\right],{\mathbf{G}}_c=\sqrt{2}\mathbb{I},\mathbf{C}=\left[\begin{array}{cc}0& 1\end{array}\right],$$

(111)

${\mathit{x}}_k=\mathit{x}\left(t_k\right)$, $t_k=k\Delta$, $\Delta =5.7471·{10}^{-3}$, $s_v=11.85$, $b_c={10}^5$. The input signal u in (109) corresponds to $\mathcal{E}$ in (108). The parameter $\theta$ in (109) is related to the Larmor frequency ${\omega }_L$ in (108), by formula $\theta ={\omega }_L{T}_2$. Since the estimation error of the parameter $\theta$ depends on the input signal u, a natural question arises as to what form this signal should take. To solve this problem, we will go to discrete time and apply the methodology described in Section 4 and Section 5. Assuming that $u\left(t\right)=u_k,t\in [t_k,t_{k+1}]$, the discrete-time system corresponding to (109) has the form

$${\mathit{x}}_{k+1}=\mathbf{A}\left(\theta \right){\mathit{x}}_k+\mathbf{B}\left(\theta \right)u_k+\mathbf{G}{w}_k,$$

(112)

where, according to the procedure given in Appendix C:

$$\begin{array}{cc}\hfill & \mathbf{A}\left(\theta \right)=e^{-\Delta }\left[\begin{array}{cc}cos\theta \Delta & sin\theta \Delta \\ -sin\theta \Delta & cos\theta \Delta \end{array}\right],\mathbf{B}\left(\theta \right)={\displaystyle \frac{b_c}{1+{\theta }^2}}\left[\begin{array}{c}\theta -e^{-\Delta }(\theta cos\theta \Delta +sin\theta \Delta )\\ 1-e^{-\Delta }(cos\theta \Delta -\theta sin\theta \Delta )\end{array}\right],\hfill \end{array}$$

(113)

$$\begin{array}{cc}\hfill & \mathbf{G}=\sqrt{1-e^{-2\Delta }}\mathbb{I}.\hfill \end{array}$$

(114)

We assume that the prior distribution of $\theta$ is Gaussian, that is, $p_0\left(\theta \right)=\mathcal{N}(\theta ,m_{\theta },s_{\theta })$ with $m_{\theta }=54.6637$, $s_{\theta }=10.76$, which corresponds to the Larmor frequency of 10 kHz and its initial uncertainty of the order of 600 Hz ($3\sigma$). At the beginning of the process, the system is in thermal equilibrium, that is, $p\left({\mathit{x}}_0\right|\mathit{\theta })=\mathcal{N}({\mathit{x}}_0,0,\mathbb{I})$. Following Lemma 2, we get ${\theta }_1=m_{\theta }-\sqrt{s_{\theta }}\approx 51$, ${\theta }_2=m_{\theta }+\sqrt{s_{\theta }}\approx 58$, $p_{0,1}=p_{0,2}=0.5$. Thus, $r=2$ in (31) and, accordingly to (35), the Bayesian optimal signal is a solution of the simplified and convex optimization problem (36) with $d_{1,2}$ defined by Lemmas 3 and 4. Moreover, since the matrices in (113), (114) do not depend on $u_k$, the last two terms in (77) can be omitted. The set of admissible signals is given by (3) with $\tilde{\mathit{U}}=0$, that is, the signal norm cannot be greater than $\varrho$.

In order to employ the classical methods described in Sect. Section 5, it is necessary to first evaluate the sensitivity of the prediction error. Similarly as in the previous example, the polynomials $A,B,C$ have the form

$$\begin{array}{cc}\hfill & A(\theta ,z)=1-2e^{-\Delta }cos(\theta \Delta )z^{-1}+e^{-2\Delta }{z}^{-2},\hfill \end{array}$$

(115)

$$\begin{array}{cc}\hfill & B(\theta ,z)={\mathbf{B}}_2\left(\theta \right)-e^{-\Delta }\left({\mathbf{B}}_1\left(\theta \right)sin(\theta \Delta )+{\mathbf{B}}_2\left(\theta \right)cos(\theta \Delta )\right)z^{-1},\hfill \end{array}$$

(116)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & C(\theta ,z)=1+\left({\mathbf{K}}_2\left(\theta \right)-2e^{-\Delta }cos(\theta \Delta )\right)z^{-1}+\hfill \\ \hfill & +e^{-\Delta }\left(e^{-\Delta }-{\mathbf{K}}_1\left(\theta \right)sin(\theta \Delta )-{\mathbf{K}}_2\left(\theta \right)cos(\theta \Delta )\right)z^{-2},\hfill \end{array}\end{array}$$

(117)

and the Kalman gain $\mathbf{K}$ and the vector $\mathbf{B}$ are given by (82) and (113), respectively. The sensitivity equation (87) is given by (106). The D-optimal signal is then obtained by maximization of the averaged D-optimal criterion (107) with constraints (3).

The optimal input signals were designed by maximizing the Bayesian criterion (44), the averaged D-optimal criterion (107), and the spectral criterion (90), subject to the constraint (3) with $\tilde{\mathit{U}}=0$. The results are presented in Figure 5 and Figure 6. Figure 6 also shows the estimation error for the step (constant) signal and harmonic signal $u\left(t\right)=acos\left(m_{\theta }t\right)$. Frequency of the harmonic signal was equal to the expected value of the a priori distribution of the parameter $\theta$. The constant, MOOSE and harmonic signals were always assigned a norm equal to $\varrho$.

Maximization of the spectral criterion (90) was performed using the MOOSE-2 solver [38], evaluated at $\mathit{\theta }={\mathit{m}}_{\theta }$ with default parameters, that is, the input spectrum was FIR-type with 20 lags and the spectrum power constraint was set to 1 (prob.spectrum.signal.power.ub=1). There were no additional constraints on the shape of the spectrum.

6.4. Bayesian Input Signal Design for Pump Laser in Optically Pumped Magnetometer

Optically pumped magnetometers operate by aligning atomic spins with a circularly polarized pump laser, after which the spins precess around the external magnetic field at the Larmor frequency. The probe laser measures this precession via polarization rotation (Faraday effect), linking the detected signal to the magnetic field [39,40]. The pump laser’s frequency strongly affect spin polarization and coherence time making precise laser control central to minimizing estimation error. The advanced control strategies can then suppress the noise and enhance sensitivity. Consequently, accurate control of the pumping laser is a key factor in achieving high-resolution and low-error magnetometer. We consider here the magnetometer model given by equation (S9) in the article [29]:

$${\displaystyle \frac{d\mathit{F}}{d\tau }}=\left(-{\gamma }_e\mathbf{B}+G{S}_3\widehat{z}\right)\times \mathit{F}-\gamma \mathit{F}+P\left(\tau \right)(\widehat{z}{F}_{max}-\mathit{F})+{\mathbf{G}}_0\left(P\left(\tau \right)\right)\mathit{w},$$

(118)

where $\mathit{F}={(F_x,F_y,F_z)}^{\mathsf{T}}$ is the collective atomic spin, ${\gamma }_e$ is electron gyromagnetic ratio, $\mathbf{B}={(B_x,B_y,B_z)}^{\mathsf{T}}$ is constant magnetic field vector, G is known positive constant, $G{S}_3\widehat{z}$ is the effective field produced by ac-Stark shifts due to the probe laser where $S_3$ is white Gaussian noise with variance ${\sigma }_3^2$. The optical pumping rate $P\left(\tau \right)⩾0$ is an input signal. The atomic spin noise ${\mathbf{G}}_0\left(P\left(\tau \right)\right)\mathit{w}$ is modelled as white Gaussian where $\mathit{w}={(w_x,w_y,w_z)}^{\mathsf{T}}$ is a vector of standard and mutually independent Wiener increments. The ${\mathbf{G}}_0$ matrix is diagonal and is given by

$${\mathbf{G}}_0\left(P\left(\tau \right)\right)=\sqrt{\frac{2}{3}F(F+1)N_A(\gamma +P\left(\tau \right))}\mathbb{I},$$

(119)

where $N_A$ is the number of atoms and F is known atomic spin number. Parameter $F_{max}=N_{A}F$ is the maximum possible polarization. The transverse relaxation rate $\gamma$ depends on the number of atoms and is given by $\gamma \left(N_A\right)=T_2{\left(N_A\right)}^{-1}={\gamma }_0+{10}^{-12}\alpha N_A$, where ${\gamma }_0$ and $\alpha$ are known positive constants and $T_2$ is the effective coherence time. The observation equation has the form

$$S_2=F_z+N_{S_2},$$

(120)

where $N_{S_2}$ denote measurement noise with variance ${\sigma }_2^2$. In the experiment, the $S_2$ component of the Stokes vector is measured at discrete time moments $t_k=k\Delta$, where $\Delta$ is sampling period. The realistic parameters of the model are given in Table 1.

In what follows, equation (118) will be interpreted in the Itö sense. Moreover, we assume that the noise $G{S}_3$ in (118) is small and can be omitted.

By introducing state variables $\mathit{\xi }=\mathit{F}\sqrt{{\displaystyle \frac{3}{F(F+1)N_A}}}$, control variable $u=P/\gamma$ and non-dimensional time $t=\gamma \tau$, and after multiplying both sides of (120) by $\sqrt{{\displaystyle \frac{3}{F(F+1)N_A}}}$, we get the following model:

$$d\mathit{\xi }=({\mathbf{A}}_c(\mathit{\theta },u)\mathit{\xi }+{\mathbf{B}}_{c}u)dt+{\mathbf{G}}_c\left(u\right)d\mathit{\eta },y_k={\mathit{\xi }}_3\left(t_k\right)+{\sigma }_v{v}_k,$$

(121)

where $\mathit{\eta }$ is the three-dimensional standard Wiener process with unit covariance and

$${\mathbf{A}}_c(\mathit{\theta },u)=\left[\begin{array}{ccc}-(1+u)& {\theta }_3& -{\theta }_2\\ -{\theta }_3& -(1+u)& {\theta }_1\\ {\theta }_2& -{\theta }_1& -(1+u)\end{array}\right],{\mathbf{B}}_c=\left[\begin{array}{c}0\\ 0\\ b_c\end{array}\right],{\mathbf{G}}_c\left(u\right)=\sqrt{2(1+u)}\mathbb{I},$$

(122)

with $b_c=\sqrt{\frac{3N_{A}F}{F+1}}$, ${\sigma }_v={\sigma }_2\sqrt{{\displaystyle \frac{3}{F(F+1)N_A}}}$. Taking the parameters from Table 1 we have $b_c=1.22·{10}^6$, ${\sigma }_v=11.85$. The parameter vector $\mathit{\theta }={({\theta }_1,{\theta }_2,{\theta }_3)}^{\mathsf{T}}$, represents the external magnetic field due to the relation $\mathit{\theta }={\gamma }_e{T}_2\mathbf{B}$. If u is a constant signal, then system (121) approaches thermodynamic equilibrium, with $\mathbb{E}\mathit{x}\left(t\right)=-{\mathbf{A}}_c{(\mathit{\theta },u)}^{-1}{\mathbf{B}}_{c}u$ and $\mathrm{cov}\left(\mathit{x}\right(t\left)\right)=\mathbb{I}$.

A closer examination of equation (121) shows that the component ${\mathit{\xi }}_3(t,\mathit{\theta })$ of its solution remains invariant under rotations of the vector $\mathit{\theta }$ about the z-axis. As a result, $\mathit{\theta }$, and hence the field $\mathbf{B}$, cannot be uniquely identified from the observations $y_0,\dots ,y_N$. The only quantities that can be uniquely identified in this model are the magnitude of the vector $\mathbf{B}$ and the angle $\eta$ between $\mathbf{B}$ and one of the coordinate axes, say the $\widehat{z}$ axis. However, to simplify the problem as much as possible, we introduce here the additional assumption that the field $\mathbf{B}$ always lies in the x–y plane, that is, $\mathbf{B}={(B_x,B_y,0)}^{\mathsf{T}}$. With this assumption, the change of variables

$$x_1={\xi }_{1}sin\phi -{\xi }_{2}cos\phi ,x_2={\xi }_3,$$

(123)

$$cos\phi ={\displaystyle \frac{{\theta }_1}{\sqrt{{\theta }_1^2+{\theta }_2^2}}},sin\phi ={\displaystyle \frac{{\theta }_2}{\sqrt{{\theta }_1^2+{\theta }_2^2}}},$$

(124)

reduces model (121) to a two-dimensional system:

$$d\mathit{x}=({\mathbf{A}}_c(\theta ,u)x+{\mathbf{B}}_{c}u)dt+{\mathbf{G}}_c\left(u\right)d\mathit{w},y_k=\mathbf{C}{\mathit{x}}_k+{\sigma }_v{v}_k,$$

(125)

where $\mathit{x}={(x_1,x_2)}^{\mathsf{T}}$, $\theta ={\gamma }_e{T}_2\sqrt{B_x^2+B_y^2}$, $\mathit{w}$ is a two-dimensional standard Wiener process with unit covariance and

$${\mathbf{A}}_c(\theta ,u)=\left[\begin{array}{cc}-(1+u)& -\theta \\ \theta & -(1+u)\end{array}\right],{\mathbf{B}}_c=\left[\begin{array}{c}0\\ b_c\end{array}\right],\mathbf{C}=\left[\begin{array}{cc}0& 1\end{array}\right],{\mathbf{G}}_c\left(u\right)=\sqrt{2(1+u)}\mathbb{I}.$$

(126)

Hence, under the assumption $B_z=0$, the observations $y_0,\dots ,y_N$, variable ${\mathit{\xi }}_3$, and the ${\mathit{F}}_z$ component of the collective spin, are fully characterized by the reduced model (125). Furthermore, within this reduced model it can be readily verified that $\theta$ is uniquely identifiable. Naturally, the accuracy of estimating $\theta$ depends on the choice of input u. To determine an input u that maximizes the information about $\theta$, we now turn to the discrete-time formulation of (125) and apply the methods described in Sect. Section 3 and Section 4. Assuming the control signal is piecewise constant, that is, $u\left(t\right)=u_k,t\in [t_k,t_{k+1}]$, the process ${\mathit{x}}_k=\mathit{x}\left(t_k\right)$ satisfies the difference equation

$${\mathit{x}}_{k+1}=\mathbf{A}(\theta ,u_k){\mathit{x}}_k+\mathbf{B}(\theta ,u_k)+\mathbf{G}\left(u_k\right){\mathit{w}}_k,$$

(127)

where ${\mathit{w}}_k\sim \mathcal{N}(0,\mathbb{I})$ and the matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{G}$ can be calculated following the procedure given in Appendix C. Upon completion of straightforward calculations, we get:

$$\begin{array}{cc}\hfill & \mathbf{A}(\theta ,u_k)=e^{-(1+u_k)\Delta }\left[\begin{array}{cc}cos\theta \Delta & -sin\theta \Delta \\ sin\theta \Delta & cos\theta \Delta \end{array}\right],\mathbf{G}\left(u_k\right)=\sqrt{1-e^{-2(1+u_k)\Delta }}\mathbb{I},\hfill \end{array}$$

(128)

$$\begin{array}{cc}\hfill & \mathbf{B}(\theta ,u_k)={\displaystyle \frac{b_c{u}_k}{(1+u_k)+{\theta }^2}}\left[\begin{array}{c}{e}^{-(1+u_k)\Delta }\left(\theta cos(\theta \Delta )+(1+u_k)sin(\theta \Delta )\right)-\theta \\ e^{-(1+u_k)\Delta }\left(\theta sin\left(\theta \right)-(1+u_k)cos(\theta \Delta )\right)+(1+u_k)\end{array}\right].\hfill \end{array}$$

(129)

At the beginning of the process, the system is in a thermal equilibrium corresponding to $u\equiv 0$. Hence, $p\left({\mathit{x}}_0\right|\mathit{\theta })=\mathcal{N}({\mathit{x}}_0,0,\mathbb{I})$. We also assume that the prior distribution of $\theta$ is Gaussian, that is, $p_0\left(\theta \right)=\mathcal{N}(\theta ,m_{\theta },s_{\theta })$ with $m_{\theta }=54.6637$, $s_{\theta }=3·{10}^{-3}$, which corresponds to the Larmor frequency of 10 kHz and its initial uncertainty of the order of 30 Hz ($3\sigma$). Similarly as in Sect. Section 6.3, we get from Lemma 2: ${\theta }_1=m_{\theta }-\sqrt{s_{\theta }}\approx 54.61$, ${\theta }_2=m_{\theta }+\sqrt{s_{\theta }}\approx 54.72$, $p_{0,1}=p_{0,2}=0.5$. Since $r=2$ in (31), then accordingly to (35), the Bayesian optimal signal is a solution of the simplified optimization problem (36) with $d_{1,2}$ calculated by Lemmas 3 and 4. Unlike in the previous examples, in this problem we maximize criterion (36) with constraints on the signal amplitude, that is, $0⩽u_k⩽u_{max}$, which is preferable in realistic scenarios.

The results are presented in Figure 7, Figure 8, and Figure 9. The optimal input signal consistently lies on the boundary of the admissible set. For small $u_{max}$, the optimal signal is rectangular, with a frequency close to the a priori Larmor frequency. Since large values of $u\left(t\right)$ strongly damp spin oscillations and increases the noise, the optimal signals for large $u_{max}$ consist of short pulses at the maximum admissible amplitude. Once the oscillations decay, the system should be re-excited by a new sequence of short pulses, repeated periodically, as illustrated in Figure 9. The harmonic signal $u\left(t\right)=0.5u_{max}(1+cos\left(m_{\theta }t\right))$ is nearly optimal for small $u_{max}$ but becomes ineffective for large $u_{max}$, as it strongly damps the oscillations (see the lower right panel of Figure 8). As a result, the measurements carry less information about the Larmor frequency, and the estimation error increases despite the higher signal amplitude. An analogous behaviour is observed for rectangular signals. More generally, let $s\left(t\right)\in [0,1]$, be any signal and define $u\left(t\right)=as\left(t\right)$ with $a⩾0$. Then, as illustrated in Figure 7, estimation error reaches a minimum for some non-zero value of the parameter a.

Extending the experimental duration from 2 to 5 ms reduces the estimation error by a factor of 2 compared to the case shown in Figure 7. For $u_{max}=200$ and an experiment duration of 5 ms, the harmonic input signal yields an estimation error of 7 mHz, while the optimal signal, shown in the lower left panel of Figure 9, reduces the error to $0.48$ mHz, that is, approximately 14 times smaller. Finally, the estimation error attains the Information-Theoretic Lower Bound (20), demonstrating that in this case the MAP estimator (10) achieves optimal performance.

It should be noted that the above models assume the Markovian environment and this condition should be checked in an experiment. To do this, one can use the criterion given in [41]. Non-Markovian models are much more complicated (see e.g. [42]) and one would need to employ a noise model with long memory. To model long-memory noise, fractional-order stochastic equations can be used instead of (118). Such models capture long-memory effects, and their noise correlation function decays slowly, for example as $t^{-1/2}$.

To implement the proposed method in real time, one can proceed as follows. First, observe that the pump signal has a simple structure, consisting of short pulses at the maximum admissible amplitude, each lasting approximately $5\mu \mathrm{s}$. These pulses should be repeated with a period of roughly $2T_2$, and each pulse should be triggered when the vector $[F_y,F_z]$ forms an angle of about ${30}^{\circ }$ with the z-axis (i.e., ${30}^{\circ }$ before the maximum of $F_z$). To estimate the unknown vector $\mathit{F}$ and the Larmor frequency, the MAP estimator is too slow for real-time applications. Instead, an Extended Kalman Filter (EKF) can be employed in a manner roughly similar to that described in [43] and [44]. This approach is considered feasible for implementation in an experimental setup.

7. Discussion and Conclusions

This paper has developed a Bayesian framework for optimal input signal design in the identification of quasi-linear stochastic dynamical systems. By an Information-Theoretic Lower Bound on estimation error and its connection to the Bayesian Cramér–Rao Bound, we showed that maximizing mutual information provides a principled alternative to Fisher-information–based criteria. The proposed method relies on the maximization of the MI lower bound (30), which produces a tractable surrogate objective for both finite parameter sets and parameter spaces of continuum cardinality. A key contribution is the algorithmic reduction of the dimension of covariance matrices required for inversion by a factor of N, making the method feasible for long-term experiments.

The comparison with the average D-optimal design highlights the practical benefits of the Bayesian approach. While classical methods are computationally efficient, they require complex differentiations to evaluate parameter sensitivities and may yield suboptimal results when parameter uncertainty is large or when the system exhibits significant non-linearities. In contrast, the proposed Bayesian method requires only the system matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, $\mathbf{G}$, together with the prior distributions of the parameter and initial conditions, without the need to calculate derivatives of prediction errors. This makes the method applicable to a much broader class of systems, while also enabling it to handle large initial parameter uncertainty.

The method also has certain limitations. The lower bound of the MI involves exponential terms that can vanish when the pairwise distance factors $d_{i,j}\left(\mathit{U}\right)$ are large, which can cause numerical problems. However, this drawback can be mitigated by appropriate scaling of the optimization problem. If we consider the simplified optimization problem (36), with only two candidate parameter values, these numerical problems never occur. The discrete approximation of the MI (29) is a potential source of problems and the weights and nodes in (38) should be carefully selected to achieve sufficient approximation accuracy. The third limitation arises from the fact that the maximized criterion is only a lower bound on the MI and is generally not tight. Consequently, there certainly exists a class of problems for which maximizing this lower bound is inefficient and may generate signals that are far from optimality in the sense of maximizing the MI (19)

In all analyzed examples, the proposed Bayesian approach, although approximate, generated signals no worse, or even better (see Figure 2 and Figure 4), than the classical methods.

The second example illustrates that a non-Gaussian prior distribution leads to increased errors in the average D-optimal method. For a Gaussian prior distribution, it was confirmed that both the average D-optimal and the proposed Bayesian method yield identical results. This observation underscores the sensitivity of the classical approach to the form of the prior distribution and highlights the necessity of employing estimation techniques that are robust to non-Gaussian priors. In the third example, the D-optimal method produces results almost identical to those of the Bayesian approach. To explain this, note that in this problem the prior distribution of the parameter $\theta$ is relatively narrow. Then the function $d_{1,2}\left(\mathit{U}\right)$, which we minimize in this task, is approximately proportional to the sensitivity of the output to changes in $\theta$. Thus, $d_{1,2}\left(\mathit{U}\right)$ can be interpreted as a quantity proportional to the Fisher information. Consequently, the resulting input signals and the corresponding estimation errors are nearly identical.

The study of atomic sensor models further demonstrates the practical relevance of the approach. The optimal signals in these examples are always better than the harmonic signal with a frequency equal to the expected natural frequency of the oscillator. The fourth example, a seemingly minor modification of the oscillator from the third example, shows that the dependence of the system matrices on the control signal is significant and leads to completely different optimal signals. In the analyzed examples, the MAP estimator achieves an Information-Theoretic Lower Bound (20), but this is not always the case and, depending on the task, there are better estimators. Unfortunately, finding them is difficult.

Since the method produces the posterior distribution $p(\mathit{\theta },{\mathit{x}}_k|{\mathit{Y}}_k)$, it can be easily converted to a sequential Bayesian Adaptive Design (BAD) algorithm [5,6]. Then the optimal strategy is a functional of the posterior, that is, ${\mathit{u}}_k={\mathit{\varphi }}_k(p\left(\mathit{\theta }\right|{\mathit{Y}}_k),{\mathit{m}}_k\left(\mathit{\theta }\right),{\mathbf{S}}_k\left(\mathit{\theta }\right))$. In the simplest case, when the matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{G}$ do not depend on ${\mathit{u}}_k$ and $\Theta =\{{\mathit{\theta }}_{\mathbf{1}},{\mathit{\theta }}_2\}$, the optimal strategy ${\mathit{\varphi }}_k$ can be determined by maximizing (30) on the trajectories of the system (71)-(75). This problem is deterministic, therefore, relatively simple, and can be solved using deterministic optimal control methods.

From a broader perspective, quasi-linear systems arise naturally in quantum mechanics, chemical engineering, and thermal processes, making the proposed method widely applicable. In conclusion, this work provides both theoretical justification and practical tools for Bayesian input design in quasi-linear stochastic systems. By bridging information-theoretic principles with efficient computational methods, it establishes a foundation for robust experimental design in a wide range of applications. The results reported here should stimulate further research at the intersection of Bayesian inference, control, and the identification of non-linear systems.