Exploiting Convexification for Bayesian Optimal Sensor Placement by Maximization of Mutual Information

TL;DR

This paper introduces a convex relaxation approach to efficiently determine optimal sensor placements for Bayesian inference in structural systems, significantly reducing computational costs compared to exhaustive search methods.

Contribution

A novel convex relaxation scheme is proposed for Bayesian sensor placement, enabling efficient and often optimal solutions without exhaustive combinatorial search.

Findings

The method finds solutions often binary and optimal for sensor placement.

It handles large configuration spaces, exemplified by over 47 trillion options.

Compared to greedy and heuristic methods, it offers improved efficiency and accuracy.

Abstract

Bayesian optimal sensor placement, in its full generality, seeks to maximize the mutual information between uncertain model parameters and the predicted data to be collected from the sensors for the purpose of performing Bayesian inference. Equivalently, the expected information entropy of the posterior of the model parameters is minimized over all possible sensor configurations for a given sensor budget. In the context of structural dynamical systems, this minimization is computationally expensive because of the large number of possible sensor configurations. Here, a very efficient convex relaxation scheme is presented to determine informative and possibly-optimal solutions to the problem, thereby bypassing the necessity for an exhaustive, and often infeasible, combinatorial search. The key idea is to relax the binary sensor location vector so that its components corresponding to all possible sensor locations lie in the unit interval. Then, the optimization over this vector is a convex problem that can be efficiently solved. This method always yields a unique solution for the relaxed problem, which is often binary and therefore the optimal solution to the original problem. When not binary, the relaxed solution is often suggestive of what the optimal solution for the original problem is. An illustrative example using a fifty-story shear building model subject to sinusoidal ground motion is presented, including a case where there are over 47 trillion possible sensor configurations. The solutions and computational effort are compared to greedy and heuristic methods.

Figures

Tables

Table 1. Table 1: Simulation results

Case				Sensor DOF # (1 $\equiv$ Base, $N_d\equiv$ Roof)
$N_d$	$N_o$	$N$	$N_k$
2	1	1000	1000	2
2	1	1000	2000	2
2	1	2000	1000	2
2	1	2000	2000	2
4	1	1000	1000	4
4	1	1000	2000	4
4	1	2000	1000	4
4	1	2000	2000	4
4	2	1000	1000	2, 4
4	2	1000	2000	2, 4
4	2	2000	1000	2, 4
4	2	2000	2000	2, 4
8	2	1000	1000	6, 8
8	2	2000	2000	6, 8

Table 2. Table 2: Description of various sensor configurations for comparative example

Case	Description
$z^{\ast }$	Optimal configuration using convex relaxation
$z_{\text{low}}$	Stories 1 through 20 instrumented
$z_{\text{high}}$	Stories 31 through 50 instrumented
$z_{\text{common}}$	Sensors evenly spaced (Stories 1 through 50 in steps of 2.5, rounded up)
$z_{\text{greedy}}$	Solution using greedy sequential placement

Table 3. Table 3: Comparison of information gain relative to optimal configuration

Case	Objective Value	# bits gain using z*
$z^{\ast }$	6.14E+01
$z_{\text{low}}$	5.75E+01	5.7
$z_{\text{high}}$	6.12E+01	0.3
$z_{\text{common}}$	6.05E+01	1.3
$z_{\text{greedy}}$	6.14E+01	0.0

Equations

$$y_n=S_o(\delta )\left(x_n({\theta }_s)+{ϵ}_n({\theta }_e)\right)$$

$$p(y_{1:N}|\theta ,\mathcal{M})=\prod _{n=1}^{N}p\left(y_n|\theta ,\mathcal{M}\right)$$

$$p(\theta |y_{1:N},\mathcal{M})=\frac{p(y_{1:N}|\theta ,\mathcal{M})p(\theta |\mathcal{M})}{p(y_{1:N}|\mathcal{M})}$$

$$p(y_{1:N}|\mathcal{M})=\int p(y_{1:N}|\theta ,\mathcal{M})p(\theta |\mathcal{M})d\theta$$

$$\mathbb{E}\left[{ϵ}_n{ϵ}_m^T\right]={\sigma }^2{I}_{N_d}{\delta }_{mn}\text{ where }{ϵ}_n\sim \mathcal{N}(\cdot |0,{\sigma }^2{I}_{N_d})\text{ i.i.d.}$$

$$p(y_{1:N}|\theta ,\mathcal{M})=\prod _{n=1}^N\mathcal{N}\left(y_n|S_o(\delta )x_n(\theta ),{\sigma }^2\right)$$

$$\log p(y_{1:N}|\theta ,\mathcal{M})=-\frac{N{N}_o}{2}\log 2\pi {\sigma }^2-\frac{1}{2{\sigma }^2}\sum _{n=1}^N\Vert y_n-S_o(\delta )x_n({\theta }_s){\Vert }^2$$

$$I(y_{1:N}(\delta ),\theta |\mathcal{M})$$

$$=\int \int p(y_{1:N}(\delta ),\theta |\mathcal{M})\left[\log \frac{p(\theta |y_{1:N}(\delta ),\mathcal{M})}{p(\theta |\mathcal{M})}\right]d{y}_{1:N}d\theta$$

$$=\int p(y_{1:N}(\delta )|\mathcal{M})\int p(\theta |y_{1:N}(\delta ),\mathcal{M})\left[\log \frac{p(\theta |y_{1:N}(\delta ),\mathcal{M})}{p(\theta |\mathcal{M})}\right]d\theta d{y}_{1:N}$$

