PDE-based Dynamic Density Estimation for Large-scale Agent Systems

Tongjia Zheng, , Qing Han, and Hai Lin This work was supported in part by the National Science Foundation under Grant IIS-1724070, and Grant CNS-1830335, and in part by the Army Research Laboratory under Grant W911NF-17-1-0072.Tongia Zheng and Hai Lin are with the Department of Electrical Engineering, University of Notre Dame, Notre Dame, IN 46556, USA (e-mail: tzheng1@nd.edu, hlin1@nd.edu.). Qing Han is with the Department of Mathematics, University of Notre Dame, Notre Dame, IN 46556, USA (e-mail: Qing.Han.7@nd.edu.).

Abstract

Large-scale agent systems have foreseeable applications in the near future. Estimating their macroscopic density is critical for many density-based optimization and control tasks, such as sensor deployment and city traffic scheduling. In this paper, we study the problem of estimating their dynamically varying probability density, given the agents’ individual dynamics (which can be nonlinear and time-varying) and their states observed in real-time. The density evolution is shown to satisfy a linear partial differential equation uniquely determined by the agents’ dynamics. We present a density filter which takes advantage of the system dynamics to gradually improve its estimation and is scalable to the agents’ population. Specifically, we use kernel density estimators (KDE) to construct a noisy measurement and show that, when the agents’ population is large, the measurement noise is approximately “Gaussian”. With this important property, infinite-dimensional Kalman filters are used to design density filters. It turns out that the covariance of measurement noise depends on the true density. This state-dependence makes it necessary to approximate the covariance in the associated operator Riccati equation, rendering the density filter suboptimal. The notion of input-to-state stability is used to prove that the performance of the suboptimal density filter remains close to the optimal one. Simulation results suggest that the proposed density filter is able to quickly recognize the underlying modes of the unknown density and automatically ignore outliers, and is robust to different choices of kernel bandwidth of KDE.

Index Terms:

Large-scale systems, estimation, Kalman filtering, stochastic systems

I Introduction

Efficiently estimating the continuously varying probability density of the states of large-scale agent systems is an important step for many density-based optimization and control tasks, such as sensor deployment [1] and city traffic scheduling. In such scenarios, the agents (such as UAVs) are built by task designers and follow the specified control commands, which means their dynamics are known. While density estimation is a fundamental problem extensively studied in statistics, we are particularly interested in the case where the samples (i.e. the agent states) are governed by known dynamics (which can be nonlinear and time-varying). We study how to take advantage of their dynamics to obtain efficient and convergent density estimates in real time.

In general, density estimation algorithms are classified as parametric and non-parametric. Parametric algorithms assume that the samples are drawn from a known parametric family of distributions with a fixed set of parameters, such as the Gaussian mixture models [2]. Performance of such estimators rely on the validity of the assumed models, and therefore they are unsuitable for estimating an evolving density. In non-parametric approaches, the data are allowed to speak for themselves in determining the density estimate. As a representative, kernel density estimation (KDE) [2] has been widely used for estimating the global density in the study of swarm robotic systems [3, 4, 5]. It is known that the performance of KDE largely depends on a smoothing parameter called the bandwidth [2]. When the samples are stationary, many bandwidth selection techniques have been proposed, such as cross-validation [2], adaptive bandwidth [6] and the plug-in technique [7]. Such techniques are heuristic in general, since any predefined optimality of bandwidth selection requires certain information of the unknown density. They are also not suitable for dynamic estimation which requires the algorithm to be computationally efficient and adapt to the density evolution. For dynamic density estimation, many adaptations of KDE have been proposed in the domain of data stream mining [8]. In such problems, the dynamics of the density are usually unknown, and therefore little can be claimed about the convergence of the algorithm.

To estimate the global distribution of large-scale agents with known dynamics, an alternative way is to formulate it as a filtering problem, for which there exists a large body of literature, such as the celebrated Kalman filters and their variants [9, 10], the more general Bayesian filters and their Monte Carlo approach, also known as particle filters [11]. Unfortunately, all these methods are known to suffer from the curse of dimensionality, especially for large-scale nonlinear systems. A potential solution is to use the so-called consensus filters [12, 13, 14], which obtain local distribution estimates of the agent system using local Kalman or Bayesian filters, and then compute the global distribution by averaging local estimates through some consensus protocols. However, few conclusions are made concerning the relationship between the true distribution and the estimated global distribution. Moreover, stability analysis is difficult when the agents’ dynamics are nonlinear and time-varying.

In summary, considering efficiency and convergence requirements, existing methods are unsuitable for estimating the time-varying density of large-scale agent systems. This motivates us to propose a dynamic and scalable density estimation algorithm that can perform online and take advantage of the dynamics to guarantee its convergence. Specifically, we show that the agents’ density is governed by a linear partial differential equation (PDE), called the Fokker-Planck equation, which is uniquely determined by the agents’ dynamics. KDE is used to construct a noisy measurement of the density (the state of this PDE). The measurement noise is approximately “Gaussian” (more precisely, asymptotically Gaussian when the agents’ population tends to infinity), so that infinite-dimensional Kalman filters can be used to design density filters. It turns out that the covariance of measurement noise depends on the true density, for which approximating the covariance is required and the density filter becomes suboptimal. We then use the notion of input-to-state stability to prove the stability of the suboptimal density filter and the associated operator Riccati equation.

