Model-Based Active Source Identification in Complex Environments

TL;DR

This paper introduces a model-based active source identification method using PDEs and robot sensors to efficiently locate sources in complex advection-diffusion environments, demonstrated through simulation and experiments.

Contribution

It presents a novel PDE-constrained optimization approach with dimension reduction and active sensing for robotic source identification in complex environments.

Findings

Effective source localization in complex AD systems.

First practical use of PDEs for robotic source identification.

Method outperforms heuristic approaches in accuracy and efficiency.

Abstract

We consider the problem of Active Source Identification (ASI) in steady-state Advection-Diffusion (AD) transport systems. Unlike existing bio-inspired heuristic methods, we propose a model-based method that employs the AD-PDE to capture the transport phenomenon. Specifically, we formulate the Source Identification (SI) problem as a PDE-constrained optimization problem in function spaces. To obtain a tractable solution, we reduce the dimension of the concentration field using Proper Orthogonal Decomposition and approximate the unknown source field using nonlinear basis functions, drastically decreasing the number of unknowns. Moreover, to collect the concentration measurements, we control a robot sensor through a sequence of waypoints that maximize the smallest eigenvalue of the Fisher Information matrix of the unknown source parameters. Specifically, after every new measurement, a SI problem is solved to obtain a source estimate that is used to determine the next waypoint. We show that our algorithm can efficiently identify sources in complex AD systems and non-convex domains, in simulation and experimentally. This is the first time that PDEs are used for robotic SI in practice.

Figures

Equations

$$-\nabla \cdot (\kappa \nabla c)+\nabla \cdot (c\mathbf{u})-s=0$$

$$c=0$$

$$a(c,v)=\ell (v;s),\forall v\in V,$$

$$a(c,v)\triangleq {\int }_{\Omega }\kappa \nabla c\cdot \nabla vd\Omega +{\int }_{\Omega }v\mathbf{u}\cdot \nabla cd\Omega ,$$

$$\ell (v;s)\triangleq ⟨\ell (s),v⟩\triangleq {\int }_{\Omega }svd\Omega ,$$

$$\mathcal{M}(c;s)\triangleq Ac-\ell (s),$$