$$={\mathbb{E}}_{y_{1:N}(\delta )|\mathcal{M}}{D}_{KL}[p(\theta |y_{1:N}(\delta ),\mathcal{M})||p(\theta |\mathcal{M})]$$

$$I(y_{1:N}(\delta );\theta |\mathcal{M})=\int \int p(y_{1:N}(\delta )|\theta ,\mathcal{M})p(\theta |\mathcal{M})\left[\log \frac{p(y_{1:N}(\delta )|\theta ,\mathcal{M})}{p(y_{1:N}(\delta )|\mathcal{M})}\right]d{y}_{1:N}(\delta )d\theta$$

$$\underset{\delta }{\mathrm{maximize}}$$

$${\delta }_i\in \{0,1\},$$

$$\text{mathds}{1}^T\delta =N_o.$$

$$I(y_{1:N}(\delta );\theta )$$

$$=H(p(\theta ))-{\mathbb{E}}_{y_{1:N}(\delta )}\left[H(p(\theta |y_{1:N}(\delta )))\right]$$

$$U(\delta )={\mathbb{E}}_{y_{1:N}(\delta )}\left[H(p(\theta |y_{1:N}(\delta )))\right]$$

$$p(\theta |y_{1:N}(\delta ))\approx \mathcal{N}(\theta |\hat{\theta }(y_{1:N}),A_N^{-1}(\hat{\theta }(y_{1:N}),\delta ))$$

$$[A_N(\theta ,\delta ){]}_{ij}=-\frac{{\partial }^2\log p(y_{1:N}(\delta )|\theta )}{\partial {\theta }_i{\theta }_j}$$

$$U(\delta )=-\frac{1}{2}{\mathbb{E}}_{y_{1:N}(\delta )}\left[\log \det A_N(\hat{\theta }(y_{1:N}),\delta )\right]$$

$$p({\theta }_s|y_{1:N},\mathcal{M})=p({\theta }_s|\hat{\theta }(y_{1:N}),\mathcal{M})$$

$$A_N({\theta }_s)=\left[\begin{array}{cc}{B}_N({\theta }_s)& 0\\ 0& C_N(\hat{\sigma })\end{array}\right]$$

$$[B_N({\theta }_s){]}_{ij}=-\left[\begin{array}{c}\frac{{\partial }^2\log p(y_{1:N}|[{\theta }_s,\hat{\sigma }],\mathcal{M})}{\partial {\theta }_{s,i}\partial {\theta }_{s,j}}\end{array}\right]$$

$$C_{N}(\hat{\sigma})=-\dfrac{\partial^{2}\log p(y_{1:N}|[\theta_{s},\sigma],\mathcal{M})}{\partial\sigma\partial\sigma}\Bigg{|}_{\sigma=\hat{\sigma}}$$

$${\hat{\sigma }}^2({\theta }_s)=\frac{1}{N{N}_o}\sum _{n=1}^N\Vert y_n-S_o(\delta )x_n({\theta }_s){\Vert }^2\stackrel{\Delta }{=}J({\theta }_s)$$

$$p(y_{1:N}|\theta =[{\theta }_s,\hat{\sigma }],\mathcal{M})=[2\pi eJ({\theta }_s){]}^{-N{N}_o/2}$$

$$p({\theta }_s|y_{1:N},\mathcal{M})=\frac{\sqrt{\det B_N({\hat{\theta }}_s)}}{(2\pi {)}^{N{N}_o/2}}\exp \left[-\frac{1}{2}({\theta }_s-{\hat{\theta }}_s{)}^T{B}_N({\hat{\theta }}_s)({\theta }_s-{\hat{\theta }}_s)\right]$$

$$B_N({\hat{\theta }}_s)=\frac{N{N}_o}{2J({\hat{\theta }}_s)}{\left[\frac{{\partial }^{2}J({\theta }_s)}{\partial {\theta }_s\partial {\theta }_s}\right]}_{{\theta }_s={\hat{\theta }}_s}$$

$$\log \det A(\hat{\theta }(y_{1:N}),\delta )=\log C_N(\hat{\sigma })+\log \det B_N({\hat{\theta }}_s)$$

$$C_N(\hat{\sigma })=\frac{2N{N}_o}{{\hat{\sigma }}^2({\hat{\theta }}_s)}$$

$$[B_{N}(\hat{\theta}_{s})]_{pq}=\dfrac{1}{\hat{\sigma}^{2}}\sum_{n=1}^{N}\left[\dfrac{\partial x_{n}}{\partial\theta_{s,p}}S_{o}S_{o}^{T}\dfrac{\partial x_{n}^{T}}{\partial\theta_{s,q}}+\epsilon_{n}\left(S_{o}\dfrac{\partial^{2}x_{n}}{\partial\theta_{s,p}\partial\theta_{s,q}}\right)^{T}\right]\Bigg{|}_{\theta_{s}=\hat{\theta}_{s}}$$

$$[B_N({\hat{\theta }}_s){]}_{pq}\approx \frac{1}{{\hat{\sigma }}^2}\sum _{n=1}^N\left[\frac{\partial x_n}{\partial {\theta }_{s,p}}{S}_o{S}_o^T\frac{\partial x_n^T}{\partial {\theta }_{s,q}}\right]$$

Full text

\pubvolume

xx \issuenum1 \articlenumber5

\history \TitleExploiting Convexification for Bayesian Optimal Sensor Placement by Maximization of Mutual Information \AuthorPinaky Bhattacharyya 1,†,‡\orcidA and James Beck 2,* \AuthorNamesPinaky Bhattacharyya and James Beck