Our contributions contain the following aspects: (i) By using a density filter, we can (largely) circumvent the problem of bandwidth selection; (ii) The density filter takes advantage of the agents’ dynamics to improve its density estimation, in the sense of minimizing the covariance of estimation error; (iii) The density filter is proved to be convergent and is scalable to the population of agents; (iv) All the results hold even if the agents’ dynamics are nonlinear and time-varying.

The rest of the paper is organized as follows. Section II introduces some preliminaries. Problem formulation is given in Section III. Section IV is our main results, in which we present a density filter and then study its stability/optimality. Section V performs an agent-based simulation to verify the effectiveness of the density filter. Section VI summarizes the contribution and points out future research.

II Preliminaries

II-A Infinite-Dimensional Kalman Filters

The Kalman filter is an algorithm that uses the system’s model, known control inputs and sequential measurements to form a better state estimate [15]. It was later extended to infinite-dimensional systems, represented using linear operators [16, 17]. Suppose the signal $x(t)$ and its measurement $y(t)$ , both in a Hilbert space, are generated by the stochastic linear differential equations

	$\displaystyle dx=A(t)xdt+B(t)q(t)dv,\quad x(t_{0})=x_{0},\quad E[x_{0}]=0$
	$\displaystyle dy=C(t)xdt+r(t)dw,\quad y(t_{0})=y_{0}$

where $dv$ and $dw$ are infinite-dimensional Wiener processes with incremental covariance operators $V$ and $W$ , respectively [18]. Assume $\operatorname{Cov}[v(t),w(\tau)]=0$ and $E[\langle v(t),w(\tau)\rangle]=0$ for all $t\neq\tau$ . Denote $Q(t)=q(t)Vq^{*}(t)$ and $R(t)=r(t)Wr^{*}(t)$ . The Kalman filter is an algorithm that minimizes the covariance of the estimation error, given by [16]

d\hat{x}=A(t)\hat{x}dt+L(t)[y-C(t)\hat{x}]dt

where $L(t)=P(t)C^{*}(t)R^{-1}(t)$ is the called the optimal Kalman gain, and $P(t)$ is the solution of the Riccati equation

\dot{P}=AP+PA^{*}-PC^{*}R^{-1}CP+BQB^{*}

with $P(t_{0})=\operatorname{Cov}[x_{0},x_{0}]$ .

II-B Input-to-state stability

Input-to-state stability (ISS) is a stability notion widely used to study stability of nonlinear control systems with external inputs [19]. Roughly speaking, a control system is ISS if it is asymptotically stable in the absence of external inputs and if its trajectories are bounded by a function of the magnitude of the input. To define the ISS concept for infinite-dimensional systems we need to introduce the following classes of comparison functions [19].

	$\displaystyle\mathcal{K}$	$\displaystyle:=\{\gamma:\mathbb{R}_{+}\to\mathbb{R}_{+}\|\gamma\text{ is continuous and strictly}$
		$\displaystyle\quad\quad\text{increasing},\gamma(0)=0\}$
	$\displaystyle\mathcal{L}$	$\displaystyle:=\{\gamma:\mathbb{R}_{+}\to\mathbb{R}_{+}\|\gamma\text{ is continuous and strictly }$
		$\displaystyle\quad\quad\text{decreasing with }\lim_{t\to\infty}\gamma(t)=0\}$
	$\displaystyle\mathcal{KL}$	$\displaystyle:=\{\beta:\mathbb{R}_{+}\times\mathbb{R}_{+}\to\mathbb{R}_{+}\|\beta\text{ is continuous},\beta(\cdot,t)\in\mathcal{K},$
		$\displaystyle\quad\quad\beta(r,\cdot)\in\mathcal{L},\forall t\geq 0,\forall r>0\}.$

Definition 1

(ISS [20]). Consider a control system $\Sigma=(X,U,\phi)$ consisting of normed linear spaces $(X,\|\cdot\|_{X})$ and $(U,\|\cdot\|_{U})$ , called the state space and the input space, endowed with the norms $\|\cdot\|_{X}$ and $\|\cdot\|_{U}$ respectively, and a transition map $\phi:\mathbb{R}_{+}\times X\times U\to X$ . The system is said to be ISS if there exist $\beta\in\mathcal{KL}$ and $\gamma\in\mathcal{K}$ , such that

\|\phi(t,x_{0},u)\|_{X}\leq\beta(\|x_{0}\|_{X},t-t_{0})+\gamma(\sup_{t_{0}\leq\tau\leq t}\|u(\tau)\|_{U})

holds $\forall x_{0}\in X$ , $\forall t\geq t_{0}$ and $\forall u\in U$ . It is called locally input-to-state stable (LISS), if there also exists constants $\rho_{x},\rho_{u}>0$ such that the above inequality holds $\forall x_{0}:\|x_{0}\|_{X}\leq\rho_{x},\forall t\geq t_{0}$ and $\forall u\in U:\|u\|_{U}\leq\rho_{u}$ .