$$\chi_{E}({\mathbf{x}})\triangleq\left\{\begin{array}[]{ll}1&{\mathbf{x}}\in E\\ 0&{\mathbf{x}}\notin E\end{array}\right.$$

$$c^m(\mathbf{x})={\chi }_E(\mathbf{x})c(\mathbf{x})(1+ϵ),$$

$${\Vert c-c^m\Vert }_{{\chi }_E}^2\triangleq {\int }_{\Omega }(c-c^m{)}^2{\chi }_{E}d\Omega$$

$$\mathcal{J}(c,s)\triangleq \frac{1}{2}{\Vert c-c^m\Vert }_{{\chi }_E}^2+\tau \mathcal{R}(s).$$

$$\underset{(c,s)\in V\times S}{\min }\mathcal{J}(c,s)s.t.\mathcal{M}(c,s)=0,$$

$$\mathcal{J}(c,s)=\frac{1}{2}{\int }_{\Omega }(c-c^m{)}^2{\chi }_{E}d\Omega +\tau {\int }_{\Omega }sd\Omega .$$

$${\mathbf{C}}_{ij}\triangleq \frac{1}{R}{\int }_{\Omega }{c}_i{c}_{j}d\Omega .$$

$${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_R\ge 0\text{ and }\mathbf{Q}=[{\mathbf{q}}^1{\mathbf{q}}^2\dots {\mathbf{q}}^R];$$

$${\psi }_k=\sum _{i=1}^R{q}_i^k{c}_i.$$

$$\frac{{\sum }_{i=1}^N{\lambda }_i}{{\sum }_{i=1}^R{\lambda }_i}\ge \eta .$$

$$c_d=\mathbf{\psi }\text{ and }{v}_d=\mathbf{\psi }\mathbf{v},$$

$$\phi_{j}({\mathbf{x}};\underline{{\mathbf{x}}}_{j},{\bar{\mathbf{x}}}_{j})\triangleq\left\{\begin{array}[]{ccl}1&&\text{if }\underline{{\mathbf{x}}}_{j}\leq{\mathbf{x}}\leq{\bar{\mathbf{x}}}_{j}\\ 0&&\text{o.w.}\end{array}\right.$$

$$s_d(\mathbf{x})=\sum _{j=1}^M{\beta }_j{\varphi }_j(\mathbf{x};{\underline{\mathbf{x}}}_j,{\bar{\mathbf{x}}}_j),$$

$$\underset{\mathbf{c},\mathbf{p}}{\min }J(\mathbf{c},\mathbf{p})$$

$$s.t.\mathbf{M}(\mathbf{c},\mathbf{p})=0,$$

$${\beta }_j\ge 0,\mathbf{l}\le {\underline{\mathbf{x}}}_j\le {\bar{\mathbf{x}}}_j\le \mathbf{u},$$

$$y_i\approx c_d({\mathbf{x}}^i;\mathbf{p})(1+ϵ),$$

$$y_i\approx c_d({\mathbf{x}}^i;\mathbf{p})+\overline{ϵ}({\mathbf{x}}^i),$$

$$\mathbf{F}=\frac{1}{{\bar{\sigma }}_{\max }^2}\sum _{i=1}^m{\left(\frac{\partial c_d({\mathbf{x}}^i;\mathbf{p})}{\partial \mathbf{p}}\right)}^T\left(\frac{\partial c_d({\mathbf{x}}^i;\mathbf{p})}{\partial \mathbf{p}}\right),$$

$$c_d(\mathbf{x})=\sum _{i=1}^N{c}_i{\psi }_i(\mathbf{x})=\mathbf{\psi }(\mathbf{x})\mathbf{c},$$

$$\frac{\partial c_d(\mathbf{x},\mathbf{p})}{\partial \mathbf{p}}=\mathbf{\psi }(\mathbf{x}){\mathbf{A}}^{-1}{\nabla }_{\mathbf{p}}\mathbf{b},$$

$$\mathbf{S}(\mathbf{p})\triangleq {\mathbf{A}}^{-1}{\mathbf{M}}_{\mathbf{p}},$$

$$\mathbf{F}({\tilde{\mathbf{x}}}_m;\mathbf{p})=\mathbf{S}(\mathbf{p}{)}^T\mathbf{X}({\tilde{\mathbf{x}}}_m{)}^T\mathbf{X}({\tilde{\mathbf{x}}}_m)\mathbf{S}(\mathbf{p}),$$

$${\mathbf{X}}({\tilde{\mathbf{x}}}_{m})=\left[\begin{array}[]{ccc}\psi_{1}({\mathbf{x}}^{1})&\dots&\psi_{N}({\mathbf{x}}^{1})\\ \vdots&\ddots&\vdots\\ \psi_{1}({\mathbf{x}}^{m})&\dots&\psi_{N}({\mathbf{x}}^{m})\end{array}\right],$$

$${\mathbf{x}}^{m+1}=\underset{\mathbf{x}\in \Omega }{\mathrm{argmax}}{\lambda }_{\min }[{\mathbf{F}}_m+\mathbf{S}({\mathbf{p}}_m{)}^T\mathbf{\psi }(\mathbf{x}{)}^T\mathbf{\psi }(\mathbf{x})\mathbf{S}({\mathbf{p}}_m)],$$

$${\Vert {\mathbf{p}}_m-{\mathbf{p}}_{m-1}\Vert }_2\le ϵ$$

Full text

I Introduction

The problem of Source Identification (SI) refers to the estimation of the properties of a source using a set of measurements of a quantity that is generated under the action of that source. The SI problem has various applications ranging from environmental protection to human safety. Locating atmospheric, underground, or underwater pollutants, finding the source of a hazardous chemical leakage, and fire detection are a few examples. In addition, SI can be an important component in higher level tasks such as search and rescue missions and crowd evacuation.

Due to its importance, the SI problem has been investigated extensively in the literature. Many approaches are heuristics that neglect the physics of the problem or model-based methods that are computationally expensive given the resources available to a mobile robot. In this work, we propose an Active Source Identification (ASI) algorithm that, to the best of our knowledge, is the first method to systematically address the SI problem for a mobile robot utilizing the physics of the underlying transport phenomenon. Particularly, we focus on the steady-state Advection-Diffusion (AD) model as the running example, although, the same principals apply to any other transport model, e.g., heat transfer or acoustics.

I-A The Source Identification Problem

The SI problem in AD transport systems is known as chemical plume tracing and odor localization in the robotic literature and has been investigated since the early $80$s. This task often entails three steps: detection, localization, and declaration, and most of the available algorithms focus on the localization stage [1]. The algorithms differ depending on the dispersal mechanism, i.e., diffusion- or turbulence-dominated, and are specialized for the particular types of sensors used to take the measurements. They are often bio-inspired and try to mimic the behavior of different bacteria [2], insects [3], or crabs [4]. Generally, the main idea is to stay in the plume and move upwind, in the concentration ascent direction, or a combination of the two. In the literature, the former approach is called anemotaxis while the latter is called chemotaxis [1].

The authors in [2] propose a controller that combines anemotaxis and chemotaxis to localize a source. This method is compared to our algorithm in Section V. Arguing that gradient based methods can get trapped in local optima and plateaus, [5] proposes a biased random walk strategy for a robotic swarm to localize multiple point sources. The authors in [6] propose fluxotaxis which uses the mass conservation principle to trace a chemical plume using a robotic swarm. In a more recent work [7], a group of mobile agents are controlled to stay in a formation centered in the plume while they move upwind and localize an ethanol source. To localize multiple sources, the authors of [8] construct a statistical model of the discovered sources allowing the robots to find the next source by subtracting the effect of the previous ones.

When the transport phenomenon is turbulent, disconnected and non-smooth concentration patches appear in the medium. In this case, gradient-based approaches are often significantly inaccurate. The authors in [9] propose a gradient-less search strategy, called infotaxis, that maximizes the expected rate of local information gain. The authors in [10] extend the infotaxis strategy to a multi-agent system. This approach shows a behavior that resembles that of a moth, i.e., casting and zigzagging, which amounts to an extensive exploration of the domain and can be energy-inefficient for a mobile sensor.

The above heuristic approaches to SI are often successful in practice but they also suffer from various drawbacks: First, they do not offer a systematic approach that can handle the localization task under a wide range of conditions. Instead, they are specialized for specific scenarios and sensors. Second, these methods often can only localize a single point source or at best multiple point sources and provide no information about the intensity of the sources. Moreover, they declare localization when the robots physically reach the source while in fact it might be unsafe to approach the source in some applications. Finally, often these heuristics are proposed for convex environments and they do not handle obstacles and non-convex domains easily. These limitations can be addressed if the underlying physics is properly incorporated in the formulation, which leads to model-based SI methods. 111Note that some of the heuristics above selectively utilize physical principles but our intention here is a dedicated, systematic formulation. These model-based SI methods are a special class of Inverse Problems (IPs) which have been studied for a long time; see, e.g., [11]. Methods to solve IPs rely on a mathematical model of the underlying transport phenomenon which often is a Partial Differential Equation (PDE) and, in the special case of SI problems, it is linear in the unknown source term.

The literature on model-based SI problems can be classified in different ways based on the state of the problem, the number of sources, and their shape. Generally, transient transport phenomena are more challenging compared to the steady-state ones, but time-dependent measurements are more informative. The localization of a single point source in steady-state is considered in [12] for a semi-infinite domain, whereas the authors in [13] address the problem for transient transport relying on the a priori knowledge of the possible point source locations. SI in the presence of multiple point sources is considered using optimization-based methods. For instance, the work in [14] addresses the detection and localization of multiple such sources using a wireless sensor network.

More general problems that involve sources of arbitrary shapes in arbitrary domains are typically solved numerically using, e.g., the Finite Element Method (FEM). Discretization of a steady-state PDE using the FEM, leads to a linear time-invariant system where the source term acts as a control input. From this perspective, the IP is similar to the problem of input reconstruction. However, one of the main assumptions in this problem is that the number of observations is no less than the dimension of the unknown input [15]. This assumption is violated in IPs which are typically ill-posed. To resolve this issue, regularization techniques can be used [11, 16]. The authors in [17] use the FEM along with total variation regularization to solve the SI problem. Similarly, in our previous work [18], we proposed the Reweighted Debiased $\ell_{1}$ algorithm, which is an iterative sparse recovery approach to the SI problem. Despite generality, numerical methods such as FEM become computationally demanding as the size of the domain grows. Furthermore, unlike heuristic methods, none of the model-based SI methods discussed above, rely on robots to collect the measurements that are needed to solve the problem.

I-B Active Sensing

Optimal measurement collection has been long studied in the robotics literature to solve state estimation problems. Given a probabilistic model of the measurement noise, information-theoretic indices, e.g., covariance [19], Fisher Information Matrix (FIM) [20], different notions of entropy [21], mutual information [22, 23], and information divergence [24], have been used for general robotic planning. For example, given an information distribution, the authors in [20] propose an optimal controller to navigate the robot through an ergodic path. We investigate the performance of this planning method for SI in Section V. A common predicament in applying some of these methods for SI is the need for the posterior distribution of the unknown source parameters. Obtaining this distribution for SI problems requires solving stochastic IPs and is computationally expensive, see, e.g., [25]. This makes the application of optimality indices that require calculation of the expected information gain, e.g., entropy, mutual information, and information divergence, intractable.

Typically in SI problems, the amount of information provided by a measurement depends on the value of the unknown parameters in addition to the measurement location. A common approach to address this point is to combine the path planning for optimal measurement collection with the solution of the SI problem in a feedback loop which leads to Active SI methods. To solve the planning problem, scalar measures of the FIM can be used, as for state estimation; see, e.g., [26, Ch. 2]. The difference is that in SI the unknown parameters cannot be obtained in closed form by a filter update, but instead they are obtained by the solution of an IP that is much more difficult to solve. Specifically, the work in [27] presents trajectory planning for an autonomous robot, utilizing the trace of the FIM, to identify parameters of a transient Advection-Diffusion model under an instantaneous gas release in an infinite domain. Similarly, [28, 29] propose continuous-time optimal control methods that utilize the determinant of the FIM for trajectory planning for IPs with a few unknowns; [28] considers the SI problem in transient state under the assumption that the noisy measurements are taken continuously, while [29] is an extension of [28] for general IPs. In a different approach, the authors of [30] propose an adaptive SI algorithm to localize a single point source emphasizing on path planning in unknown, possibly non-convex, environments. Common in the above literature on ASI is that the proposed methods avoid the solution of complex SI problems by either assuming very simple mathematical models for the SI problem that can be efficiently solved, or by assuming that the solution of the SI problem is provided a priori and the goal of planning is to collect measurements that are then examined to find a source term that matches those measurements. Therefore, these methods do not apply to general SI problems and for this reason they have also not been demonstrated in practice. To the best of our knowledge, our work is the first ASI method that considers active measurement collection for realistic SI problems.

I-C Proposed Method

In this paper, we consider the problem of Active Source Identification in Advection-Diffusion (AD) transport systems in steady-state. As most path planning methods for state estimation, we propose a method that combines SI and path planning in a feedback loop. The difference is that here the estimation problem is not solved by a closed form filter update, but instead it requires the solution of a complex PDE-constrained optimization problem.

Particularly, given a set of noisy measurements, we formulate the SI problem as a variational regularized least squares optimization problem subject to the AD-PDE. To obtain a tractable solution to this problem, we employ Proper Orthogonal Decomposition [31] to approximate the concentration field using a set of optimal basis functions. Moreover, we model the source term using nonlinear basis functions, which decreases the dimension of the parameter space significantly, although at the expense of introducing nonlinearity. Using these parameterizations, we approximate the functional formulation of the SI problem with a low dimensional, nonlinear, constrained optimization problem, which we solve iteratively utilizing the gradient and Hessian information that we explicitly provide. To initialize this nonlinear optimization problem, we rely on the point-source Sensitivity Analysis of the SI objective function [32].

Assuming a small number of measurements are available to initialize the identification process, we determine a sequence of waypoints from where a mobile robot sensor can acquire further measurements by formulating a path planning problem that maximizes the minimum eigenvalue of the FIM of the unknown source parameters with respect to the noisy concentration measurements. The integrated algorithm, alternates between the solution of the SI and path planning problems. In particular, with every new measurement the solution of the SI problem produces a new source estimate, which is used in the path planning problem to determine a new location from where a new measurement should be taken, and the process repeats. By appropriately decomposing the domain, we show that the proposed algorithm can identify multiple sources in complex AD systems that live in non-convex environments.

I-D Contributions

To the best of our knowledge, this is the first model-based ASI framework that has been successfully demonstrated in practice. As discussed in Section I-A, existing literature on ASI includes either model-free bio-inspired heuristic methods or model-based approaches that employ simplifying assumptions, e.g., point sources in infinite domains, to mitigate the complexity of solving real SI problems. Compared to the heuristic approaches, i.e., chemotaxis, amenotaxis, fluxotaxis, and infotaxis, our proposed model-based SI method combines all these bio-inspired behaviors systematically in a general identification framework. Specifically, the concentration readings are explicitly modeled in the least squares objective while the gradient information (chemotaxis), velocity information (anemotaxis), and the first principals (fluxotaxis) are rigorously encapsulated in the AD-PDE. Finally, the information content of the measurements (infotaxis) are incorporated in the solution of the planning problem. On the other hand, compared to model-based approaches that rely on oversimplified models, our method can solve more general and realistic problems. We have shown that our method outperforms existing approaches both in simulation and experimentally.

Like the path planning methods for state estimation, our method combines SI and path planning in a feedback loop. However, as discussed in Section I-B, the difference is that unlike the estimation problem that is solved by a closed-form filter update, the SI problem requires the solution of a complex PDE-constrained optimization problem. The key ideas that enable a tractable solution to this problem are: (a) a suitable integration of model order reduction, point-source Sensitivity Analysis, and domain decomposition methods, (b) a nonlinear representation of the source term that reduces the dimension of the parameter space, and (c) an information theoretic metric to measure the value of measurements for identifying unknown source parameters. The result is a set of techniques, insights, and methodological advancements that show how to efficiently design a model-based SI method that can be implemented onboard robots. Form a technical standpoint, the proposed framework bridges the gap in the rich but disconnected literature on source localization and active sensing that was discussed before; see Figure 1.

A significant contribution of this work is that it is the first to demonstrate applicability of robotic model-based SI methods to real-world problems. Real-world SI problems present major practical challenges related to modeling and estimation of the flow properties, which serve as the data in the AD-PDE, instability of the AD-PDE itself, and the uncertainties that are present in the parameters and boundary conditions; see the discussion in Section V-C. We show that our algorithm is robust to uncertainties and performs well despite various simplifying assumptions made to model the real world; see [33, 34] for more details about these assumptions.

A preliminary version of this work can be found in the conference paper [35]. The conference version only discusses the SI problem in convex domains and many theoretical details are absent due to space limitations; these details are included in Section II. In addition, here we discuss the integrated ASI algorithm that combines SI and planning and we also consider non-convex environments. Finally, we present experimental results that illustrate our method for real-world SI problems.

The rest of this paper is organized as follows. In Section II we introduce the SI problem and discuss the proposed approach to solve it. Section III is devoted to the formulation and discussion of the path planning algorithm. We discuss the integration of the SI and path planning algorithms along with their extension to non-convex domains in Section IV. Section V contains the numerical simulations and experimental results and finally Section VI concludes the paper.

II Source Identification

II-A Advection-Diffusion Transport

Let $\Omega\subset{\mathbb{R}}^{d}$ denote the domain of interest ($1\leq d\leq 3$), and assume the presence of sources is modeled by a nonnegative function, $s:\Omega\rightarrow{\mathbb{R}}_{+}$. 222For the problem considered here, we assume that sources are strictly positive functions. In general, sources can also be negative in the case of sinks. Sinks can appear, e.g., in the presence of chemical reactions that consume a contaminant. The treatment of the problem in that case is similar. Let $c:\Omega\rightarrow{\mathbb{R}}_{+}$ be the measurable quantity, such as concentration, generated by this source function. Moreover, let the velocity at which this quantity is transported via advection be ${\mathbf{u}}\in{\mathbb{R}}^{d}$ and $\kappa\in{\mathbb{R}}_{+}$ be the diffusivity of the medium. Under steady-state assumption and applying a zero-valued Dirichlet condition to the boundaries $\Gamma$ of the domain, we arrive at the following Boundary Value Problem (BVP) [36, ch. 2]

[TABLE]

We consider Dirichlet conditions for the sake of simplicity; more general boundary conditions can be considered without any additional complications [37]. In order for the BVP (1) to have a solution we assume that $s\in L^{2}(\Omega)$, i.e., $s$ is square integrable over $\Omega$, and define the feasible set for the source term as $S=\left\{s\in L^{2}(\Omega)\ |\ s\geq 0\right\}$.

The BVP (1) can be equivalently represented in variational form as follows. Consider the set $V\subset H^{1}_{0}(\Omega)$, i.e., the set of functions that themselves and their first weak derivatives are square integrable and have compact supports. Thus every $v\in V$ satisfies the boundary condition (1b). Multiplying equation (1a) by the trial function $v\in V$, integrating over the domain, and using Green’s theorem, we obtain the variational formulation of the Advection-Diffusion PDE as

[TABLE]

where $a:V\times V\to{\mathbb{R}}$ is a non-symmetric continuous positive-definite bilinear form defined as

[TABLE]

and $\ell(s):V\to{\mathbb{R}}$ is a continuous linear functional defined as

[TABLE]

where the notation $\langle\ell(s),v\rangle$ indicates the operation of $\ell(s)$ on the function $v$. Given $s\in S$, we define the linear functional ${\mathcal{M}}(c;s):V\to{\mathbb{R}}$ as

[TABLE]

where the operator $A:V\to V^{\prime}$ is defined by $\langle Ac,v\rangle=a(c,v),\ \forall v\in V$. The notation $V^{\prime}$ denotes the dual space of $V$, i.e., the space of linear functionals acting on $V$. Using this definition, the VBVP (2) is equivalent to the operator equation ${\mathcal{M}}(c;s)=0$ where ${\mathcal{M}}:V\times S\to V^{\prime}$. Note that the functions $c$ and $v$ in the VBVP (2) have to be differentiable once. Moreover, it can be shown that for $s\in S$ the BVP (1) and VBVP (2) are equivalent and we can use them interchangeably. For further theoretical details, see [37, ch. 8, 9].

II-B The Source Identification Problem

In this section, we formulate the SI problem as a constrained optimization problem subject to the AD transport model (2). Specifically, consider $m$ stationary sensors deployed in the domain $\Omega$ that take measurements of the concentration $c$, and let $E\subset\Omega$ be the set of $m$ compactly supported measurement areas enclosing the sensor locations. 333Note that the compact measurement area around any given sensor can be made arbitrarily small so that this sensing model approximates point measurements. Define, further, the indicator function $\chi_{E}:\Omega\to\left\{0,1\right\}$ for the set $E$ as

[TABLE]

and let $c^{m}:\Omega\rightarrow{\mathbb{R}}_{+}$ be a function that assigns to ${\mathbf{x}}\in\Omega$ the noisy concentration measurement at that location, i.e.,

[TABLE]

where $\epsilon\sim{\cal{N}}(0,\ \sigma^{2})$ and the measurement noise is proportional to the signal magnitude. Then, the SI problem that we consider in this paper consists of determining an estimate $s$ of the true source term ${\bar{s}}$, given a set $E$ of $m$ noisy measurements in the domain $\Omega$, so that the AD model ${\mathcal{M}}(c;s)=0$ defined in (5) predicts the measurements $c^{m}$ as close as possible in the least squares sense.

The main challenges in solving the SI problem arise due to the following two reasons. First, generally the number of measurements $m$ is considerably smaller than the number of parameters that are used to describe the unknown source term. Second, the measurements are contaminated with noise. To address these two challenges, we follow a standard approach and formulate the SI problem as a regularized least squares optimization problem subject to the AD model (2). Let

[TABLE]

be a measure of discrepancy between the measurements and concentration field predicted by the AD model and define the cost functional ${\mathcal{J}}(c,s):V\times S\to{\mathbb{R}}_{+}$ to be optimized by

[TABLE]

In this equation, $\tau$ is the regularization parameter and ${\mathcal{R}}(s)$ is a functional that specifies the characteristics of the source $s$ that is selected as the solution of the SI problem. In this work, we select ${\mathcal{R}}(s)\triangleq\left\|s\right\|_{L^{1}}=\int_{\Omega}\left|s\right|\ d\Omega=\int_{\Omega}s\ d\Omega,$ where the last equality holds since $s$ is nonnegative. This choice of regularization penalizes the size of the source term. Optimization of the objective functional ${\mathcal{J}}(c,s)$ subject to the AD model (2) gives rise to the following problem

[TABLE]

where the functional ${\mathcal{M}}(c,s)$ is defined by (5) and

[TABLE]

To solve the SI problem (9), the gradient of the cost functional ${\mathcal{J}}(c,s)$ is needed. We obtain this gradient using the so called Adjoint Method. This method allows us to solve (9) directly in the reduced space $S$ of source functions rather than in the full space $V\times S$ of the concentration and source functions. This is possible by using the model ${\mathcal{M}}(c,s)=0$ to represent the concentration $c$ as a function of the source term $s$, i.e., $c={\mathcal{F}}(s)$ where ${\mathcal{F}}:S\to V$. 444As discussed in Section II-A such a representation exists and is unique. Using this gradient information we can minimize the cost functional $\bar{{\mathcal{J}}}(s)={\mathcal{J}}({\mathcal{F}}(s),s)$ and determine the source term $s$ that solves the original problem (9). See Appendix A for the details of the Adjoint Method.

II-C Finite Dimensional Approximation

The variables $c$ and $s$ of the optimization problem (9) are functions that live in the infinite dimensional function spaces $V$ and $S$, respectively. Therefore, in order to solve this problem numerically, it is necessary to approximate $V$ and $S$ by finite dimensional subspaces $V_{d}\subset V$ and $S_{d}\subset S$ determined by appropriate sets of basis functions. This approximation allows us to parametrize the concentration and source functions by a finite number of parameters that depend on the basis functions that constitute $V_{d}$ and $S_{d}$.

The key idea to obtain the finite dimensional subspace $V_{d}$ of the concentration function space $V$ is to use Proper Orthogonal Decomposition (POD) to reduce the order of the model. The POD method is easy to implement and gives an optimal set of basis functions that can be readily used in our formulation to parameterize $c$. For a survey of popular model order reduction methods, see, e.g., [38]. At the same time, we use a nonlinear representation of the source term $s$ as a combination of compactly supported tower functions. This representation reduces the dimension of $S_{d}$ drastically, compared to classical approaches that utilize the Finite Element method.

II-C1 Model Order Reduction

To reduce the order of a model using POD we need to solve the AD-PDE (1) for all values of the unknown source term and build a set of basis functions that span the solution of the AD model. We refer to the solutions as the snapshots of the problem. Let $C=\left\{c_{i}({\mathbf{x}})\right\}_{i=1}^{R}$ denote a set of $R$ snapshots obtained by solving the AD-PDE (1) for different realizations of the source term, i.e., each $c_{i}({\mathbf{x}})\in V$ corresponds to a given $s_{i}({\mathbf{x}})\in S$. The objective of POD is to generate a set of optimal basis functions that maximize the averaged projection of the snapshots over these basis functions; see, e.g., [39]. This optimization problem is equivalent to an eigenvalue problem for the covariance matrix ${\mathbf{C}}\in{\mathbb{R}}^{R\times R}$ defined by

[TABLE]

The details of this procedure are presented in Algorithm 1, which yields $V_{d}=\text{span}\left\{\psi_{k}\right\}_{k=1}^{N}$ for $N<R$ where $\psi_{k}$ are the POD basis functions.

In line 2 of this algorithm $\boldsymbol{\Lambda}$ is the diagonal matrix of the eigenvalues.

As shown in [39, thm. 1], the $i$-th eigenvalue $\lambda_{i}$ of matrix ${\mathbf{C}}$ contains the average energy in the $i$-th mode. Moreover for a given number $N<R$ of basis functions, the POD bases have the maximum possible energy and are optimal. Thus, for a given fraction $\eta$, we can select the number $N$ of required bases as the smallest number such that

[TABLE]

II-C2 Parameterization

Using the basis functions $\psi_{k}$ that constitute $V_{d}=\left\{\psi_{k}\right\}_{k=1}^{N}$ we can represent the functions $c$ and $v$ by a finite number of parameters, that can be used for numerical optimization. Specifically, we define

[TABLE]

where $\boldsymbol{\psi}=[\psi_{1}\ \dots\ \psi_{N}]$ and ${\mathbf{c}},{\mathbf{v}}\in{\mathbb{R}}^{N}$.

To parametrize the source function $s$ we follow a different approach. Specifically, we propose a nonlinear representation of this term as a combination of compactly supported tower functions. The motivation for this representation is that each compactly supported source area can be approximately described by a very small number of parameters corresponding to the intensity and shape of the source. In this paper we focus on rectangular sources, although other geometric shapes can also be used for this purpose.

In particular, let $M$ be the number of basis functions used to approximate the source term in the domain $\Omega\subset{\mathbb{R}}^{d}$ and consider two parameters $\left\{\underline{{\mathbf{x}}}_{j},{\bar{\mathbf{x}}}_{j}\right\}$ for each basis function, where $\underline{{\mathbf{x}}}_{j},{\bar{\mathbf{x}}}_{j}\in{\mathbb{R}}^{d}$ and $j\in\left\{1,\dots,M\right\}$. We define the compactly supported tower functions as

[TABLE]

where the inequalities are considered component-wise and $\underline{{\mathbf{x}}}_{j}\leq{\bar{\mathbf{x}}}_{j}$; cf. Figure 2. Then, for practical purposes we can approximate the desired source term $s\in S$ by

[TABLE]

where we require $\beta_{j}\geq 0$ so that $s_{d}\in S$. We denote by ${\mathbf{p}}=(\beta_{1},\underline{{\mathbf{x}}}_{1},{\bar{\mathbf{x}}}_{1},\dots,\beta_{M},\underline{{\mathbf{x}}}_{M},{\bar{\mathbf{x}}}_{M})$ the vector of parameters associated with the source term $s_{d}$. Thus for $\Omega\in{\mathbb{R}}^{d}$, ${\mathbf{p}}\in{\mathbb{R}}^{p}$ where $p=M(2d+1)$.

Substituting the approximations $\left\{c_{d},s_{d}\right\}$ of the concentration and source terms $c$ and $s$ in equation (9), we obtain a finite dimensional counterpart of the SI problem as

[TABLE]

where $j\in\left\{1,\dots,M\right\}$ and ${\mathbf{l}},{\mathbf{u}}\in{\mathbb{R}}^{d}$ are the lower and upper bounds on the coordinates of the domain. 555We can include additional convex constraints to confine the feasible region that might contain the true source. The extension to non-convex domains is considered in Section IV. Moreover, the objective $J:{\mathbb{R}}^{N\times p}\to{\mathbb{R}}_{+}$ is defined as $J({\mathbf{c}},{\mathbf{p}})={\mathcal{J}}(c_{d},s_{d})$ and the finite dimensional model ${\mathbf{M}}:{\mathbb{R}}^{N\times p}\to{\mathbb{R}}^{N}$ is defined as ${\mathbf{M}}({\mathbf{c}};{\mathbf{p}})={\mathcal{M}}(c_{d};s_{d})$.

The optimization problem (17) can be solved by a variety of available nonlinear optimization algorithms. Any such algorithm requires the first and possibly second order information, i.e., the gradient and Hessian of the objective function, as well as a proper initialization since the problem is nonlinear. In Appendix B, we derive explicit expressions for the gradient and Hessian of the objective function in reduced space $S_{d}$. Consequently, only the bound constraints in (17) need to be considered explicitly for numerical optimization. In the same appendix, we also discuss the Sensitivity Analysis (SA) method for the initialization of problem (17).

III Mobile Robot Path Planning

In Section II we developed a way to efficiently solve the SI problem provided a set of stationary measurements is available. In this section we propose a method to plan the path of a robot so that it collects the required measurements in a way that is optimal with respect to a desired optimality measure. Specifically, we employ the Fisher Information Matrix (FIM). Since the concentration field depends nonlinearly on the source parameters ${\mathbf{p}}$ via the AD model (31), derived in Appendix B-A, the FIM depends on the unknown source parameter. Thus, selecting an optimal set of measurements requires feedback from the SI problem (17) and couples the SI and planning problems.

The path planning problem is initialized with an initial set of ${\bar{m}}$ measurements covering the domain of interest, that are necessary to obtain an initial estimate of the source locations and calculate the initial value of the FIM. These measurements can be obtained by a stationary sensor network that can detect the presence or not of a concentration by monitoring the domain of interest. Optimal selection of these measurement locations is beyond the scope of this paper and here we assume that they are given; see, e.g., [40]. After initialization, the robot collects new measurements sequentially so that the minimum eigenvalue of the FIM is maximized, providing in this way worst-case performance guarantees.

Let ${\tilde{\mathbf{x}}}_{m}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{m})$ for $m>{\bar{m}}$ denote the sequence of measurement locations that belong to the set $E$, introduced in Section II-B, and define by ${\mathbf{y}}_{m}\in{\mathbb{R}}^{m}$ the associated vector of measurements, where from equation (7) we have $y_{i}=c^{m}({\mathbf{x}}^{i})$ for $i\in\left\{1,\dots,m\right\}$. Approximating the concentration function $c$ with its finite-dimensional counterpart $c_{d}$ from Section II-C2, we get

[TABLE]

where $\epsilon\sim{\mathcal{N}}(0,\sigma^{2})$ and we include ${\mathbf{p}}$ to emphasize that the concentration and consequently the measurements depend on the value of the source parameters. We can equivalently represent this equation as

[TABLE]

where $\bar{\epsilon}({\mathbf{x}}^{i})\sim{\mathcal{N}}(0,\bar{\sigma}^{2})$ and $\bar{\sigma}({\mathbf{x}}^{i})=\sigma c_{d}({\mathbf{x}}^{i})$. In order to derive a closed form representation of the FIM, we make the conservative assumption that $\bar{\epsilon}({\mathbf{x}}^{i})\sim{\mathcal{N}}(0,\bar{\sigma}_{\max}^{2})$ where $\bar{\sigma}_{\max}=\max_{{\mathbf{x}}\in\Omega}\bar{\sigma}({\mathbf{x}})$. Then, the additive noise $\bar{\epsilon}({\mathbf{x}}^{i})$ is Normal, spatially independent, and identically distributed and the FIM is given explicitly as [26]

[TABLE]

where ${\mathbf{F}}\in{\mathbb{S}}_{+}^{p}$, i.e., ${\mathbf{F}}$ belongs to the space of $p\times p$ symmetric positive-definite matrices and $p=M(2d+1)$ is the number of unknown parameters. The FIM provides a measure of the amount of information that the measurement vector ${\mathbf{y}}$ contains about the source parameters ${\mathbf{p}}$. Note that the information values for independent observations are additive.

Recall from Section II-C1 that we can construct a set of POD bases whose linear combination spans the finite dimensional concentration field $c_{d}$ of the AD model (2) as

[TABLE]

where $\boldsymbol{\psi}=[\psi_{1}\ \dots\ \psi_{N}]$ and ${\mathbf{c}}\in{\mathbb{R}}^{N}$. Given a set of values for the source parameters ${\mathbf{p}}$, we can calculate the coefficients for this linear representation as ${\mathbf{c}}={\mathbf{A}}^{-1}{\mathbf{b}}({\mathbf{p}})$, where matrix ${\mathbf{A}}$ and vector ${\mathbf{b}}$ are defined in equation (30) in Appendix B-A. Thus $c_{d}({\mathbf{x}};{\mathbf{p}})=\boldsymbol{\psi}({\mathbf{x}})\,{\mathbf{A}}^{-1}{\mathbf{b}}({\mathbf{p}}).$ Therefore, we can calculate the desired derivative in the definition of the FIM (18) as

[TABLE]

where $\nabla_{{\mathbf{p}}}{\mathbf{b}}=-{\mathbf{M}}_{{\mathbf{p}}}$ and ${\mathbf{M}}_{{\mathbf{p}}}$ is the derivative of the finite dimensional AD model (31) with respect to ${\mathbf{p}}$ and is given via equation (B-B) in Appendix B-B. Let

[TABLE]

be a function of ${\mathbf{p}}$ with ${\mathbf{S}}\in{\mathbb{R}}^{N\times p}$, and without loss of generality assume $\bar{\sigma}_{\max}=1$. Note that as long as the variance of the noise is constant, its value is irrelevant for planning. Then, from equation (18), we get ${\mathbf{F}}({\tilde{\mathbf{x}}}_{m};{\mathbf{p}})=\sum_{i=1}^{m}{\mathbf{S}}({\mathbf{p}})^{T}\boldsymbol{\psi}({\mathbf{x}}^{i})^{T}\boldsymbol{\psi}({\mathbf{x}}^{i}){\mathbf{S}}({\mathbf{p}})$ or in matrix form

[TABLE]

where

[TABLE]

is the $m\times N$ design matrix.

Given the sequence of waypoints ${\tilde{\mathbf{x}}}_{m}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{m})$ at step $m>{\bar{m}}$ and the corresponding vector of noisy measurements ${\mathbf{y}}_{m}$, we solve the SI problem (17) to obtain the estimation ${\mathbf{p}}_{m}$ of the unknown source parameters at current step. Then, the Path Planning problem consists of finding the next best waypoint ${\mathbf{x}}^{m+1}$ from where if a new measurement is taken, it will maximize the minimum eigenvalue of the FIM. In mathematical terms

[TABLE]

where ${\mathbf{F}}_{m}={\mathbf{F}}({\tilde{\mathbf{x}}}_{m},{\mathbf{p}}_{m})\in{\mathbb{S}}_{+}^{p}$ is a constant FIM, defined by equation (20), that contains the information from the current $m$ measurements. The second term in the right-hand-side of (22) is a rank-one update capturing the information added by measuring at a new location ${\mathbf{x}}$. Given the solution of (22), we use a motion planner to navigate the robot from its current position to the next measurement location ${\mathbf{x}}^{m+1}$ while avoiding obstacles; see Section V-B for more details. Note that since $\Omega\subset{\mathbb{R}}^{d}$, the dimension of (22) is very small which makes it particularly attractive for online implementation on a mobile robot. The proposed planning algorithm is presented in Algorithm 2.

Line 6 corresponds to solving the Next Best Measurement Problem (22). To do so, we reformulate (22) into a Semi-Definite Programming (SDP) problem and solve it using the Sequential SDP (SSDP) method; see Appendix C for details.

IV Active Source Identification

in Complex Domains

The developments of Sections II and III relied on the assumption that the domain of interest $\Omega$ is represented by a set of convex box constraints. Although extension to handle a set of affine constraints defining $\Omega$ or even any other set of convex constraints is straightforward, this is not the case if $\Omega$ is non-convex. To solve the SI problem discussed in Section II in a non-convex domain $\Omega$, we first decompose this domain into convex subdomains. Then, using the Sensitivity Analysis (SA) initialization method discussed in Appendix B-C, we select the largest subdomains that contain the candidate source locations and solve the SI problem (17) in those subdomains. To solve the planning problem discussed in Section III in non-convex domains, we follow a similar approach. Particularly to determine every new waypoint ${\mathbf{x}}^{m+1}$ of the robot, we define a subdomain of $\Omega$ around the initialization point ${\mathbf{x}}_{0}$, given by equation (65) in Appendix C-D, and solve the nonlinear SDP (C-A) in this region. Note that this initialization scheme ensures existence of a local optimum in the selected subdomain and preserves the global convergence property of the SSDP Algorithm 7 for non-convex domains.

Integrating the solution of the SI problem with the planning problem, discussed in Sections II and III, respectively, in a feedback loop and incorporating the proposed domain decomposition method to handle optimization in non-convex domains, we obtain the proposed Active Source Identification (ASI) method. Specifically, given a set of initial measurements, the robot solves the SI problem as discussed in Section II over the subdomains that contain possible source locations as indicated by the SA method in Algorithm 6. Then, given the solution of the SI problem, the robot plans its next measurement according to the procedure developed in Section III, and the process repeats. The proposed ASI algorithm terminates when

[TABLE]

for some $0<\epsilon\ll 1$, where ${\mathbf{p}}_{m}$ is the solution of the SI problem (17) at step $m$. The proposed integrated method is summarized in Algorithm 3 and illustrated in Figure 3.

In lines 3 and 6 of Algorithm 3, $E_{m}$ denotes the set of infinitesimal areas enclosing measurement locations; cf. Section II-B. Note that via successive solutions of the SI problem in line 6, the solver corrects its previous estimates of the source parameters ${\mathbf{p}}_{m}$ taking into account the most recent concentration measurement $y_{m}$. The initialization of the SI problem with the previous solution results in faster convergence.

V Results

In this section we provide numerical simulations and real-world experiments to illustrate the ASI Algorithm 3. We solve the constrained nonlinear optimization problem (17) utilizing the fmincon($\cdot$) function in the MATLAB optimization toolbox that employs an interior-point algorithm which accepts the Hessian-vector multiplication information; cf. Appendix B-B. We also use the clusterdata($\cdot$) function in MATLAB to perform the clustering required in Algorithm 6 of Appendix B-C. Moreover, we solve the SDP (C-A) with a primal-dual method using cvx, a package for specifying and solving convex programs [41]. For more details about interior-point and primal-dual algorithms, see, e.g., [42].

In order to quantify the performance of our method, we report four different error metrics, namely the uncovered source ratio $e_{\text{un}}$, the false detection ratio $e_{\text{fd}}$, the normalized intensity error $e_{\text{int}}$, and the normalized localization error $e_{\text{loc}}$. In mathematical terms $e_{\text{un}}={\left\|{\bar{s}}_{d}-s_{d}\right\|_{\chi_{F}}}/{\left\|{\bar{s}}_{d}\right\|_{L^{2}}}\ \ \text{and}\ \ e_{\text{fd}}={\left\|s_{d}\right\|_{\chi_{\Omega\setminus F}}}/{\left\|{\bar{s}}_{d}\right\|_{L^{2}}},$ where $F$ is the support set of the true source ${\bar{s}}_{d}$, $\chi_{F}$ denotes the indicator function of $F$ defined in equation (6), and $\left\|\cdot\right\|_{\chi_{F}}$ is introduced in (8). The error term $e_{\text{un}}$ measures the fraction of the true source ${\bar{s}}_{d}$ that is left out by the estimated source $s_{d}$ and the error term $e_{\text{fd}}$ considers the parts of the estimated source $s_{d}$ that do not overlap with the true source ${\bar{s}}_{d}$. Note that any value $e_{\text{un}}<1$ indicates an overlap between the true and estimated sources. Finally, the errors $e_{\text{int}}$ and $e_{\text{loc}}$ are defined for a single source as $e_{\text{int}}=\left|\bar{\beta}-\beta\right|/\beta_{\max}$ and $e_{\text{loc}}=\left\|{\bar{\mathbf{z}}}-{\mathbf{z}}\right\|_{2}/l$, where $\beta_{\max}$ is an upper bound on the source intensity, ${\mathbf{z}}\in\Omega$ is the center of the rectangular source support, and $l$ is the characteristic length of the domain $\Omega$. We also define the signal to noise ratio for simulations in dBs as $\text{SNR}=20\log\left({\left\|c^{m}({\mathbf{x}})\right\|_{\chi_{E}}}/{\left\|\epsilon({\mathbf{x}})\right\|_{\chi_{E}}}\right),$ where $\epsilon({\mathbf{x}})$ denotes the noise field.

We study the performance of the ASI Algorithm 3 as a function of the dimensionless Peclet number, which is a measure of the relative dominancy of advection over diffusion and is defined as $Pe={u\,l}/{\bar{\kappa}}$, where $l$ is the characteristic length, $u$ is the magnitude of the inlet velocity, and $\bar{\kappa}$ is the average diffusivity of the medium.

To generate the snapshots for POD Algorithm 1, and to solve problem (17) numerically, we need to obtain the solution of AD model (1) for a given source function. To this end, we discretize the domain $\Omega$ and use the FE method with standard Galerkin scheme [43]. Let $n$ denote the size of the required FE mesh which we generate using cubit [44]. We construct the discrete FE models using an in-house FE code based on the DiffPack C++ library [45]. Moreover, to approximate the first- and second-order derivatives over FE meshes, we use finite difference for structured meshes and polynomial interpolation for unstructured meshes.

We select the thresholding parameter of the SA Algorithm 6 as $\alpha=0.7$; see Appendix B-C. To build the POD basis functions via Algorithm 1, we need to generate snapshots of the AD-PDE (1). Since the relationship between the magnitudes of the source and concentration functions is linear, cf. equation (1a), we cover the domain $\Omega$ with tower functions (15) with unit intensity for this purpose. Finally, we set the regularization parameter in (10) to $\tau=10^{-8}$ and the stopping tolerance in the ASI Algorithm 3 to $\epsilon=10^{-3}$.

V-A Numerical Simulations

In this section, we study the performance of the ASI Algorithm 3 in a large non-convex domain $\Omega$. We assume that the air flows into the domain through the sides, i.e., constant velocity inlet boundary conditions. Then, an in-house fluid dynamics code is utilized to simulate the steady-state velocity in the domain as depicted in Figure 4, where a FE mesh with $n=15034$ points is used.

We consider Peclet numbers $Pe=2.5,Pe=25$, and $Pe=250$. Moreover to reduce the order of the model, we utilize $R=597$ snapshots with $\eta=0.97$ corresponding to $N=154,N=183$, and $N=205$ basis functions for each Peclet number, respectively.

In our first simulation study, we compare the planning method presented in Section III to placement over a lattice and an ergodic placement method that uses the determinant of the FIM as the information metric. Both of these planning approaches employ the algorithm developed in Section II to solve the SI problem (9) and, therefore, the purpose of this comparison is to showcase the relative performance of the planning method proposed here. Specifically, given an information distribution, the ergodic planner aims at designing paths where the time spent by the robot at any region is proportional to the information at that region. This approach is proposed in [46] and used for active sensing in [20]. Specifically, we use the normalized determinant of the FIM (20) as the information distribution in the ergodic planner, where we assign zero information value to the points on obstacles. We implement the controller proposed in [46] and, similar to the ASI Algorithm 3, we update the information distribution at every step as newer estimates of the source parameters become available.

To highlight the advantages of the proposed ASI framework against the heuristic methods discussed in Section I-A, we also compare our algorithm to a heuristic approach that drives the robot along the normalized concentration gradient ascent and upwind directions, as proposed in [2]. The robot uses the initial ${\bar{m}}$ measurements to detect the plume and initializes its path from the highest measured concentration point. The velocity field is known exactly to the robot and the concentration gradient at each point is approximated by taking two additional measurements in orthogonal directions. Since the heuristic approach only provides a location estimate, it is compared to the other methods in terms of $e_{\text{loc}}$.

In the following simulations, we use ${\bar{m}}=28$ initial measurements for the ASI, ergodic, and heuristic methods and set the maximum number of steps to $m_{\max}=42$. The lattice placement uses $m_{\max}=42$ measurements obtained by sensors located on an equidistant grid. For the ergodic and heuristic methods, we use a first-order model for the dynamics of the robot. The results are plotted in Figure 5 where we average over $50$ randomly generated sources.

It can be seen that the proposed planning Algorithm 2 outperforms the lattice method in all cases and performs more consistently in terms of standard deviation. Particularly, the success rate of the ASI algorithm, i.e., the number of instances that the algorithm finds an overlapping source estimate, is considerably higher for all three Peclet numbers. Note that the average false detection error of the lattice placement for $Pe=250$ is smaller than the ASI algorithm but since this method often fails to find an overlapping source estimate for this Peclet number, $e_{\text{fd}}$ only indicates that the falsely detected sources have smaller volumes than the true sources on average.

The performance of the ergodic approach is close to the proposed planning method since it uses a similar optimality index to collect the measurements. Note that as the domain becomes larger, the computation of the information distribution required by this approach becomes expensive rendering this planning method intractable. Furthermore, since the performance of the ergodic method depends on the combination of exploration and exploitation [20], we allow the robot to travel through obstacles and take a measurement every three steps so that it can reach more informative regions of the domain more often. To the contrary, the behavior of the proposed ASI method indicates that given an initial set of measurements, necessary to detect the unknown sources, the most informative measurements are obtained close to the location of sources as opposed to points farther away. Therefore, the better performance of the ASI algorithm, i.e., its smaller and more consistent false detection error values, can be attributed to this fact.

Considering the last subfigure in Figure 5, we observe that the ASI algorithm provides more accurate localization for $Pe=2.5$ and $Pe=25$ but the heuristic method performs better for $Pe=250$. The reason for this is that for very high Peclet numbers, for which advection is the primary means of transport, the reduced order AD model (31) becomes inaccurate resulting in poor localization for model-based methods; see Section V-C. Nevertheless, the heuristic approach does not provide any information about the size or intensity of the source and in the case of multiple sources, it localizes at most one source or fails altogether.

In our second simulation study, we use the same settings as before and consider an AD transport with $Pe=25$ and two sources, specifically, a circular source centered at $(2.5,0.25)$ with radius of $0.08$ and intensity of $0.25$ and a rectangular source with parameters ${\bar{\mathbf{p}}}=(0.2,3.85,3.95,0.8,0.95)$ creating the concentration filed given in Figure 6(a). Since the two sources are not located in one convex domain, decomposing $\Omega$ into convex subdomains and following the procedure described in Section IV is necessary to recover both sources. We note that the ASI algorithm has no a priori knowledge of the number of sources. It solves the problem in $1951$sec in $22$ steps, which amounts to solving $22$ instances of the SI problem (17). Time required to solve the planning problem (22) is negligible. The final error values are $e_{\text{un}}=0.67$ and $e_{\text{fd}}=0.61$ with SNR $=19.02$ dB. The waypoints of the robot are given in Figure 6(a) by white stars. Note the accumulation of the measurements around the high concentration regions of the domain, i.e., the hot spots [22]. Note also that to cover both sources simultaneously, the robot needs to move back and forth between them. We can minimize the travelled distance, by adding a penalty term in the planning problem (22) to encourage more measurements before moving to the next source but this would be suboptimal from an information perspective. A more viable option is to use multiple robots, which is part of our future work. The result of the SA initialization Algorithm 6 and the final solution with the true source overlaid on it are plotted in Figure 6(b).

Note that the performance of the ASI Algorithm 3 only depends on the dimensionless Peclet number. The units of the other quantities are arbitrary as long as they are selected consistently. Particularly, given a unit for concentration $c$, the unit for source term $s$ is defined as concentration per unit time; cf. AD-PDE (1). See the next section for a specific example.

V-B Experimental Results

In this section, we demonstrate the performance of the proposed ASI Algorithm 3 experimentally for the identification of an ethanol source in air. Particularly, we consider a non-convex domain with dimensions $2.2\times 2.2\times 0.4$m3, cf. Figure 7.

We connect a fan to the domain through a duct that blows air into domain with an axial speed of $0.9$m/s and a tangential speed of $0.2$m/s, creating a turbulent flow. We utilize ANSYS-Fluent to obtain the desired flow properties namely, the velocity and diffusivity fields. Since turbulence is a $3$D phenomenon, we build a $3$D mesh of the domain with $1.94\times 10^{6}$ elements. We note that determining the velocity and diffusivity fields for turbulent flow is non-trivial. Turbulent flow is recognized with high Reynolds numbers and is characterized by severe fluctuations in the flow properties. These fluctuations enhance the mixing in the flow and facilitate the transport of the concentration. This enhanced mixing is often encapsulated in an effective turbulent diffusivity which is proportional to effective turbulent viscosity $\mu$ and the proportionality constant is the Schmidt number Sc. Then, the total diffusivity is the sum of laminar $\kappa_{0}$ and turbulent diffusivities [34]. Mathematically,

[TABLE]

where $\rho\in{\mathbb{R}}_{+}$ denotes the density of the transport medium.

We construct a $2$D discretized AD model using the FEM with $n=12121$ grid points located on the plane of the robot concentration sensor at height of $0.27$m. Using the $2$D model instead of a $3$D model is an approximation since we ignore the transport in $x_{3}$-direction but considerably decreases the computational cost. The turbulent flow properties at the nodes of the AD FE mesh are interpolated from the ANSYS-Fluent $3$D model and are given in Figure 8. The corresponding Peclet number is Pe $=213$.

Regarding the diffusivity field, the following points are relevant. For molecular diffusion of ethanol in air, $\kappa_{0}\approx 1.1\times 10^{-5}$m2/s. Then, from Figure 8(b) it can be seen that the turbulent diffusivity is considerably larger than the laminar diffusivity. This contributes to the numerical stability of the AD model by decreasing the Peclet number. To further increase the stability, we add an artificial diffusion lower bounding the total diffusivity (24) by $10^{-3}$m2/s; see the discussion of Section V-C for more details. Given the discretized model, we use $N=900$ basis functions to construct the reduced model (31).

The ethanol source is located at $x_{1}=1.8$m, $x_{2}=1.8$m, and $x_{3}=0.3$m, across from the velocity inlet and releases ethanol at a steady rate. To collect the measurements, we use a custom built differential drive mobile robot equipped with a MiCS-5524 concentration sensor. In order to eliminate the effect of velocity field on the sensor readings, we place the sensor in a confined box and utilize an air pump to deliver air to the sensor with an approximately constant flow rate, cf. Figure 9.

As mentioned in Section I-A, non-smooth concentration patches appear when the flow is turbulent. To ensure the detection of these patches, we need to allow enough instantaneous readings at each location. Furthermore, to minimize the effect of intervals of low concentration between the detections, we calculate the final concentration value as the average of the readings at the highest quartile. In the experiment, we record $1000$ instantaneous readings at each measurement location. We perform the computations off-board and communicate the commands to the robot via radio communications. This allows us to use a very small robot minimizing the interference with the flow field. The localization needed for motion control of the robot, is provided by an OptiTrack motion capture system and a simple controller is implemented for line tracking. We utilize the VisiLibity toolbox to generate the geodesic path between each pair of waypoints given by the planning Algorithm 2, taking into account the obstacle; see [47].

The robot collects ${\bar{m}}=16$ initial measurements; these measurements are shown in Figure 10(a). Note that the pattern of readings are in agreement with the model prediction verifying the overall accuracy of the numerical solutions of the flow and the AD-PDE (1). Figure 10(b) shows the waypoints where we set the maximum number of measurements to $m_{\max}=25$.

The final solution is ${\mathbf{p}}=(3140,1.69,1.76,1.77,1.85)$ resulting in a location error $e_{\text{loc}}=0.03$. This small error is due to the different heights of the planes at which ethanol is released ($x_{3}=0.30$m) and measurements are taken ($x_{3}=0.27$m). This causes the peak to be somewhat displaced downstream. The intensity of the source is predicted to be approximately $3\times 10^{3}$ppm/s for the estimated source area of $5.6\times 10^{-3}$m2. The video of the SI process is given in [48]. It can be seen from the video that the estimation of the source location approaches the true value immediately after the initial measurements are collected and the rest of the measurements correct the solution for the newer information that becomes available. Particularly, as the measurements get closer to the source, the intensity of the estimated source spikes to account for much larger observed concentration measurements.

V-C Discussion

An important predicament in application of model-based ASI Algorithm 3 is handling advection-dominated problems. Given a transport medium, e.g., air, high advection translates to turbulent flow which is non-trivial to model and an active area of research. Currently viable approaches are based on the Reynolds-Averaged Navier-Stokes (RANS) models that provide time-averaged properties. These models often suffice for engineering applications but major assumptions used in their derivation, technicalities pertaining to mesh generation and boundary layer treatment, and oftentimes conflicting predictions from different models affect their predictive ability so that additional experimental studies may be necessary to validate them; see [33] for more details.

On the other hand, advection-dominated AD models are also challenging and an active area of research. The reason is usually the presence of numerical instabilities. In the previous section, we added a constant artificial diffusion which is a common practice in the relevant literature [49]. Nevertheless, this might in general lead to forward solutions that are inconsistent with the solution of the original AD-PDE [50]. Numerous more advanced stabilization techniques exist that artificially introduce diffusion in a consistent manner; see e.g. [51]. Note that very high Peclet numbers, e.g., $Pe\approx 10^{3}$, are reported in the literature for the forward solution of the AD-PDE (1) but solving the Inverse Problem using the AD model is considerably more challenging.

Specifically, instability of the AD model adversely affects the POD method, the SA initialization, and consequently the nonlinear optimization problems (9) and (22). After extensive simulation and experimental studies, we have observed that our method works well for Peclet numbers up to approximately $Pe\approx 250$. While here we employ the standard Galerkin scheme for simplicity, more sophisticated FEM could be employed to improve this bound. Note that an important feature of the ASI Algorithm 3 is that it is highly modular, meaning that different components, i.e., the formulation of the SI problem (9), planning problem (22), model reduction Algorithm 1, and SA initialization Algorithm 6, can be independently improved for better SI performance. For instance in Section V-A we replaced the planning module with the ergodic placement for the purpose of comparison.

Finally, for very low Peclet numbers, i.e., diffusion dominated cases, high concentration regions are local, thus it is possible for the SA method to miss some sources if there are no measurements close enough to those sources. Moreover, in the case of multiple sources, if the intensities differ considerably, the SA technique typically detects the high intensity ones. In these cases, using higher numbers of initial measurements ${\bar{m}}$ and tuning the thresholding parameter $\alpha$ in the SA Algorithm 6 can improve the initialization. We also note that, the SA initialization can be combined with possible prior knowledge about the sources to improve the performance of the proposed ASI algorithm.

VI Conclusions

In this paper we considered the problem of Active Source Identification (ASI) in Advection-Diffusion (AD) transport systems in steady-state. Specifically, given a set of noisy concentration measurements, we developed a novel feedback control scheme that combines a SI problem and a planning problem to guide a mobile robot through an optimal sequence of measurements allowing it to estimate the desired source. We employed model reduction and source parametrization techniques to reduce dimensionality and, therefore, the size of the SI and planning problems, and proposed a domain decomposition method to handle non-convex domains and multiple sources. We illustrated our proposed ASI algorithm in simulation and real-world experiments.

Acknowledgements

The authors would like to thank Eric Stach, Yihui Feng, and Yan Zhang for their help with the design of the experimental setup.

Appendix A Adjoint Method

In this section, we discuss the details of the Adjoint Method to obtain the gradient of the SI problem (9) that we formulated in Section II-B. Particularly, the Lagrangian function of this constrained optimization problem is given as

[TABLE]

where $w\in V^{\prime\prime}$ is the adjoint variable. From reflexivity of the Hilbert space $V$, we get $V^{\prime\prime}=V$. Then, referring to the definition of the AD model (5), we have

[TABLE]

Thus, we can rewrite the Lagrangian as

[TABLE]

where $w\in V$ is the adjoint variable.

In what follows, we use the notion of a Gâteaux derivative to differentiate the Lagrangian (25); see, e.g., [37, sec. 9.4].

Definition A.1** (Gâteaux derivative)**

A functional ${\mathcal{T}}:V\to{\mathbb{R}}$ on a normed space $V$ is Gâteaux-differentiable at $c\in V$ if there exists an operator $D_{c}{\mathcal{T}}:V\to V^{\prime}$ defined by

[TABLE]

for all $h\in V$. We use the two notations $\langle D_{c}{\mathcal{T}},h\rangle$ and $\langle{\mathcal{T}}^{\prime}_{c},h\rangle$ interchangeably whenever one of them is clearer.

The Adjoint Method consists of the following three steps that yield an organized procedure for the calculation of the desired gradient; see, e.g., [52, sec. 4]. First, in order to satisfy the AD constraint in the SI problem (9), we set the Gâteaux derivative of the Lagrangian (25) with respect to the adjoint variable $w$ and in an arbitrary direction $v$ equal to zero. The bilinear form $a(c,w)$ and the functional $\ell(w;s)$ are the terms in the Lagrangian that contain $w$. Gâteaux differentiating $a(c,w)$ with respect to $w$ we get

[TABLE]

where we have used linearity of the bilinear operator $a(c,w)$ in each argument. Similarly, Gâteaux differentiating $\ell(w;s)$ with respect to $w$ we get $\langle D_{w}\ell(w;s),v\rangle=\ell(v;s)$. Therefore, the first equation of the Adjoint Method is given as

[TABLE]

Note that this equation is identical to the VBVP (2) and for the function $c$ satisfying this equation, i.e., $c={\mathcal{F}}(s)$, we get ${\mathcal{L}}(c,s,w)=\bar{{\mathcal{J}}}(s)$. Consequently, we can differentiate the Lagrangian (25) to get the desired derivative $\bar{{\mathcal{J}}}^{\prime}_{s}$.

Since $c={\mathcal{F}}(s)$, in order to calculate $D_{s}{\mathcal{L}}(c,s,w)$ we need the derivative ${\mathcal{F}}^{\prime}_{s}$. We can avoid this calculation by setting the Gâteaux derivative of the Lagrangian (25) with respect to the concentration $c$ equal to zero for any arbitrary direction $h$. The two terms containing $c$ are the objective functional ${\mathcal{J}}(c,s)$ and the bilinear form $a(c,w)$. From Definition A.1, the Gâteaux derivative of ${\mathcal{J}}(c,s)$ with respect to $c$ can be calculated explicitly using equation (10) as

[TABLE]

Moreover, similar to the previous case the Gâteaux derivative of the bilinear form $a(c,w)$ with respect to $c$ is given by $\langle D_{c}a(c,w),h\rangle=a(h,w)=a^{*}(w,h),$ where $a^{*}(w,h)$ is the adjoint operator of the bilinear form $a(h,w)$. Therefore, the second equation of the Adjoint Method is given as

[TABLE]

Because of the appearance of the adjoint operator, this equation is called the adjoint equation and the procedure of calculating the desired gradient is referred to as Adjoint Method. Given the concentration $c$ obtained from (26), the solution of equation (28) yields the corresponding adjoint variable $w$.

From the definition of the Lagrangian (25), for the functions $c$ and $w$ satisfying equations (26) and (28), we have $D_{s}{\mathcal{L}}(c,s,w)=\bar{{\mathcal{J}}^{\prime}_{s}}$. Thus, we can calculate the desired gradient of the objective functional $\bar{{\mathcal{J}}}(s)$ with respect to the source term $s$ in a given direction $q$ by Gâteaux differentiating the Lagrangian (25) as $\langle{\mathcal{L}}^{\prime}_{s},q\rangle=\langle{\mathcal{J}}^{\prime}_{s},q\rangle-\langle\ell^{\prime}_{s}(w;s),q\rangle.$ Combining equations (26) and (28) with this equation, we summarize the Adjoint Method to calculate the gradient of $\bar{{\mathcal{J}}}(s)$ with respect to $s$ in a given direction $q$ as:

[TABLE]

Appendix B Numerical Solution of the

Source Identification Problem

B-A First Order Information

In Appendix A we discussed the Adjoint Method to obtain the gradient of the SI problem (9) when the variable $s$ is a function that lives in the infinite dimensional function space $S$. Here, we employ the approximations $V_{d}$ and $S_{d}$ defined in Section II-C to obtain a finite dimensional counterpart of the Adjoint Method equations (29) that is needed to solve the finite dimensional SI problem (17) numerically.

First, we substitute the finite dimensional representations (14) into equation (29a) to get

[TABLE]

Writing the equations for all $i\in\left\{1,\dots,N\right\}$ in matrix form, we obtain the following linear system of equations

[TABLE]

where ${\mathbf{A}}\in{\mathbb{R}}^{N\times N}$ and ${\mathbf{b}}$ is a fixed vector for a given ${\mathbf{p}}$. Using equation (30), we define the finite dimensional model in equation (17) explicitly as

[TABLE]

As explained in Section II-A, the AD model (2) has a unique solution that translates to the invertibility of matrix ${\mathbf{A}}$ in (31).

Similar to approximations (14), we can write $w_{d}=\boldsymbol{\psi}\,{\mathbf{w}}\text{ and }\,h_{d}=\boldsymbol{\psi}\,{\mathbf{h}}$ where ${\mathbf{w}},{\mathbf{h}}\in{\mathbb{R}}^{N}$. Substituting these definitions into the adjoint equation (29b), we get

[TABLE]

where the derivative $\langle{\mathcal{J}}^{\prime}_{c},\cdot\rangle$ is defined by equation (27). Again writing the equations for all $i\in\left\{1,\dots,N\right\}$ in matrix form, we obtain

[TABLE]

where the transpose sign appears in (32) because of the adjoint operator in the equations.

Given values for the source parameters ${\mathbf{p}}$, the linear systems (30) and (32) can be used to obtain the corresponding concentration ${\mathbf{c}}$ and adjoint variable ${\mathbf{w}}$. This information can then be used in (29c) to calculate the desired gradient $\nabla_{{\mathbf{p}}}{\bar{J}}$ of the objective function ${\bar{J}}({\mathbf{p}})=\bar{{\mathcal{J}}}(s_{d})$ with respect to ${\mathbf{p}}$. In order to simplify the notation and without loss of generality, we assume a single source in a $2$D domain given by $s_{d}({\mathbf{x}})=\beta\,\phi({\mathbf{x}};\underline{{\mathbf{x}}},{\bar{\mathbf{x}}}),$ where $\underline{{\mathbf{x}}}=(\underline{x}_{1},\underline{x}_{2})$ and ${\bar{\mathbf{x}}}=({\bar{x}}_{1},{\bar{x}}_{2})$. Substituting the approximations (14) in the Lagrangian (25), we get ${\mathcal{L}}(c_{d},s_{d},w_{d})={\mathcal{J}}(c_{d},s_{d})+a(c_{d},w_{d})-\ell(w_{d};s_{d}).$ To obtain the finite dimensional counterpart of equation (29c), we need to take the derivative of this Lagrangian with respect to the parameters ${\mathbf{p}}$ of the source term $s_{d}$. The terms that contain $s_{d}$ are ${\mathcal{J}}(c_{d},s_{d})$ and $\ell(w_{d};s_{d})$. For the objective functional ${\mathcal{J}}(c_{d},s_{d})$ from equation (10), the only part involving $s_{d}$ is the regularization term $\int_{\Omega}s_{d}\ d\Omega=\beta\ ({\bar{x}}_{1}-\underline{x}_{1})({\bar{x}}_{2}-\underline{x}_{2}).$ From this expression we can calculate the derivatives of ${\mathcal{J}}(c_{d},s_{d})$ with respect to ${\mathbf{p}}$, e.g.,

[TABLE]

For the functional $\ell(w_{d};s_{d})$, substituting $s_{d}$ into the definition (4) we get $\ell(w_{d};s_{d})=\int_{\underline{x}_{1}}^{{\bar{x}}_{1}}\int_{\underline{x}_{2}}^{{\bar{x}}_{2}}\beta\ w_{d}({\mathbf{x}})\ dx_{2}\ dx_{1}.$ The derivative with respect to $\beta$ is straightforward and for the other parameters we use the Leibniz rule, e.g.,

[TABLE]

Then by equation (29c), combining the two derivatives for $\underline{x}_{1}$ we get $\partial{\bar{J}}/\partial\underline{x}_{1}=\partial{\mathcal{J}}/\partial\underline{x}_{1}-\partial\ell/\partial\underline{x}_{1}.$ The other derivatives can be calculated exactly the same way and we get the following values for the desired gradient

[TABLE]

where ${\bar{J}}({\mathbf{p}})=\bar{{\mathcal{J}}}(s_{d})$ and ${\mathbf{p}}=(\beta,\underline{x}_{1},\underline{x}_{2},{\bar{x}}_{1},{\bar{x}}_{2})$. The process for calculating the desired gradient $\nabla_{{\mathbf{p}}}{\bar{J}}$ given a set of values for the parameters ${\mathbf{p}}$ is described in Algorithm 4.

Note that if there are multiple sources, i.e., if $M>1$, then we calculate the gradients for each basis function separately. This follows from the rule for differentiating sums. Moreover if $\Omega\subset{\mathbb{R}}^{3}$, we can exactly follow the same steps to calculate the gradient.

B-B Second Order Information

Including second order information in the optimization algorithm can make the solution of the SI problem (17) more efficient and accurate. Such information can be in the form of the Hessian ${\mathbf{H}}=\nabla_{{\mathbf{p}}{\mathbf{p}}}{\bar{J}}$ of the objective function itself, or in the form of a Hessian-vector product ${\mathbf{H}}{\mathbf{v}}$, for some vector ${\mathbf{v}}$, that is used in the optimization algorithm; see, e.g., [42, ch. 7]. The procedure to calculate the Hessian-vector multiplication is an attractive choice for large-scale problems but we use it here since it provides an organized approach to incorporate the AD model (31) into the Hessian calculations. Specifically, using the finite dimensional approximation of the Lagrangian (25) given as

[TABLE]

we can devise a procedure to calculate the product ${\mathbf{H}}{\mathbf{v}}$ for a given vector ${\mathbf{v}}$. The details of this derivation can be found in [53] and it results in the following equation

[TABLE]

where the subscripts denote differentiation and the process to calculate the vectors ${\mathbf{h}}_{1}$ and ${\mathbf{h}}_{4}$ is explained in Algorithm 5.

In what follows, we discuss all the second order terms needed in Algorithm 5 starting with the derivative of the AD model (31) with respect to the parameters denoted by ${\mathbf{M}}_{{\mathbf{p}}}$. Recalling equation (30) and using the Leibniz rule, row $i$ of matrix ${\mathbf{M}}_{{\mathbf{p}}}\in{\mathbb{R}}^{N\times 5}$ is given as

[TABLE]

Using equation (31), the derivative of the AD model with respect to ${\mathbf{c}}$ is given as ${\mathbf{M}}_{{\mathbf{c}}}={\mathbf{A}}.$ In addition, since there are no terms containing the multiplication of the concentration and source parameters, ${\mathbf{c}}$ and ${\mathbf{p}}$, in $L({\mathbf{c}},{\mathbf{p}},{\mathbf{w}})$, we have

[TABLE]

Finally we need to calculate the second order derivatives of the Lagrangian with respect to each of ${\mathbf{c}}$ and ${\mathbf{p}}$. Note that from equation (34), $\nabla^{2}_{{\mathbf{c}}{\mathbf{c}}}L=\nabla^{2}_{{\mathbf{c}}{\mathbf{c}}}J$ and the value $\langle{\mathcal{J}}^{\prime}_{c},\psi_{i}\rangle$ is basically the directional derivative in the direction $\psi_{i}$ or the derivative $\partial J/\partial c_{i}$. Thus we can Gateaux differentiate equation (27) once more to get the element in row $i$ and column $j$ of $\nabla^{2}_{{\mathbf{c}}{\mathbf{c}}}L$ as

[TABLE]

where $i,j\in\left\{1,\dots,N\right\}$. Note that this expression is independent of the parameters ${\mathbf{p}}$ and can be calculated offline.

In order to calculate $\nabla^{2}_{{\mathbf{p}}{\mathbf{p}}}L$ note that the terms $J({\mathbf{c}},{\mathbf{p}})$ and ${\mathbf{w}}^{T}{\mathbf{b}}$ in the Lagrangian (34) contribute to this derivative. The calculation for $\nabla^{2}_{{\mathbf{p}}{\mathbf{p}}}J$ can be done by differentiating the result of Section B-A for $\nabla_{{\mathbf{p}}}J$ once more to get

[TABLE]

For the second term we have ${\mathbf{w}}^{T}{\mathbf{b}}=\ell(w_{d};s_{d})$, since the Lagrangians (25) and (34) are equivalent. Thus, we can differentiate the expression for $\nabla_{{\mathbf{p}}}\ell$ from Section B-A once more, using the Leibniz rule, to get $\nabla^{2}_{{\mathbf{p}}{\mathbf{p}}}\ell$ as is shown in equation (B-B).

Putting the two terms given by equations (39) and (B-B) together, we have

[TABLE]

Notice that we basically have differentiated equation (B-A) once more in this process.

The case of multiple sources only affects the terms ${\mathbf{M}}_{{\mathbf{p}}}$ and $\nabla^{2}_{{\mathbf{p}}{\mathbf{p}}}L$ given by equations (B-B) and (46), respectively. Since the source term $s_{d}$ defined in equation (16) is the summation of nonlinear basis functions, for ${\mathbf{M}}_{{\mathbf{p}}}$ we need to append more columns using equation (B-B) corresponding to each basis function. On the other hand, for $\nabla^{2}_{{\mathbf{p}}{\mathbf{p}}}L$ we have to add blocks of matrices given by equation (46) corresponding to each basis function to the diagonal of $\nabla^{2}_{{\mathbf{p}}{\mathbf{p}}}L$.

B-C Initialization

Appropriate initialization is critical for the solution of nonlinear optimization problems, such as (17), since otherwise the solution can get trapped in undesirable local minima. In this paper, we employ a result on the point-source Sensitivity Analysis (SA) of the SI cost functional, presented in [32], for initialization of our method. The idea is to determine the sensitivity of the objective functional ${\mathcal{J}}(c,s)$ to the appearance of a point source in $\Omega$, i.e., we calculate the derivative of the objective with respect to the point-source term. The regions with highest sensitivity represent the potential areas where the support of the true source function ${\bar{s}}$ is nonzero. Note that by linearity of the AD-PDE (1), we only need to consider the infinitesimal deviations of the point-source from zero for a source-free domain, i.e., we calculate the derivative for the constant source function $s=0$.

In [32] it is shown that the adjoint variable is a measure of the sensitivity of the cost functional to these infinitesimal changes. Thus given the set of measurements $E$ introduced in Section II-B, we can obtain an approximation to the source locations via a solution of the adjoint equation. Specifically, we solve ${\mathbf{A}}^{T}{\bar{\mathbf{w}}}=-{\bar{\mathbf{d}}}$ with ${\bar{d}}_{i}=\int_{\Omega}c^{m}\ \psi_{i}\ d\Omega$ for $i\in\left\{1,\dots,N\right\}$, to get the desired finite dimensional adjoint function as ${\bar{w}}_{d}=\boldsymbol{\psi}\,{\bar{\mathbf{w}}}$. Then an approximate localization of the source is obtained through thresholding as

[TABLE]

where ${\bar{w}}^{\min}_{d}=\min_{{\mathbf{x}}\in\Omega}{\bar{w}}_{d}({\mathbf{x}})$ and $\alpha\in(0,1)$.

The thresholding parameter $\alpha$ determines the size of the support of ${\hat{w}}_{d}({\mathbf{x}})$ and thus, the number of compact regions that indicate candidate source locations. In order to separate these compact regions, we utilize the Single Linkage Agglomerative Clustering (SLAC) algorithm; see, e.g., [54]. Specifically, given the nodal values ${\hat{\mathbf{w}}}_{d}$ of ${\hat{w}}_{d}({\mathbf{x}})$ over the FE-mesh, we cluster the nonzero nodal values into sets $C_{k}$ for $k\in\left\{1,\dots,K\right\}$. Then, we initialize the SI problem (17) by placing a basis function with a small area at the point with highest sensitivity, given by equation (47), in each cluster; see Algorithm 6 for details.

Appendix C Sequential Semi-Definite Programming

for the Next Best Measurement Problem

C-A The Next Best Measurement Problem

In this section we discuss the details of the numerical solution for the path planning developed in Section III. Let ${\mathbf{F}}({\mathbf{x}})={\mathbf{S}}^{T}[{\mathbf{X}}^{T}{\mathbf{X}}+\boldsymbol{\psi}({\mathbf{x}})^{T}\boldsymbol{\psi}({\mathbf{x}})]{\mathbf{S}}$ denote the FIM at step $m$, where to simplify notation we have dropped the subscripts. Then introducing an auxiliary variable $z$ we can rewrite the optimization problem (22) as

[TABLE]

where $\lambda_{i}$ denotes the $i$-th eigenvalue of the FIM. Problem (C-A) can be equivalently written as

[TABLE]

where the notation $\succ$ denotes a matrix inequality. The optimization problem (C-A) is a nonlinear Semi-Definite Program (SDP) that can be solved using nonlinear optimization techniques; see, e.g., [55]. In this paper we employ the Sequential SDP (SSDP) method which is the extension of sequential quadratic programming; see, e.g., [56]. Defining

[TABLE]

where $0<\bar{\epsilon}\ll 1$ is a very small positive number added to eliminate the strict inequality constraint, we can rewrite problem (C-A) in standard form as

[TABLE]

The Lagrangian corresponding to problem (C-A) is given as

[TABLE]

where $\boldsymbol{\Lambda}\geq{\mathbf{0}}$ is the Lagrange multiplier matrix and the inner-product of two $r\times t$ real matrices ${\mathbf{B}}$ and $\boldsymbol{\Lambda}$ is defined as

[TABLE]

Note that ${\mathbf{B}}(z,{\mathbf{x}}):{\mathbb{R}}^{d+1}\to{\mathbb{S}}^{p}$ in (50b) is a negative-semidefinite symmetric matrix function. Since the Karush-Kuhn-Tucker (KKT) optimality conditions of the nonlinear SDP (C-A) are locally identical to the second-order approximation around any point $({\bar{z}},{\bar{\mathbf{x}}},\bar{\boldsymbol{\Lambda}})$, we can solve a sequence of convex SDPs to build the solution of the nonlinear problem (C-A) iteratively. Under certain conditions that are satisfied for the functions in (50), the SSDP approach converges to a local minimum of the nonlinear SDP (C-A); see, e.g., [56].

Specifically, at each iteration $k$, we construct a second-order convex approximation of (C-A) at point $(z_{k},{\mathbf{x}}_{k},\boldsymbol{\Lambda}_{k})$ as

[TABLE]

where ${\mathbf{v}}_{k}=(z_{k},{\mathbf{x}}_{k})$ is the primal variable at iteration $k$ and ${\mathbf{d}}=(d_{z},{\mathbf{d}}_{{\mathbf{x}}})$ is a vector of directions, where $d_{z}\in{\mathbb{R}}$ and ${\mathbf{d}}_{{\mathbf{x}}}\in{\mathbb{R}}^{d}$ are the directions corresponding to $z_{k}$ and ${\mathbf{x}}_{k}$, respectively. The matrix ${\mathbf{H}}_{k}$ is a positive semidefinite approximation of the Hessian $\nabla_{{\mathbf{v}}{\mathbf{v}}}^{2}L(z_{k},{\mathbf{x}}_{k},\boldsymbol{\Lambda}_{k})$ of the Lagrangian (52) and $D_{{\mathbf{v}}}{\mathbf{B}}({\mathbf{v}}_{k})[{\mathbf{d}}]$ is the directional derivative of the matrix function (50b) at point ${\mathbf{v}}_{k}$ and direction ${\mathbf{d}}$ that is used to linearize the matrix inequality constraint around the current iterate ${\mathbf{v}}_{k}$. This quantity along with the Hessian of the Lagrangian are derived in Appendix C-B. We assume that the domain $\Omega$ is convex and defined by a set of affine constraints so that the linear constraint ${\mathbf{x}}_{k}+{\mathbf{d}}_{{\mathbf{x}}}\in\Omega$ can be directly incorporated in the SDP (C-A). This assumption holds for the box constrained domain that we considered in the SI problem (17).

The solution of the SDP (C-A) denoted by ${\mathbf{d}}_{k}\in{\mathbb{R}}^{d+1}$ determines the descent direction for the nonlinear problem (C-A). Using this solution, we update the primal variables as

[TABLE]

where $\alpha_{k}$ is a step-size whose selection is explained in Appendix C-C. Note that by the last constraint in the SDP problem (C-A), we implicitly assume that the maximum step-size is equal to one, i.e., $\alpha_{\max}=1$. We update the dual variable $\boldsymbol{\Lambda}_{k+1}$ directly as the optimal dual of the tangent problem (C-A). The details of the SSDP to solve the optimization problem (22) are presented next.

C-B First and Second Order Information

To solve problem (C-A) we need the gradient and Hessian information. For this, we first define the expressions $D_{{\mathbf{v}}}{\mathbf{B}}({\mathbf{v}})[{\mathbf{d}}]$ and $\nabla_{{\mathbf{v}}{\mathbf{v}}}^{2}L(z,{\mathbf{x}},\boldsymbol{\Lambda})$ that appear in the convex second-order SDP (C-A). For the first term, we have

[TABLE]

where $d_{i}$ is the $i$-th element of the vector of directions ${\mathbf{d}}$ and

[TABLE]

for $i\in\left\{1,\dots,d+1\right\}$. The operator $D_{{\mathbf{v}}}{\mathbf{B}}({\mathbf{v}}):{\mathbb{R}}^{d+1}\to{\mathbb{S}}^{p}$ is linear in ${\mathbf{d}}$ and $D_{{\mathbf{v}}}{\mathbf{B}}({\mathbf{v}})[{\mathbf{d}}]\in{\mathbb{S}}^{p}$. Therefore the corresponding constraint in the SDP (C-A) is a linear matrix inequality.

For the second term, i.e., the Hessian $\nabla_{{\mathbf{v}}{\mathbf{v}}}^{2}L(z,{\mathbf{x}},\boldsymbol{\Lambda})$ of the Lagrangian (52), since the objective function $f(z,{\mathbf{x}})$ defined by equation (50a) is linear, we have

[TABLE]

where

[TABLE]

and

[TABLE]

Recalling the definition of the matrix function ${\mathbf{B}}(z,{\mathbf{x}})$, given in equation (50b), we calculate the required derivatives (56) and (61) for, e.g., the $2$D case in which ${\mathbf{x}}=(x_{1},x_{2})$. These derivatives then are used in equation (C-A) to build quadratic SDPs that we solve sequentially to find the local optimum of the nonlinear SDP (C-A).

For the first order derivatives used in equation (55), we have

[TABLE]

The value for ${\mathbf{B}}^{(3)}$ is exactly the same as ${\mathbf{B}}^{(2)}$, except that the differentiation variable is $x_{2}$. Similarly for the second-order derivatives used in equation (60), we have

[TABLE]

[TABLE]

${\mathbf{B}}^{(3,3)}$ can be calculated exactly the same way. Finally, for the cross-derivative we have

[TABLE]

After calculating the Hessian $\nabla_{{\mathbf{v}}{\mathbf{v}}}^{2}L(z,{\mathbf{x}},\boldsymbol{\Lambda})$ of the Lagrangian (52), we construct a positive-definite approximation ${\mathbf{H}}$ of it so that the SDP (C-A) is strictly convex with a unique global minimizer ${\mathbf{d}}$. Such an approximation of ${\mathbf{H}}$ can be obtained in different ways; see, e.g., [57]. Here, we add a multiple of the identity matrix so that the minimum eigenvalue is bounded from zero by a small amount $\delta$, i.e., we set

[TABLE]

where $\mu=\max(0,\delta-\lambda_{\min}(\nabla_{{\mathbf{v}}{\mathbf{v}}}^{2}L))$. The positive-definite matrix ${\mathbf{H}}$ is the closest to the Hessian $\nabla_{{\mathbf{v}}{\mathbf{v}}}^{2}L$ measured by the induced Euclidean norm. Note that since the Hessian is a low dimensional matrix, i.e., $\nabla_{{\mathbf{v}}{\mathbf{v}}}^{2}L\in{\mathbb{S}}^{d+1}$, we can easily calculate its minimum eigenvalue.

C-C Step-Size Selection

Necessary for the solution of the nonlinear SDP (C-A) is an effective line-search strategy that connects the successive solutions of the quadratic SDPs (C-A). In this paper, we utilize the results from [58] to select an appropriate step-size $\alpha_{k}$ for the iterations of the SSDP defined by equation (54). The final SSDP algorithm to solve the nonlinear SDP (C-A) is presented in Algorithm 7.

In this algorithm, we define the penalty function $\theta_{\gamma}({\mathbf{v}})$ for the selection of the step-size $\alpha_{k}$ as

[TABLE]

where $\gamma>0$ is the penalty parameter, $\lambda_{\max}({\mathbf{B}})_{+}=\max\left\{0,\lambda_{\max}({\mathbf{B}})\right\}$, and the functions $f({\mathbf{v}})$ and ${\mathbf{B}}({\mathbf{v}})$ are defined in equation (50). The upper bound $\Delta_{k}$ on the directional derivative $\theta^{\prime}_{\gamma_{k}}({\mathbf{v}}_{k};{\mathbf{d}}_{k})$ of the penalty function in a direction ${\mathbf{d}}_{k}$ is given as

[TABLE]

where ${\mathbf{H}}_{k}$ is the positive-definite approximation given by equation (62) and we have included the index $k$ to emphasize that we use the dual variable $\boldsymbol{\Lambda}_{k+1}$ to select the step-size $\alpha_{k}$. The upper bound $\Delta_{k}$ is used in order to satisfy the Armijo condition in the backtracking line-search corresponding to line 8 in Algorithm 7. For theoretical details see [58].

Note that since the domain of interest $\Omega$ is represented by a set of affine constraints that require no further linearization, the constraint ${\mathbf{x}}\in\Omega$ does not appear in the penalty function (63). Essentially, the constraint ${\mathbf{x}}\in\Omega$ is never violated and thus we do not penalize it in (63).

C-D Initialization and Stopping

Since the eigenvalue optimization problem (22) is nonlinear, appropriate initialization of Algorithm 7 is critical to obtain a reasonable solution. Moreover, addition of a new measurement, reshapes the objective function (22) and makes it flat around that measurement location. In other words, adding more measurements in that vicinity does not provide more information about the unknown source parameters compared to farther locations. Therefore without global knowledge of the objective function, the algorithm gets trapped in undesirable local minima where the objective function does not change no matter how many measurements are taken in that region.

In order to generate new informative measurements, we sample the objective function of the Next Best Measurement Problem (22), $g({\mathbf{x}})=\lambda_{\min}[{\mathbf{F}}({\mathbf{p}})+{\mathbf{S}}({\mathbf{p}})^{T}\boldsymbol{\psi}({\mathbf{x}})^{T}\boldsymbol{\psi}({\mathbf{x}}){\mathbf{S}}({\mathbf{p}})],$ over a coarse set of points ${\mathbf{z}}_{i}\in{\mathbb{R}}^{d}$ from the FE mesh, where $i\in\left\{1,\dots,Z\right\}$ for some $Z\ll n$ and $n$ is the number of FE grid points. Then to initialize the primal variable ${\mathbf{v}}_{0}=(z_{0},{\mathbf{x}}_{0})$ in Algorithm 7 for step $m+1$ of the robot, we calculate the values of the objective function $g_{i}=g({\mathbf{z}}_{i})$ over this set of points and we set

[TABLE]

where $z_{0}$ is the auxiliary variable introduced in (C-A). Note that each evaluation of the function $g({\mathbf{x}})$ amounts to solving a minimum eigenvalue problem for a $p\times p$ matrix where $p$ is the number of unknown parameters. The computational cost of such sampling procedure is comparable to a single backtracking line-search step of Algorithm 7 in line 8.

Finally, to determine whether Algorithm 7 has reached a local minimum we evaluate bounds on the gradient of Lagrangian (52), the nonlinear matrix inequality constraint violation, and the complementarity condition as follows

[TABLE]

See [58] for theoretical results supporting this selection of stopping criteria.