\corresCorrespondence: pinaky_b@alumni.caltech.edu; Tel.: +1-626-399-8070

1 Introduction

In system identification and structural health monitoring, a common goal is to inform the parameters governing a predictive model of a structure using relevant data collected from sensors on it. For this purpose, it is often desirable to incorporate a Bayesian model updating framework to quantify uncertainty surrounding the model parameters based on sensor data. (Beck and Katafygiotis, 1998; Beck, 2010; Vanik et al., 2000; Zhang and Au, 2016). For example, in a structural dynamics setting, the sensors would typically measure the physical response of the structure in terms of accelerations, strains, etc. under natural, induced or ambient excitation. The data is then used to inform the parameters of a model, which is typically a linear dynamical finite-element model with viscous damping but need not be. Any available prior information about the model parameters that is available can be incorporated in the Bayesian model updating framework, where prior and posterior probability distributions over parameters of interest are used to quantify uncertainty in their values before and after utilization of the data, respectively.

The optimal sensor placement problem may be viewed as optimally distributing sensors over a structure to maximize the information gained from the data about the model parameters, given a fixed sensor budget (Beck et al., 1998; Papadimitriou et al., 2000; Yuen et al., 2001; Papadimitriou, 2004; Yuen and Kuok, 2015; Sun and Büyüköztürk, 2015; Argyris et al., 2018). For structures that are not simple geometric shapes, which is typically the case, the optimal sensor placement problem over all possible locations becomes difficult to specify. Therefore, only a pre-determined subset of locations is considered for instrumentation, for instance, the degrees of freedom of a finite-element model of the system. The optimization problem therefore becomes an inherently discrete configuration selection problem. In this case, heuristics such as genetic algorithms, greedy selection, swarm optimization etc., (Papadimitriou et al., 2000) (Papadimitriou, 2004) (Sun and Büyüköztürk, 2015) are usually needed in order to obtain a satisfactory sub-optimal solution, as the computational cost of a brute-force combinatorial search becomes prohibitive for realistic problems.

The optimality criterion is defined according to the purpose of collecting the sensor data. In this paper, it is assumed that the data would be used in a Bayesian model updating framework. In this setting, optimally placed sensors would collect data to maximize information gained about the model parameters. Then the mutual information (Cover and Thomas, 2005) between the model parameters and predicted data is a natural choice for the objective function for this optimization problem. Equivalently, the expected entropy of the posterior of the model parameters with respect to the predicted data can be chosen. This entropy-based objective function for optimal sensor placement was introduced in (Beck et al., 1998) (Papadimitriou et al., 2000).

Under the assumptions that imply applicability of the Laplace approximation to the posterior distribution, and assuming sufficiently small prediction errors, the entropy-based optimal sensor placement problem can be simplified to the optimization of the log-determinant of the Fisher information matrix (FIM) of the model parameters (Beck et al., 1998). This gives us a rigorous information-theoretic foundation for one of the proposed choices of the objective function in (Udwadia, 1994), namely, the aforementioned log-determinant, rather than the trace, of the FIM. The main focus of this paper is to show how this log-determinant formulation can lend itself to a fast sub-optimal, and possibly optimal, solution using convex optimization (Boyd and Vandenberghe, 2004) (Joshi and Boyd, 2009), bypassing the need for a discrete optimization over all possible sensor configurations. The method is applied to an illustrative example of a 50-story uniform shear building, where there are over 47 trillion possible sensor configurations and the optimal sensor placement is found using fewer than 100 objective function evaluations.

2 Bayesian model updating

This section lays out the framework for Bayesian model updating for structural dynamics which forms the basis for the optimal sensor placement problem. The measured data on which inference is to be performed is assumed to be available in the form of dynamic test data in the time domain. It is also assumed here that the complete system input can be determined given the parameters. In a typical structural dynamics setting, one observes the acceleration, velocity or displacement time history at a certain number of locations. The locations that are selected, could, for instance, correspond to a subset of the degrees-of-freedom of a finite-element model of the structure.

Denote the uncertain system model parameters by $\theta_{s}$. These could correspond, for instance, to sub-structure stiffnesses, Rayleigh damping parameters, etc. in the finite element model used to predict the response of the structure.

Denote by $y_{n}$ the stochastic prediction of the quantity to be observed at the time instant $t_{n}$. Note that $y_{n}\in\mathbb{R}^{N_{o}}$, where $N_{o}$ is the actual number of observed degrees-of-freedom out of a total of $N_{d}$ potential degrees-of-freedom for observation. The set of stochastic predictions for the observations over the whole duration of the measurements, $[t_{1},t_{N}]$, is denoted by $y_{1:N}$.

The equation relating stochastic predictions of the observations, $y_{n}$, to the deterministic predictions from the model is:

[TABLE]

where the vector of deterministic predictions at each time-step of interest is denoted by $x_{n}(\theta_{s})\in\mathbb{R}^{N_{d}}$; the uncertain prediction-error is denoted by $\epsilon_{n}\in\mathbb{R}^{N_{d}}$; $\delta\in\{0,1\}^{N_{d}}$ is a binary vector that indicates the presence or absence of a sensor at a degree-of-freedom, so that the sum of entries of $\delta$ equals the number of observed DOFs, $N_{o}$; and the sensor selection matrix, $S_{o}(\delta)\in\mathbb{R}^{N_{o}\times N_{d}}$, selects the components that are to be observed from the full stochastic prediction. The uncertain prediction-errors are modeled by a probability distribution that depends on uncertain prediction-error parameters denoted by $\theta_{e}$; for example, using a zero-mean, discrete-time Gaussian white-noise process with a stationary and isotropic co-variance matrix having diagonal entries $\sigma^{2}$, $\theta_{e}$ would just be the scalar $\sigma$. For an anisotropic model for the co-variance matrix without prediction error correlation between different locations, $\theta_{e}$ would be a vector of variances corresponding to the diagonal entries of the co-variance matrix.