III Problem formulation

This paper studies the problem of estimating the dynamically varying probability density of large-scale stochastic agents. The dynamics of the agents are assumed to be known and satisfy the following stochastic differential equations

dX_{i}=\boldsymbol{v}(X_{i},t)dt+\boldsymbol{\sigma}(X_{i},t)dB_{t},\quad i=1,\dots,n,

(1)

where $X_{i}\in\mathbb{R}^{N}$ represents the state of the $i$ -th agent (e.g. positions), $\boldsymbol{v}=(v_{1},\ldots,v_{N})\in\mathbb{R}^{N}$ is the deterministic dynamics (e.g. velocity fields), $B_{t}$ is an $M$ -dimensional standard Wiener process, and $\boldsymbol{\sigma}=[\sigma_{jk}]\in\mathbb{R}^{N\times M}$ represents the stochastic dynamics. We also assume the collection of states $\{X_{i}(t)\}_{i=1}^{n}$ are observable.

The probability density $p(x,t)$ of agent states $\{X_{i}(t)\}_{i=1}^{n}$ is known to satisfy the Fokker-Planck equation [21]:

\begin{array}[]{cl}&\displaystyle\frac{\partial p(x,t)}{\partial t}=-\sum_{i=1}^{N}\frac{\partial}{\partial x_{i}}[v_{i}(x,t)p(x,t)]\\ &\quad\quad\quad\quad\quad\displaystyle+\sum_{i=1}^{N}\sum_{j=1}^{N}\frac{\partial^{2}}{\partial x_{i}\partial x_{j}}[D_{ij}(x,t)p(x,t)],\\ &p(\cdot,0)=p_{0},\end{array}

(2)

where $p_{0}$ is the initial density and

D_{ij}(x,t)=\frac{1}{2}\sum_{k=1}^{M}\sigma_{ik}(x,t)\sigma_{jk}(x,t).

If the agent states are confined within a bounded domain $\Omega\subseteq\mathbb{R}^{N}$ , one can impose a reflecting boundary condition

\mathbf{n}\cdot(\boldsymbol{g}-\boldsymbol{v}p)=0,\quad\text{on }\partial\Omega,

(3)

where $\boldsymbol{g}=(\sum_{j=1}^{N}\frac{\partial}{\partial x_{j}}(D_{1j}p),\ldots,\sum_{j=1}^{N}\frac{\partial}{\partial x_{j}}(D_{Nj}p))$ , $\partial\Omega$ is the boundary of $\Omega$ and $\mathbf{n}$ is the outward normal to $\partial\Omega$ .

Remark 1

Note that the dynamics of agents’ density (2) are uniquely determined by their individual dynamics (1) and are therefore known. This relationship holds even if (1) is nonlinear and time-varying. Therefore, the density filter to be presented applies to a very wide class of systems. It is worth pointing out that PDE (2) is always linear even if (1) is not. (It is time-varying if (1) is time-varying.)

Now, we can formally state the problem to be solved as follows:

Problem 1

Given the density dynamics (2) and agent states $\{X_{i}(t)\}_{i=1}^{n}$ , we want to estimate their density $p(x,t)$ .

IV Main Results

In this section, we design a density filter to estimate the state (i.e. the density) of (2) by combining KDE and infinite-dimensional Kalman filters. Specifically, we use KDE to construct a noisy measurement of the unknown density and show that the measurement noise is approximately “Gaussian”, which enables us to use infinite-dimensional Kalman filters to design a density filter. Stability analysis is also presented.

IV-A Density Filter Design

In this section, we use KDE to construct a measurement with Gaussian noise and design a density filter.

KDE is a non-parametric way to estimate an unknown probability density function (pdf) [2]. The agents’ states $\{X_{i}(t)\}_{i=1}^{n}\subseteq\mathbb{R}^{d}$ at time $t$ can be seen as a set of $n$ independent samples drawn from the pdf $p(x,t)$ . Fix $t$ and denote $f(x)=p(x,t)$ . The density estimator is given by

f_{n}(x)=\frac{1}{nh^{d}}\sum_{i=1}^{n}K\Big{(}\frac{1}{h}(x-X_{i})\Big{)},

(4)

where $K(x)$ is a kernel function [2] and $h$ is the bandwidth, usually chosen as a function of $n$ such that $\lim_{n\rightarrow\infty}h=0$ and $\lim_{n\rightarrow\infty}nh=\infty$ . The Gaussian kernel is frequently used due to its infinite order of smoothness, given by

K(x)=\frac{1}{(2\pi)^{d/2}}\exp\Big{(}-\frac{1}{2}x^{\intercal}x\Big{)}.

It is known that the $f_{n}(x)$ is asymptotically normal and that $f_{n}(x_{i})$ and $f_{n}(x_{j})$ are asymptotically uncorrelated for any $x_{i}\neq x_{j}$ . These two properties are summarized as follows.

Lemma 1

(Asymptotic normality [22]) Under conditions $\lim_{n\rightarrow\infty}h=0$ and $\lim_{n\rightarrow\infty}nh=\infty$ , if the bandwidth $h$ tends to zero faster than the optimal rate, i.e.,

h^{*}=o\left(\frac{1}{n}\right)^{1/(d+4)}.

then as $n\to\infty$ , we have

\sqrt{nh^{d}}(\hat{f}(x)-f(x))\rightarrow\mathcal{N}\Big{(}0,f(x)\int[K(u)]^{2}du\Big{)}.

(5)

Lemma 2

(Asymptotic uncorrelatedness [22]) Let $x_{i}$ and $x_{j}$ be two distinct continuity points of $f$ . Then under $\lim_{n\rightarrow\infty}h=0$ , as $n\rightarrow\infty$ , the asymptotic covariance of $f_{n}(x_{i})$ and $f_{n}(x_{j})$ satisfies

nh^{p}\operatorname{Cov}[f_{n}(x_{i}),f_{n}(x_{j})]\rightarrow 0.

In practice, $n$ is finite. Hence, the kernel density estimator is biased (i.e. $E[\hat{f}_{n}(x)]\neq f(x)$ ). Also, $f_{n}(x_{i})$ and $f_{n}(x_{j})$ are correlated. However, one can always choose a smaller bandwidth $h$ to reduce the bias and the covariance [22, 2], which means that $f_{n}(x)-f(x)$ can be made approximately Gaussian with independent components when $n$ is large. It is known that the performance of KDE largely depends on its bandwidth [2]. Any predefined optimality of the bandwidth selection requires certain information of the unknown density $f$ . For example, the one that balances the estimation bias and variance depends on the second-order derivatives of $f$ . An excessively large (small) bandwidth may cause the problem of oversmoothing (undersmoothing). For the density filter to be presented, we tend to choose a smaller bandwidth. We recall that “filters” essentially combine past outputs to produce better estimates, which can ease the problem of undersmoothing. In this regard, by using a density filter, we can largely circumvent the problem of optimal bandwidth selection. Simulation studies will show that the performance of the proposed density filter is robust to different choices of bandwidth.

We are now ready to present a density filter using the infinite-dimensional Kalman filters. We rewrite the PDE of density evolution (2) in the form of an evolution equation and use KDE to construct a noisy measurement $y(t)$ :

	$\displaystyle\dot{p}(t)$	$\displaystyle=A(t)p(t)$		(6)
	$\displaystyle y(t)$	$\displaystyle=p_{\text{KDE}}(t)=p(t)+w(t)$		(6)

where $A(t)=-\sum_{i=1}^{N}\frac{\partial}{\partial x_{i}}(v_{i}\cdot)+\sum_{i=1}^{N}\sum_{j=1}^{N}\frac{\partial^{2}}{\partial x_{i}\partial x_{j}}(D_{ij}\cdot)$ is a linear operator, $p_{\text{KDE}}(t)$ represents a kernel density estimator using the states $\{X_{i}(t)\}_{i=1}^{n}$ at time $t$ , and $w(t)$ is the measurement noise which is approximately Gaussian with covariance operator $R(t)=k\operatorname{diag}(p(t))$ where $k>0$ is a constant depending on $n$ and $h$ .

According to Section II-A, the optimal density filter can be designed as

\dot{\hat{p}}=A(t)\hat{p}+L(t)(y-\hat{p}),\quad\hat{p}(t_{0})=p_{\text{KDE}}(t_{0}),

(7)

where $L(t)=P(t)R^{-1}(t)$ is the optimal Kalman gain and $P(t)$ is a solution of the following operator Riccati equation

\dot{P}=AP+PA^{*}-PR^{-1}P.

(8)

The Riccati equation (8) actually depends on the unknown state/density $p(t)$ because $R(t)=k\operatorname{diag}(p(t))$ , which means we have to approximate $R(t)$ at the same time. Intuitively, $\bar{R}(t)=k\operatorname{diag}(\hat{p}(t))$ would be a reasonable approximation. However, this would make (7) and (8) strongly coupled, for which it is very difficult to analyze the stability. Therefore, we use $\bar{R}(t)=\bar{k}\operatorname{diag}(p_{\text{KDE}}(t))$ , where $\bar{k}$ is computed as $\bar{k}=(\int[K(u)]^{2}du)/(nh^{d})$ according to (5). In this way, we can treat the approximation error as an external disturbance and use the concept of ISS to study its stability. We note that (5) should be understood to be accurate in the limit ( $n\to\infty$ ). When $n$ is finite, the estimate may be inaccurate. However, we shall point out that the estimation error of $\bar{k}$ is also included as part of the approximation error of $\bar{R}(t)$ .

By approximating $R(t)$ with $\bar{R}(t)=\bar{k}\operatorname{diag}(p_{\text{KDE}}(t))$ , the “suboptimal” density filter is correspondingly given by

\dot{\hat{p}}=A(t)\hat{p}+\bar{L}(t)(y-\hat{p}),\quad\hat{p}(t_{0})=p_{\text{KDE}}(t_{0})

(9)

where $\bar{L}(t)=\bar{P}(t)\bar{R}^{-1}(t)$ is the suboptimal Kalman gain and $\bar{P}(t)$ is a solution of the approximated Riccati equation