Denote the complete parameter vector by $\theta$, which is the collection of uncertain system parameters and prediction-error parameters, i.e., $\theta=[\theta_{s},\theta_{e}]$. Prior information on the parameter values is expressed in terms of a prior distribution, $p(\theta|\mathcal{M})$. The stochastic forward model, typically denoted by $p(y_{1:N}|\theta,\mathcal{M})$, prescribes the probability of observing a particular data set given the parameters. The probability model class, $\mathcal{M}$, is used to indicate a specific pairing of stochastic forward model and prior. The stochastic forward model with actual data $\hat{y}_{1:N}$ substituted for $y_{1:N}$, when viewed simply as a function of the uncertain parameters, $\theta$, is called the likelihood function. Our prediction-error modeling assumptions imply:

[TABLE]

Bayes’ rule implies that the posterior probability density, $p(\theta|y_{1:N},\mathcal{M})$, is proportional to the product of the likelihood function and prior probability density, thereby allowing an update of the prior:

[TABLE]

Here, the normalizing factor $p(y_{1:N}|\mathcal{M})$ is called the evidence for the model class $\mathcal{M}$ provided by the data $y_{1:N}$. It can be determined using the Total Probability Theorem as:

[TABLE]

For a given data set, the parameter vector that maximizes the likelihood function is known as the maximum likelihood estimate (MLE), denoted by $\hat{\theta}_{MLE}$. Similarly, the parameter vector that maximizes the posterior probability density is known as the maximum a posteriori (MAP) estimate of the parameters, denoted by $\hat{\theta}_{\text{MAP}}$.

2.1 The prediction-error model

Based on the Principle of Maximum Information Entropy (Beck, 2010) (Cover and Thomas, 2005), we model the prediction errors as a discrete-time Gaussian process with equal variance at every degree of freedom, and uncorrelated across time and location, so that:

[TABLE]

This implies that knowing the prediction-error at one time and location provides no information a priori about the errors at other times and locations. This Gaussian distribution on the prediction errors gives the maximum entropy (most uncertainty) for a specified zero mean and common variance $\sigma^{2}$ for all times and locations. Combining Equations (1), (2) and (5), the stochastic forward model then becomes:

[TABLE]

The logarithm of this is needed later:

[TABLE]

3 Mutual information-based optimal sensor placement

The mutual information between two stochastic variables is a measure of about how much they imply about each other’s possible values. It is a fundamental information-theoretic quantification of their correlation (Cover and Thomas, 2005). In our case, the mutual information between the stochastic predictions for the observed data, $y_{1:N}$, and the parameter vector, $\theta$, is defined by:

[TABLE]

where $D_{KL}(\cdot||\cdot)$ is the Kullback-Leibler (KL) divergence (the inner integral over $\theta$ in Equaton (10)).

In the definition of mutual information in Equation (8), the dependence of the stochastic prediction of the data $y_{1:N}(\delta)$ on the sensor configuration $\delta$ is made explicit but will sometimes be omitted for convenience. Equations (9) and (10) come from the product rule in probability. These equations show that the mutual information can be expressed as the expected Kullback-Leibler (KL) divergence, or expected relative entropy, of the posterior with respect to the prior.

In the optimal sensor placement framework, we wish to maximize the mutual information between the predicted sensor data and the model parameters over all possible configurations of placing $N_{o}$ sensors at $N_{d}$ locations. This is equivalent to maximizing the expected KL divergence in Equation (11), which implies that the optimal sensor configuration has the largest expected reduction in entropy of the posterior distribution from the prior. This corresponds to the greatest expected information gain about the parameters from the data during Bayesian updating. In other words, considering all sensor configurations, the optimal one would result in the maximum reduction in parameter uncertainty, on average, upon obtaining the data.

The mutual information in Equation (9) can also be expressed in terms of two specified distributions specified by the model class $\mathcal{M}$, i.e., the prior and the stochastic forward model, by using the product rule again:

[TABLE]

The denominator in the logarithm term is what one would associate with the evidence for the model class $\mathcal{M}$, based on given data, as defined in Equation (4). Here, the data is not available and this term is simply a function obtained upon marginalizing out the parameters from the joint on $y_{1:N}(\delta)$ and $\theta$.

Choosing the mutual information between the uncertain parameters and uncertain data as the objective function to be maximized over the sensor configurations, the optimal sensor placement problem becomes:

[TABLE]

Equation (12) can be used for the objective function for the problem in (13). Recall that $\delta$ is a binary $N_{d}$-dimensional vector specifying the sensor configuration over the degrees of freedom in the system and so it has $N_{o}$ unit components and $N_{d}-N_{o}$ zero components. The explicit dependence on $\mathcal{M}$ has been suppressed in the notation in (13).

In principle, to perform the optimization in Equation (13), we need to search over all possible sensor configurations. Instead, we employ a convex optimization technique that greatly reduces the computational complexity of this search.

4 Equivalent entropy-based optimal sensor placement