\dot{\bar{P}}=A\bar{P}+\bar{P}A^{*}-\bar{P}\bar{R}^{-1}\bar{P}.

(10)

The optimal density filter (7) essentially combines linearly all past density estimations $\{\hat{p}(\tau)\}_{t_{0}\leq\tau<t}$ and the most recent KDE measurement $p_{\text{KDE}}(t)$ to produce a density estimate $\hat{p}(t)$ with minimum estimation error covariance. The approximated density filter (9) is called “suboptimal” in the sense that the suboptimal gain $\bar{L}$ remains “close” to the optimal gain $L$ (which will be proved in Corollary 1).

IV-B Stability Analysis of the Suboptimal Density Filter

In this section, we study the stability of the suboptimal density filter (9) and the associated Riccati equation (10).

Let $P_{*}>0$ be the minimum (but unknown) covariance at $t=t_{0}$ . We denote by $\Pi(t)$ the solution of (8) with the initial condition $P_{*}$ , which represents the flow of minimum covariance and also generates the optimal gain. Let $P_{0}>0$ be a “guessed” initial condition. Denote by $\bar{\Pi}(t)$ the solution of (10) with the initial condition $P_{0}$ . We will show that $\bar{\Pi}(t)$ converges to and remains close to $\Pi(t)$ in the presence of approximation error on $R$ .

Define $\Gamma=\bar{\Pi}-\Pi$ . Using (8) and (10) we have

		$\displaystyle\dot{\Gamma}=A\Gamma+\Gamma A^{*}-\bar{\Pi}\bar{R}^{-1}\bar{\Pi}+\Pi R^{-1}\Pi,$		(11)
		$\displaystyle\Gamma(t_{0})=P_{0}-P_{*}.$		(11)

Define $\tilde{p}=\hat{p}-p$ . Then along $\bar{\Pi}(t)$ we have

\dot{\tilde{p}}=(A-\bar{\Pi}\bar{R}^{-1})\tilde{p}+\bar{\Pi}\bar{R}^{-1}w.

(12)

Our idea is to show that (11) is ISS with respect to the approximation error on $R$ (more precisely $\|\bar{R}^{-1}-R^{-1}\|$ , which equals $0$ if and only if $R=\bar{R}$ ). In this way, the suboptimal gain $\bar{L}$ remains close to $L$ when there is any approximation error, and converges to $0$ when the approximation error vanishes. We also show that the suboptimal density filter (9) is stable even though we use an approximation for $R$ .

Remark 2

To avoid confusions, we point out that if we use (8) to compute the gain for (7), then the estimation error covariance $\operatorname{Cov}[\tilde{p},\tilde{p}]$ along $\Pi(t)$ also satisfies (8), which is a property of Kalman filters [15]. However, if we use (10) to compute the gain for (9), then the covariance $\operatorname{Cov}[\tilde{p},\tilde{p}]$ along $\bar{\Pi}(t)$ , denoted by $Q(t)$ in this case, actually satisfies

\dot{Q}=AQ+QA^{*}-\bar{\Pi}(t)\bar{R}^{-1}R\bar{R}^{-1}\bar{\Pi}(t).

(13)

It would be desirable to prove that the solutions of (8) and (13) remain close, which is much harder because it involves three Riccati equations and will be left as our future work.

The following assumption is required for proving stability.

Assumption 1

Assume that $\|\Pi(t)\|$ and $\|\bar{\Pi}(t)\|$ are uniformly bounded, and that there exist positive constants $c_{1}$ and $c_{2}$ such that for all $t\geq t_{0}$ ,

0<c_{1}I\leq R^{-1}(t),\bar{R}^{-1}(t),\Pi^{-1}(t),\bar{\Pi}^{-1}(t)\leq c_{2}I.

(14)

Remark 3

The assumption for $R^{-1}(t)$ is an imposed requirement for $p(t)$ considering $R(t)=k\operatorname{diag}(p(t))$ , which roughly speaking, requires that $p(t)$ has positive upper bounds and lower bounds. The assumption for $\bar{R}^{-1}(t)$ can be easily satisfied because $\bar{R}(t)=\bar{k}\operatorname{diag}(p_{\text{KDE}}(t))$ and we construct $p_{\text{KDE}}(t)$ . For finite-dimensional systems, the assumptions for $\|\Pi(t)\|$ , $\|\bar{\Pi}(t)\|$ , $\Pi^{-1}(t)$ and $\bar{\Pi}^{-1}(t)$ follow the assumption for $R^{-1}(t)$ and $\bar{R}^{-1}(t)$ . This is because when $\Pi(t)$ and $\bar{\Pi}(t)$ are large, the solutions of (8) and (10) will be dominated by the negative second-order term and decay. Although intuitively correct, we find it much harder to prove for the infinite-dimensional case and thus leave it as our future work.

Stability results for (11) and (12) are given as follows.

Theorem 1

Under Assumption 1, the unforced part of (12), given by

\dot{\tilde{p}}=(A-\bar{\Pi}\bar{R}^{-1})\tilde{p},

(15)

is uniformly exponentially stable, and (11) is locally input-to-state stable (LISS) with respect to the approximation error in the form of $\|\bar{R}^{-1}-R^{-1}\|$ .

Proof:

To prove the first statement, consider a Lyapunov functional $V_{1}=\langle\bar{\Pi}^{-1}\tilde{p},\tilde{p}\rangle$ . We have

	$\displaystyle\dot{V}_{1}$	$\displaystyle=\big{\langle}\bar{\Pi}^{-1}\dot{\tilde{p}},\tilde{p}\big{\rangle}+\big{\langle}\bar{\Pi}^{-1}\tilde{p},\dot{\tilde{p}}\big{\rangle}-\big{\langle}\bar{\Pi}^{-1}\dot{\bar{\Pi}}\bar{\Pi}^{-1}\tilde{p},\tilde{p}\big{\rangle}$
		$\displaystyle=\big{\langle}\bar{\Pi}^{-1}(A-\bar{\Pi}\bar{R}^{-1})\tilde{p},\tilde{p}\big{\rangle}+\big{\langle}(A^{*}-\bar{R}^{-1}\bar{\Pi})\bar{\Pi}^{-1}\tilde{p},\tilde{p}\big{\rangle}$
		$\displaystyle\quad-\big{\langle}\bar{\Pi}^{-1}(A\bar{\Pi}+\bar{\Pi}A^{*}-\bar{\Pi}\bar{R}^{-1}\bar{\Pi})\bar{\Pi}^{-1}\tilde{p},\tilde{p}\big{\rangle}$
		$\displaystyle=-\langle\bar{R}^{-1}\tilde{p},\tilde{p}\rangle.$

In view of (14), we conclude that (15) is uniformly exponentially stable. Similarly, by considering a Lyapunov functional defined by $V_{2}=\langle\Pi^{-1}\tilde{p},\tilde{p}\rangle$ , one can show that the following system along $\Pi(t)$ is also uniformly exponentially stable:

\dot{\tilde{p}}=(A-\Pi R^{-1})\tilde{p}.

(16)

We also note that the inner products in $V_{1}$ and $V_{2}$ are equivalent because of the assumption (14). To prove the second statement, we rewrite (11) as

	$\displaystyle\dot{\Gamma}$	$\displaystyle=A\Gamma+\Gamma A^{*}-\bar{\Pi}\bar{R}^{-1}\Gamma-\bar{\Pi}\bar{R}^{-1}\Pi-\Gamma R^{-1}\Pi+\bar{\Pi}R^{-1}\Pi$		(17)
		$\displaystyle=(A-\bar{\Pi}\bar{R}^{-1})\Gamma+\Gamma(A^{*}-R^{-1}\Pi)-\bar{\Pi}(\bar{R}^{-1}-R^{-1})\Pi$		(17)

We note that (17) is essentially a linear equation. Now fix $q$ with $\|q\|=1$ . Since (15) and (16) are uniformly exponentially stable, and $\|\bar{\Pi}\|$ and $\|\Pi\|$ are assumed to be uniformly bounded, there exist constants $\lambda,c>0$ such that

	$\displaystyle\quad\\|\Gamma(t)q\\|$	$\displaystyle\leq e^{-\lambda(t-t_{0})}\\|\Gamma(t_{0})q\\|$
		$\displaystyle\quad+\int_{t_{0}}^{t}e^{-\lambda(t-\tau)}c\\|\bar{R}^{-1}(\tau)-R^{-1}(\tau)\\|\\|q\\|d\tau$

	$\displaystyle\leq e^{-\lambda(t-t_{0})}\\|\Gamma(t_{0})q\\|$
	$\displaystyle\quad+c\sup_{t_{0}\leq\tau\leq t}\\|\bar{R}^{-1}(\tau)-R^{-1}(\tau)\\|\int_{t_{0}}^{t}e^{-\lambda(t-s)}ds$
	$\displaystyle\leq e^{-\lambda(t-t_{0})}\\|\Gamma(t_{0})q\\|+\frac{c}{\lambda}\sup_{t_{0}\leq\tau\leq t}\\|\bar{R}^{-1}(\tau)-R^{-1}(\tau)\\|$
	$\displaystyle=:\beta(\\|\Gamma(t_{0})q\\|,t-t_{0})+\gamma(\sup_{t_{0}\leq\tau\leq t}\\|\bar{R}^{-1}(\tau)-R^{-1}(\tau)\\|),$

where $\beta\in\mathcal{KL}$ and $\gamma\in\mathcal{K}$ . Under the uniform boundedness assumption of $\Pi$ and $\bar{\Pi}$ , we conclude that (11) is LISS. ∎

We can conclude from Theorem 1 that the suboptimal gain $\bar{L}$ also remains close to the optimal gain $L$ , given as follows.

Corollary 1

Under Assumption 1, there exist functions $\beta_{1}\in\mathcal{KL}$ and $\gamma_{1}\in\mathcal{K}$ such that

\|\bar{L}-L\|\leq\beta_{1}(\|\Gamma(t_{0})\|,t-t_{0})+\gamma_{1}(\|\bar{R}^{-1}-R^{-1}\|).

Proof:

Observe that

	$\displaystyle\\|\bar{L}-L\\|$	$\displaystyle=\\|\bar{\Pi}\bar{R}^{-1}-\Pi R^{-1}\\|$
		$\displaystyle\leq\\|\bar{\Pi}\bar{R}^{-1}-\Pi\bar{R}^{-1}+\Pi\bar{R}^{-1}-\Pi R^{-1}\\|$
		$\displaystyle\leq\\|\bar{R}^{-1}\\|\\|\bar{\Pi}-\Pi\\|+\\|\Pi\\|\\|\bar{R}^{-1}-R^{-1}\\|$

In view of Definition 1, Theorem 1 and the uniform boundedness of $\|\Pi\|$ and $\|\bar{R}^{-1}\|$ , we obtain the desired result. ∎

V Simulation Studies

In this section, we study the performance of the proposed density filter. Consider a collection of $300$ agents given by

dX_{i}=\nabla\cdot\frac{D\nabla f(x)}{f(x)}dt+DdB_{t},\quad i=1,\dots,300,

(18)

where the states $\{X_{i}\}_{i=1}^{300}$ are restricted within $\Omega=[0,1]^{2}$ , $D=0.05$ and $f(x)$ is a continuous pdf over $\Omega$ to be specified. The initial conditions $\{X_{i}(t_{0})\}_{i=1}^{300}$ are drawn from the uniform distribution over $\Omega$ . Therefore, the ground truth density of the agents satisfies

\displaystyle\begin{split}&\partial_{t}p(x,t)=-\nabla\cdot\frac{Dp(x,t)\nabla f(x)}{f(x)}+\frac{1}{2}D^{2}\Delta p(x,t),\\ &p(\cdot,t_{0})=1,\end{split}

with a reflecting boundary condition (3). For illustration purpose, we let $f(x)$ be a Gaussian mixture of two components with the common covariance matrix $\text{diag}(0.015,0.015)$ but different time-varying means $[0.5+0.35\cos(0.2t),0.5+0.35\sin(0.2t)]^{\intercal}$ and $[0.5+0.35\cos(0.2t+\pi),0.5+0.35\sin(0.2t+\pi)]^{\intercal}$ . Under this design, the agents are nonlinear and time-varying and their states will converge to the two “spinning” Gaussian components.

We use the finite difference method [23] to numerically solve the density filter (9) and the operator Riccati equation (10). Specifically, $\Omega$ is equally divided into a $30\times 30$ grid. The pdfs $\hat{p}$ and $p_{\text{KDE}}$ are represented as $900\times 1$ vectors. The operators $A$ , $\bar{P}$ and $\bar{R}$ are represented as $900\times 900$ matrices. We set the initial conditions to be $\bar{P}(t_{0})=I$ and $\hat{p}(t_{0})=p_{\text{KDE}}(t_{0})$ . The time difference is $dt=0.1s$ and the bandwidth is $h=0.05$ . Note that $A$ is highly sparse and $\bar{R}$ is diagonal, so the computation is very fast in general.

Simulation results are given in Fig. 2. We see that the estimate by KDE is usually rugged, especially in the early stage when the true density is nearly uniform. (If we had known that it is nearly uniform, we would have chosen a much larger bandwidth $h$ for better smoothing.) In the late stage, the estimate by KDE has two major components since the samples are more concentrated. But it still has undesired modes due to outliers. The density filter, by taking advantage of the dynamics, quickly ”recognizes” the two underlying Gaussian components and gradually catches up with the evolution of the ground truth density. In Fig. 1, we compare the $L^{2}$ norms of estimation errors of the KDE and the density filter under different choices of bandwidth. It shows that the estimation error of the density filter quickly converges and the performance is robust to different choices of bandwidth.

Refer to caption — Figure 1: Estimation errors of the KDE (solid line) and the density filter (dashed line) with different choices of bandwidth.

VI Conclusion

We have presented a density filter for estimating the dynamic density of large-scale agent systems with known dynamics by combining KDE with infinite-dimensional Kalman filters. The density filter improved its estimation using the dynamics and the real-time states of the agents. It was scalable to the population of agents and was proved to be convergent. All results held even if the agents’ dynamics are nonlinear and time-varying. This algorithm can be used for many density-based optimization and control tasks of large-scale agent systems when density feedback information is required, and can be potentially extended to other dynamic density estimation problems when certain prior knowledge of its spatiotemporal evolution is available, such as population migration. Our future work includes addressing the remaining problems noted in Remark 2 and 3, decentralizing the density filter and integrating it into density feedback control for large-scale agent systems.

References