In this section, an equivalent entropy-based optimal sensor placement problem is formulated and an efficient solution using convex optimization techniques is presented. The Laplace approximation to the posterior in a Bayesian system identification problem is described and used to express the entropy of the posterior.

Using the expression in Equation (9) for mutual information:

[TABLE]

where $H(p(x))=-\int p(x)\log p(x)dx$ is the entropy corresponding to $p(x)$. Here, the prior entropy term is irrelevant in the sensor placement problem and may be discarded to yield an objective function that depends only on the posterior entropy term:

[TABLE]

Thus, the objective function maximization in Equation (13) can be replaced by an equivalent minimization of the expected posterior entropy of the model parameters $\theta$. The idea of choosing the locations of the $N_{o}$ sensors so that the posterior uncertainty of the parameters is minimized was first introduced in (Papadimitriou et al., 2000) (Beck et al., 1998). In these papers, it is shown that this entropy-based criterion can be simplified to a more manageable log-determinant criterion when the following assumptions are made about the probability model class $\mathcal{M}$:

• •

The model class, upon collection of the data, is globally identifiable and the posterior distribution has a single peak (Beck and Katafygiotis, 1998).

• •

The Laplace approximation holds; that is, a sufficient number of data points are measured so that the posterior is approximately Gaussian (Beck and Katafygiotis, 1998).

Applying the Gaussian approximation to the posterior, and using $\hat{\theta}(y_{1:N})$ to denote the MAP value of $\theta$ based on data $y_{1:N}(\delta)$:

[TABLE]

where the precision matrix $A_{N}$ here is given by the negative of the Hessian of the logarithm of the stochastic forward model (Beck and Katafygiotis, 1998):

[TABLE]

The objective function in Equation (15) is now simply related to the entropy of this approximately Gaussian posterior which may be expressed in terms of the log-determinant of its precision matrix as:

[TABLE]

Note that the optimal sensor location problem is to be solved typically before any response data is available from the real structure to be instrumented, so predictions of the data must be made based on a probability model. This can be constrained based on available information about the structure, either from its design or from preliminary tests, in the form of a structural model and a prior distribution over the system parameters, as shown in Section 2. A nice feature of the log-determinant formulation above, is that the dependence on the data of the approximate posterior distribution over the uncertain parameters is only through the optimal parameters, $\hat{\theta}=\hat{\theta}(y_{1:N})$, that is,

[TABLE]

The precision matrix $A_{N}$ in Equation (17) of the multivariate Gaussian in Equation (16) can be partitioned into two parts corresponding to the MAP values of the system parameters $\hat{\theta}_{s}\in\mathbb{R}^{N_{m}}$ and the prediction-error parameter, $\hat{\sigma}>0$:

[TABLE]

where

[TABLE]

and

[TABLE]

The off-diagonal terms in the Hessian of Equation (20) are zero because of our choice of stochastic forward model and prior. Also, in Equations (20) - (22), $\hat{\sigma}^{2}$ for fixed system parameters $\theta_{s}$ is the MAP prediction-error variance, assuming a uniform prior on $\sigma^{2}$ over a large positive interval $(0,\sigma_{\max}^{2})$:

[TABLE]

Equation (23) can be used to re-express the stochastic forward model at the MAP variance parameter:

[TABLE]

Also, the posterior distribution of the system parameters, with the prediction-error parameter marginalized out, can be expressed as:

[TABLE]

where from Equations (7) and (21):

[TABLE]

Furthermore, using Equation (20), the log-determinant in Equation (18) can be expressed as:

[TABLE]

where from Equations (7) and (22):

[TABLE]

Expanding the second-order derivative in Equation (26) using the expression for $J(\theta_{s})$ from Equation (23):

[TABLE]

It is now assumed that the term involving the second derivative in Equation (29) can be neglected. This can be justified if at least one of the following is true: the prediction errors, $\epsilon_{n}$, are small in magnitude; or the deterministic predictions, $x_{n}$, vary slowly with respect to the system parameters and therefore have small second derivatives. The curvature assumption can be assessed before instrumenting the structure, while the prediction-error assumption can only be validated after checking the agreement between the deterministic predictions and the data.

With this assumption, the relevant portion of the Hessian matrix can be expressed in terms of the prediction sensitivity coefficients at the observed degrees of freedom. The resulting approximation is then:

[TABLE]

In addition, notice the dependence on the sensor configuration becomes $S_{o}^{T}S_{o}=\text{diag}(\delta)$, so Equation (30) can be transformed to a double-sum:

[TABLE]

Equation (31) is a linear sum of contributions of terms from each sensor location. For convenience of notation in the formulation of the optimal sensor placement problem, we introduce the $N_{m}\times N_{m}$ matrix $Q=\hat{\sigma}^{2}B_{N}(\hat{\theta}_{s})$, so:

[TABLE]

where

[TABLE]

Having developed an approximate expression for the posterior distribution over the system parameters, we are now in a position to tackle the optimal sensor location problem.

To greatly simplify calculations, we use nominal values $\theta_{s0}$ and $\sigma_{0}^{2}$ for the system and variance parameters respectively, in place of the optimal parameters $\hat{\theta}_{s}$ and $\hat{\sigma}^{2}$. Then, from Equation (25):

[TABLE]

Since the designer is uncertain about the nominal values, a prior distribution, $p(\theta_{0})$, where $\theta_{0}=[\theta_{s,0},\sigma_{0}^{2}]$, is specified for the nominal parameters. Therefore, Equation (34) replaces the dependence of the posterior on unavailable data with prior information that is already available in terms of the nominal parameters. These parameters are either available in the form of structural design information, or in the form of information obtained from preliminary tests of the structure.

Thus, the objective function $U(\delta)$ in Equation (18) becomes, after using Equations (27), (32) and dropping the constant terms and scale factor of $1/2$:

[TABLE]

where the expression for the matrix $Q(\delta,\theta_{s,0})$ is given in Equations (32) and (33).

Thus, as in (Papadimitriou et al., 2000), the entropy-based optimal sensor location problem requires a minimization of the expected posterior entropy of the system parameters over the binary vector, $\delta$, to give the optimal sensor configuration, $\delta^{\ast}$:

[TABLE]

This is a combinatorial optimization problem, that is, we wish to minimize $h(\delta)$ by choosing $N_{o}$ sensor locations among $N_{d}$ possible locations. Depending upon the problem, this could become prohibitively expensive. Each sensor configuration requires one evaluation of the expectation of the log-determinant of a computationally expensive prediction. For example, selecting 20 locations from among 50 locations to instrument gives a number of possible combinations of about 47 trillion.

Clearly, a brute-force search for the optimal solution over all possible sensor configurations will often not be a feasible approach. However, it is usually the case that it is not possible to guarantee an optimal solution without an exhaustive search. Heuristic methods such as genetic algorithms can produce an acceptable sub-optimal value if run for long enough. Incremental greedy algorithms are also guaranteed to produce a sub-optimal value within $(1-1/e)$ of the optimal (Guestrin et al., 2005), relative to the range from the least to the most optimal configuration. Here, we present an alternative approach along the lines of (Joshi and Boyd, 2009) that applies a convex relaxation technique to provide a sub-optimal solution to the problem.

5 Convex relaxation of the combinatorial optimization

In this section, a relaxed version of the optimization problem is set up and argued to be convex. The relevant partial derivatives of the objective function, that is, its gradient and its Hessian matrix, are derived for the purpose of computing the solution. The final step involves replacing the expectation integral over the prior distribution on the uncertain parameters by its Monte Carlo approximation so that the problem may be solved by a generic convex solver.

The original optimization problem in Equation (36) is relaxed to allow the components of the binary vector $\delta$ to take on any real values between zero and one. The problem thus specified turns out to be a convex optimization for which there are efficient algorithms to find the unique solution, as shown in the next subsection. However, we note that the solution to this new optimization problem need not be a binary vector, because it may contain entries between 0 and 1. The solution is still very meaningful, however. If the solution is indeed a binary vector, then it is the optimal solution to the original combinatorial problem. If not, then it is an upper bound to the value of the objective function that also gives guidance to what the optimal solution is likely to be.

We replace the original binary sensor placement vector, $\delta$, by a vector, $z$, in the hypercube $[0,1]^{N_{d}}$, resulting in the relaxed optimization problem:

[TABLE]

where the objective function is a straightforward extension of the previous one in Equation (35):

[TABLE]

by replacing $\delta_{i}$ with $z_{i}$ in the definition of $Q$ is Equation (32). A major advantage of the relaxed problem is that it allows one to use continuous, rather than discrete, optimization packages to solve the problem.

5.1 Convex nature of the relaxed optimization

Equation (37) describes a convex optimization problem where the objective function is convex in $z$, and the equality and inequality constraints are affine. The log-determinant function in Equation (38) is convex in $z$ (Joshi and Boyd, 2009) and the expectation operator preserves the convexity of its argument. Thus, the problem has a unique global optimum. This can be determined computationally, as long as the objective function can be computed at every $z$.

For larger problems, it is useful to efficiently supply the gradient and Hessian of the objective function with respect to $z$. This avoids expensive computations of their finite-difference counterparts. Fortunately, their analytical expressions are tractable and, along with their Monte Carlo approximations, are described in Appendix A.

5.2 Solver for the relaxed optimization

Since the problem in Equation (37) is convex, it may be solved using Newton’s method. In order to apply this in practice, numerical approximations to the objective and its derivatives may be used, as in Equations (68) through (70). The MATLAB routine fmincon was used for the examples in this paper. Since the optimization problem has equality and inequality constraints, a constrained convex optimization solver needs to be employed. The expression for the Newton step is not as straightforward as in the unconstrained case and is given in (Boyd and Vandenberghe, 2004).

6 Application to multistory structures

6.1 Problem description

The problem of optimally placing a fixed number $N_{o}$ of sensors over a structure for system identification is considered. The structure is assumed to be well modeled by linear dynamics with classical viscous damping. For optimal design of the sensor layout, the structure is imagined to be subject to a sinusoidal ground acceleration at the base specified by an amplitude $a_{0}$ and frequency $\omega$. Therefore, the governing differential equation for its displacement vector $x(t)$ as a function of time is given by,

[TABLE]

where $M$, $C$ and $K$ are the $N_{d}\times N_{d}$ mass, damping and stiffness matrices respectively, and $\mathds{1}$ is a vector of 1’s.

In this problem, a nominal natural frequency parameter, the Rayleigh damping parameters, input frequency and input amplitude are all taken to be uncertain, i.e. $\theta_{s}=[\omega_{0},\alpha,\beta,\omega,a_{0}]$. The prediction-error variance $\sigma^{2}$ for a stationary, isotropic prediction error co-variance matrix, is the uncertain prediction-error parameter. The designer specifies a prior distribution on these uncertain parameters.