[1] G. Foderaro, P. Zhu, H. Wei, T. A. Wettergren, and S. Ferrari, “Distributed optimal control of sensor networks for dynamic target tracking,” IEEE Transactions on Control of Network Systems, vol. 5, no. 1, pp. 142–153, 2016.
[2] B. W. Silverman, Density Estimation for Statistics and Data Analysis. CRC Press, 1986, vol. 26.
[3] G. Foderaro, S. Ferrari, and M. Zavlanos, “A decentralized kernel density estimation approach to distributed robot path planning,” in Proceedings of the Neural Information Processing Systems Conference, Lake Tahoe, NV, 2012.
[4] V. Krishnan and S. Martínez, “Distributed control for spatial self-organization of multi-agent swarms,” SIAM Journal on Control and Optimization, vol. 56, no. 5, pp. 3642–3667, 2018.
[5] M. H. de Badyn, U. Eren, B. Açikmeşe, and M. Mesbahi, “Optimal mass transport and kernel density estimation for state-dependent networked dynamic systems,” in 2018 IEEE Conference on Decision and Control (CDC). IEEE, 2018, pp. 1225–1230.
[6] S. R. Sain, “Multivariate locally adaptive density estimation,” Computational Statistics & Data Analysis, vol. 39, no. 2, pp. 165–186, 2002.
[7] S. J. Sheather and M. C. Jones, “A reliable data-based bandwidth selection method for kernel density estimation,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 53, no. 3, pp. 683–690, 1991.
[8] A. Qahtan, S. Wang, and X. Zhang, “Kde-track: An efficient dynamic density estimator for data streams,” IEEE Transactions on Knowledge and Data Engineering, vol. 29, no. 3, pp. 642–655, 2016.
[9] S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” Proceedings of the IEEE, vol. 92, no. 3, pp. 401–422, 2004.
[10] I. Arasaratnam and S. Haykin, “Cubature kalman filters,” IEEE Transactions on automatic control, vol. 54, no. 6, pp. 1254–1269, 2009.
[11] Z. Chen et al., “Bayesian filtering: From kalman filters to particle filters, and beyond,” Statistics, vol. 182, no. 1, pp. 1–69, 2003.
[12] R. Olfati-Saber, “Distributed kalman filtering for sensor networks,” in 2007 46th IEEE Conference on Decision and Control. IEEE, 2007, pp. 5492–5498.
[13] S. Bandyopadhyay and S.-J. Chung, “Distributed estimation using bayesian consensus filtering,” in 2014 American control conference. IEEE, 2014, pp. 634–641.
[14] G. Battistelli and L. Chisci, “Kullback–leibler average, consensus on probability densities, and distributed state estimation with guaranteed stability,” Automatica, vol. 50, no. 3, pp. 707–718, 2014.
[15] R. E. Kalman and R. S. Bucy, “New results in linear filtering and prediction theory,” Journal of basic engineering, vol. 83, no. 1, pp. 95–108, 1961.
[16] P. L. Falb, “Infinite-dimensional filtering: The kalman-bucy filter in hilbert space,” Information and Control, vol. 11, no. 1-2, pp. 102–137, 1967.
[17] A. Bensoussan, Filtrage optimal des systèmes linéaires. Dunod, 1971.
[18] G. Da Prato and J. Zabczyk, Stochastic equations in infinite dimensions. Cambridge university press, 2014.
[19] E. D. Sontag and Y. Wang, “On characterizations of the input-to-state stability property,” Systems & Control Letters, vol. 24, no. 5, pp. 351–359, 1995.
[20] S. Dashkovskiy and A. Mironchenko, “Input-to-state stability of infinite-dimensional control systems,” Mathematics of Control, Signals, and Systems, vol. 25, no. 1, pp. 1–35, 2013.
[21] G. A. Pavliotis, Stochastic processes and applications: diffusion processes, the Fokker-Planck and Langevin equations. Springer, 2014, vol. 60.
[22] T. Cacoullos, “Estimation of a multivariate density,” Annals of the Institute of Statistical Mathematics, vol. 18, no. 1, pp. 179–189, 1966.
[23] J. Chang and G. Cooper, “A practical difference scheme for fokker-planck equations,” Journal of Computational Physics, vol. 6, no. 1, pp. 1–16, 1970.

	$\displaystyle\leq e^{-\lambda(t-t_{0})}\\|\Gamma(t_{0})q\\|$
	$\displaystyle\quad+c\sup_{t_{0}\leq\tau\leq t}\\|\bar{R}^{-1}(\tau)-R^{-1}(\tau)\\|\int_{t_{0}}^{t}e^{-\lambda(t-s)}ds$
	$\displaystyle\leq e^{-\lambda(t-t_{0})}\\|\Gamma(t_{0})q\\|+\frac{c}{\lambda}\sup_{t_{0}\leq\tau\leq t}\\|\bar{R}^{-1}(\tau)-R^{-1}(\tau)\\|$
	$\displaystyle=:\beta(\\|\Gamma(t_{0})q\\|,t-t_{0})+\gamma(\sup_{t_{0}\leq\tau\leq t}\\|\bar{R}^{-1}(\tau)-R^{-1}(\tau)\\|),$

	$\displaystyle\\|\bar{L}-L\\|$	$\displaystyle=\\|\bar{\Pi}\bar{R}^{-1}-\Pi R^{-1}\\|$
		$\displaystyle\leq\\|\bar{\Pi}\bar{R}^{-1}-\Pi\bar{R}^{-1}+\Pi\bar{R}^{-1}-\Pi R^{-1}\\|$
		$\displaystyle\leq\\|\bar{R}^{-1}\\|\\|\bar{\Pi}-\Pi\\|+\\|\Pi\\|\\|\bar{R}^{-1}-R^{-1}\\|$