Before converting this system into its modal co-ordinates, $q_{j}(t)$, with corresponding natural frequency and damping ratio $\omega_{j}$ and $\zeta_{j}$ respectively, the problem is simplified to take into account uncertainties in modeling the system. This is done largely through a substructuring scheme. The assumptions result in a simple expression, relative to the uncertain parameters, for each mode.

The following simplifications are assumed:

$M=mM^{\ast}$ and $K=kK^{\ast}$, where the mass and stiffness matrices are factored into an uncertain scalar and a known matrix part. Scalars $k$ and $m$ will be called the uncertain nominal stiffness and mass of the system respectively. 2. 2.

$\omega_{0}=\sqrt{\dfrac{k}{m}}$ is the uncertain nominal natural frequency parameter. Although $k$ and $m$ are uncertain, only their square-root ratio $\omega_{0}$ is relevant and its uncertainty may be treated explicitly. 3. 3.

$C=\alpha M+\beta K=\alpha mM^{\ast}+\beta kK^{\ast}$, is the Rayleigh damping matrix with uncertain parameters $\alpha$ and $\beta$.

As a consequence of these simplifying assumptions, the following conditions hold:

The modal frequencies $\omega_{j}=c_{j}\omega_{0}$, where $c_{j}^{2}$ are the known eigenvalues of the $(K^{\ast},M^{\ast})$ system. 2. 2.

The eigenvector matrix $\Phi$ that diagonalizes the $(K^{\ast},M^{\ast})$ system also diagonalizes the $(K,M)$ system and can be used to solve the original dynamical system in Equation (39). 3. 3.

$\Phi^{T}M^{\ast}\Phi=\textup{diag}\{\mu_{j}\}$ and $\Phi^{T}K^{\ast}\Phi=\textup{diag}\{\kappa_{j}\}$, such that $\dfrac{\kappa_{j}}{\mu_{j}}=c_{j}^{2}$.

Given these simplifications, a modal solution to the problem may be calculated analytically.

6.2 Modal solution

We let $x(t)=\Phi q(t)$, where vector $q(t)$ contains the modal displacement co-ordinates, so the original equation may be simplified to:

[TABLE]

Under the assumptions in Section 6.1, this simplifies to a decoupled set of modal ordinary differential equations:

[TABLE]

where

[TABLE]

Each modal equation corresponds to a driven, damped oscillator whose solution is:

[TABLE]

The prediction equation for the displacements then becomes,

[TABLE]

6.3 Sensor placement algorithm

In order to solve the relaxed optimal sensor placement problem in this case, the following operations need to be performed:

Determine the eigenvalues $c_{j}^{2}$ and eigenvector matrix $\Phi$ of the $(K^{\ast},M^{\ast})$ system 2. 2.

Choose an initial guess $z_{0}$ for the sensor positions that satisfies the constraints 3. 3.

Generate $N_{k}$ samples of $\theta_{s}=[\omega_{0},\alpha,\beta,a_{0},\omega]$ from the designer-specified prior distribution $p(\theta_{s})$ 4. 4.

For each sample $\theta_{s}^{(k)}$, calculate the gradient of $q_{j}(t)$ with respect to $\theta_{s}$ for $t_{n}=n\Delta t\text{ for }n=1,\ldots,N$ 5. 5.

Calculate the gradient ${\nabla_{\theta_{s}}}x_{i}(t_{n},\theta_{s}^{(k)})$ 6. 6.

Calculate the elementary matrices $Q^{(i)}(\theta_{s}^{(k)})$ for every $i$ and $k$ and store them. 7. 7.

Compute the objective function, its gradient and Hessian at the current vector $z_{m}$ 8. 8.

Update $z_{m}$ to $z_{m+1}$ and repeat steps 7 and 8 until satisfaction of the convergence criterion

6.4 Illustrative example of uniform shear buildings

Consider the problem of placing a fixed number $N_{o}$ of displacement sensors on a $N_{d}$-DOF uniform shear building, as shown in Figure 1.

In this case, the mass and stiffness matrices are given by:

[TABLE]

The prior parameters are taken to be independent and are distributed as:

[TABLE]

These choices are all reasonable and what is expected to be encountered in typical situations, although they have not been chosen to model any specific physical system. While the stochastic forward model is different for structures with a different number of stories, the same prior is used for each structure considered in this example.

Table 1 shows the results of solving the relaxed, convex optimization problem for structures with up to 8 degrees-of-freedom on which up to 2 sensors are to be installed. Once the $Q^{(i)}(\theta_{s}^{(k)})$ matrices have been computed and stored for samples $\theta_{s}^{(k)}$ of $\theta_{s}$, the optimization algorithm always converges to the same solution irrespective of the initial starting point, $z_{0}$. This is expected, since the problem is convex and is guaranteed to have a global minimum. These are exact optimal solutions since they give $z^{\ast}_{i}\in\{0,1\}$ in each case. For each structure (value of $N_{d}$), the sensor configuration remains stable to changes in the number of time-steps $N$ or prior samples $N_{k}$ used (each chosen as either 1000 or 2000). There is a sensor on the roof in every case.

For a more challenging problem we consider the sensor location scheme for a 50-story shear building ($N_{d}=50$) to be instrumented with $N_{o}=20$ sensors. Some typical results are presented in Figures 2 and 3. Since this problem involves a much larger number of degrees-of-freedom, there is a higher variance associated with the computation of the expectation of the log-determinant of the sensitivity matrix using prior samples. For this example, this shows up as a difference in the value of the non-binary sensor configuration at stories 26 through 28. It is to be noted, however, that there does not appear to be any ambiguity upon rounding the numbers to their binary values. For example, no sensors should be placed at stories 26 and 27 in Figure 2, or at story 27 in Figure 3. The results show that all floors above the 27th, except floors 40, 43, 46 and 49, and no floor below the 28th except floor 2 should be instrumented.

As a reminder, an exhaustive search for the optimal configuration would have required over 47 trillion objective function evaluations. The gradient-based method requires fewer than 100 evaluations with a tolerance of $10^{-6}$ in MATLAB’s fmincon function. Of course, heuristic methods might have provided a good sub-optimal solution but in this example, optimal convergence is achieved by our method because the implied binary solution given by rounding the results in Figures 2 and 3 give essentially the same objective value as the optimum for the relaxed optimization problem.

6.5 Comparison of differently chosen sensor configurations

The number of bits of information gain of the optimal configuration, $z^{\ast}$, obtained using the convex technique over some other configurations is of interest here. Table 2 describes the sensor configurations that are considered for the uniform 50-story structure:

The configuration $z_{\text{greedy}}$ is the sensor configuration obtained using a greedy sequential placement algorithm that picks the best sensor location for the $(n+1)^{\text{th}}$ sensor after $n$ sensors are in place, starting with just one sensor. This algorithm is called the forward sequential sensor placement algorithm in (Papadimitriou, 2004). The number of objective function evaluations is $N_{o}(2N_{d}-N_{o})=1600$ in this example, compared with fewer than 100 evaluations for the convex algorithm.

Table 3 shows that it is inefficient to place sensors in the lower stories relative to the higher ones. It turns out in this case that $z_{\text{greedy}}=z^{\ast}$. With the convex relaxation scheme, however, it is immediately understood that the solution is optimal because it came out as a binary result (while the greedy solution is always binary but not necessarily optimal).

The commonly-used sensor configuration given by $z_{\text{common}}$ in Table 2 appears to be significantly sub-optimal for the given model class. This configuration is sometimes popular, however, for considerations beyond the model class being considered; for instance, when it is expected that the behavior of the structure is not well understood, and so not well modeled.

6.6 Variation of optimal configuration with parameter sample set

For a Monte Carlo sample size, $N_{k}$, of 1000 or 2000, the optimal configuration is consistent for smaller structures ($N_{d}<10$). For larger structures ($N_{d}\approx 50$), however, the optimal configuration may differ from sample set to sample set, presumably because the two sample sizes are not large enough for accurate estimation of the expectation. The severity of the problem increases with the size of the structure and the number of parameters used. It turns out that for this example, the majority of the sensor locations agree. Those in agreement are typically found in the upper half of the structure. This may change with the specific model choices made in formulating the problem, however.

Possible solutions to this problem:

• •

Select a sufficiently large sample size to avoid variation between sample sets

• •

Avoid naive Monte Carlo integration and use an integral that provides an estimate with smaller variance. The method would depend on the nature of the integrand. An algorithm similar to simulated annealing could be used to compute the expectation value, such as in (Beck and Zuev, 2013)

• •

Obtain better design information about the uncertain parameters to obtain a possibly more peaked prior with less variation between samples within a sample set, but this is not easy to accomplish

6.7 Variation of optimal configuration with number of observation time-steps

For a fixed parameter sample set, there is a variation in the optimal configuration when the number of observation time steps $N$ is changed from 1000 to 2000 for structures with $N_{d}>10$. Sensor locations near the roof typically remain the same. However, below some story, the sensor locations can be quite different. This may be because of the significantly different amount of information provided about the dynamic response when the duration of measurement is doubled.

7 Concluding remarks

It was noted that the optimal sensor placement problem for system identification using any of the objective functions of (i) mutual information between the predicted sensor data and the system model parameters; (ii) the expected information gain about the parameters from the data; or (iii) the expected posterior uncertainty (entropy) of the parameters, give equivalent results.

Then the original combinatorial entropy-based minimization problem for optimal sensor placement was cast into a relaxed convex optimization problem and solved efficiently to yield sub-optimal, and sometimes optimal, solutions. Locations for which a sensor was unambiguously present or absent were typically obtained, across differing Monte Carlo samples for performing a required expectation and differing number of sample time steps for the data. Some partly ambiguous locations occurred where sample variability led to slightly different results. Ambiguous locations could be resolved using a brute-force combinatorial approach since the problem size for optimizing over just the ambiguous locations would be reduced and manageable.

In future work, it would be beneficial to explore the effects on the optimal sensor locations of spatial correlation between the predicted data at the possible sensor locations. Such a study could build on the preliminary investigation in Chapter IV of (Bhattacharyya, 2016).

\appendixtitles

no \appendixsectionsmultiple

Appendix A Derivatives of the objective function

A.1 Gradient of the objective

[TABLE]

[TABLE]

A.2 Hessian of the objective

[TABLE]

A.3 Numerical approximations

Given $N_{k}$ samples $\theta_{s}^{(k)}$ distributed according to the prior $p(\theta_{s})$ specified by the designer of the sensor installation, the expectation integrals in Equations (38), (62) and (67) may be approximated by their corresponding Monte Carlo estimates.

For the objective function,

[TABLE]

For the gradient of the objective function,

[TABLE]

Finally, for the Hessian of the objective function,

[TABLE]

Some computational effort may be spared by noting that,

[TABLE]

Hence, stored values of $Q^{(i)}(\theta_{s})$ may be used to determine $Q(z,\theta_{s})$.

\externalbibliography

yes