Distributed Kalman-filtering: Distributed optimization viewpoint

TL;DR

This paper reformulates distributed Kalman filtering as a consensus optimization problem and introduces a new algorithm based on distributed dual ascent, enabling sensor networks to perform filtering without a central fusion center.

Contribution

It presents a novel reformulation of distributed Kalman filtering as a consensus optimization problem and develops a new dual ascent based algorithm for it.

Findings

The new algorithm performs effectively in numerical experiments.

Distributed optimization techniques can be applied to Kalman filtering.

The approach eliminates the need for a central fusion center.

Abstract

We consider the Kalman-filtering problem with multiple sensors which are connected through a communication network. If all measurements are delivered to one place called fusion center and processed together, we call the process centralized Kalman-filtering (CKF). When there is no fusion center, each sensor can also solve the problem by using local measurements and exchanging information with its neighboring sensors, which is called distributed Kalman-filtering (DKF). Noting that CKF problem is a maximum likelihood estimation problem, which is a quadratic optimization problem, we reformulate DKF problem as a consensus optimization problem, resulting in that DKF problem can be solved by many existing distributed optimization algorithms. A new DKF algorithm employing the distributed dual ascent method is provided and its performance is evaluated through numerical experiments.

Figures

Equations

$$x_{k+1}$$

$$y_k$$

$${\hat{x}}_{k|k-1}$$

$$P_{k|k-1}$$

$$=F\mathbb{E}\{e_{k-1}{e}_{k-1}^{\top }\}{F}^{\top }+\mathbb{E}\{w_k{w}_k^{\top }\}$$

$$=F{P}_{k-1}{F}^{\top }+Q$$

$$\mathfrak{L}({\xi }_c)=\frac{1}{\sqrt{(2\pi {)}^{(m+n)}|S_k|}}{e}^{-\frac{1}{2}(z_k-{\bar{H}}_c{\xi }_c{)}^{\top }{S}_k^{-1}(z_k-{\bar{H}}_c{\xi }_c)}$$

$${\hat{x}}_k:=\underset{{\xi }_c}{\mathrm{argmax}}(\mathfrak{L}({\xi }_c)).$$

$$\begin{array}{rl}{\hat{x}}_k& =\underset{{\xi }_c}{\mathrm{argmin}}(f_c({\xi }_c))\\ & ={\hat{x}}_{k|k-1}+K(y_k-H{\hat{x}}_{k|k-1})\end{array}$$

$$P_k$$

$$=P_{k|k-1}-(H^{\top }{R}^{-1}H+P_{k|k-1}^{-1}{)}^{-1}{H}^{\top }{R}^{-1}H{P}_{k|k-1}.$$

$${\hat{x}}_{i,k|k-1}$$

$$P_{i,k|k-1}$$

$$\sum _{i=1}^N{f}_i({\xi }_c)=\sum _{i=1}^N\frac{1}{2}({\bar{z}}_{i,k}-{\bar{H}}_i{\xi }_c{)}^{\top }{\bar{S}}_{i,k}^{-1}({\bar{z}}_{i,k}-{\bar{H}}_i{\xi }_c)$$

$$\sum _{i=1}^N{f}_i({\xi }_i)$$

$${\xi }_1=\cdots ={\xi }_N.$$

$$\mathcal{L}(\xi ,\lambda )$$

$${\mathcal{L}}_i({\xi }_i,{\lambda }_i)$$

$${\nabla }_{\xi }\mathcal{L}(\xi ,\lambda )$$

$${\nabla }_{\lambda }\mathcal{L}(\xi ,\lambda )$$

$$\left[\begin{array}{cc}-{\bar{H}}^{\top }{\bar{S}}_k^{-1}\bar{H}& -\bar{L}\\ \bar{L}& 0\end{array}\right]\left[\begin{array}{c}{\xi }^{\ast }\\ {\lambda }^{\ast }\end{array}\right]=\left[\begin{array}{c}-{\bar{H}}^{\top }{\bar{S}}_k^{-1}{\bar{z}}_k\\ 0\end{array}\right]$$

$$(1_N^{\top }\otimes I_n){\bar{H}}^{\top }{\bar{S}}_k^{-1}\bar{H}{\xi }^{\ast }=(1_N^{\top }\otimes I_n){\bar{H}}^{\top }{\bar{S}}_k^{-1}{\bar{z}}_k.$$

$$(1_N^{\top }\otimes I_n){\bar{H}}^{\top }{\bar{S}}_k^{-1}\bar{H}(1_N\otimes I_n){\xi }^{†}=(1_N^{\top }\otimes I_n){\bar{H}}^{\top }{\bar{S}}_k^{-1}{\bar{z}}_k.$$

$$\displaystyle\Big{(}P^{-1}_{k|k-1}$$

$$=P_{k|k-1}^{-1}{\hat{x}}_{k|k-1}+\sum _{i=1}^N{H}_i^{\top }{R}_i^{-1}{y}_{i,k}.$$

$${\xi }^{†}={\hat{x}}_{k|k-1}+K_k(y_k-H{\hat{x}}_{k|k-1})$$

$$(L\otimes I_n){\lambda }^{\ast }={\bar{H}}^{\top }{\bar{S}}_k^{-1}({\bar{z}}_k-\bar{H}(1_N\otimes I_n){\xi }^{†})$$

$$\left(\left[\begin{array}{cc}0& \\ & \bar{\Lambda }\end{array}\right]\otimes I_n\right)\left(\left[\begin{array}{c}{U}_1^{\top }\\ {\bar{U}}^{\top }\end{array}\right]\otimes I_n\right){\lambda }^{\ast }=\left(\left[\begin{array}{c}{U}_1^{\top }\\ {\bar{U}}^{\top }\end{array}\right]\otimes I_n\right)b$$

$$\sum _{i=1}^N{h}_i({\eta }_i)$$

$${\eta }_1=\cdots ={\eta }_N$$

$$h_i({\eta }_i)=\frac{1}{2}$$

$$\displaystyle\hskip 25.6073pt+y_{i}^{\top}R_{i}^{-1}y_{i}+\frac{1}{N}\tau_{i,k|k-1}^{\top}\Omega_{i,k|k-1}^{-1}\tau_{i,k|k-1}\Big{)}$$

$${\mathcal{L}}_{\eta }(\eta ,\lambda )$$

$$\left[\begin{array}{cc}-({\bar{H}}^{\top }{\tilde{S}}_k^{-1}\bar{H}{)}^{-1}& -\bar{L}\\ \bar{L}& 0\end{array}\right]\left[\begin{array}{c}{\eta }^{\ast }\\ {\nu }^{\ast }\end{array}\right]=\left[\begin{array}{c}-{\bar{H}}^{\top }{\tilde{S}}_k^{-1}{\tilde{z}}_k\\ 0\end{array}\right]$$

$${\nabla }_{{\eta }_i}{h}_i({\eta }_i)={\Phi }_i^{-1}{\eta }_i-{\Phi }_i^{-1}(H_i^{\top }{R}_i^{-1}{y}_i+\frac{1}{N}{\tau }_{i,k|k-1})$$

$${\xi }_{l+1}$$

Full text

I Introduction

It goes without saying that the Kalman-filter, an optimal state estimator for dynamic systems, has had a huge impact on various fields such as engineering, science, economics, etc. [1, 2, 3, 4]. Basically, the filter predicts the expectation of the system state and its covariance based on the dynamic model and the statistical information on the model uncertainty or process noise, and then correct them using new measurement, sensor model, and the information on measurement noise. When multiple sensors possibly different types are available, we can just combine the sensor models to process the measurements altogether.

Thanks to the rapid development of sensor devices and communication technology, we are now able to monitor large scale systems or environments such as traffic network, plants, forest, sea, etc. In those systems, sensors are geometrically distributed, may have different types, and usually not synchronized. To process the measurements, the basic idea would be to deliver all the data to one place, usually called fusion center, and do the correction step as in the case of multiple sensors. This is called the centralized Kalman-filtering (CKF). As expected, CKF requires a powerful computing device to handle a large number of measurements and sensor models, is exposed to a single point of failure, and is difficult to scale up. In order to overcome these drawbacks, researchers developed the distributed Kalam-filtering (DKF) in which each sensor in the network solves the problem by using local measurements and communicating with its neighbors. Compared with CKF, DKF has advantageous in terms of the scalability, robustness to component loss, computational cost, and thus the literature on this topic is expanding rapidly [5, 6, 7, 8, 9, 10, 11, 12]. For more details on DKF, see the survey [13] and references therein.

Some relevant results are summarized as follows. In [5], the author proposed scalable distributed Kalman-Bucy filtering algorithms in which each node only communicates with its neighbors. An algorithm with average consensus filters using the internal models of signals being exchanged is proposed in [7]. It is noted that the algorithm works in a single-time scale. In the work [11], the authors proposed a continuous-time algorithm that makes each norm of all local error covariance matrices be bounded, thus overcomes a major drawback of [5]. In [10], an algorithm with a high gain coupling term in the error covariance matrix is introduced and it is shown that the local error covariance matrix approximately converges to that of the steady-state centralized Kalman-filter. An in-depth discussion on distributed Kalman-filtering problem has been provided in [14, 15], and the algorithms that exchange the measurements themselves, or exchange certain signals instead of the measurements are proposed, respectively.

Although each of the existing algorithms has own novel ideas and advantages, to the best of the authors’ knowledge, we do not have a unified viewpoint for DKF problem. Motivated by this, it is the aim of this paper to provide a framework for the problem from the perspective of distributed optimization.

We start by observing that the correction step of Kalman-filtering is basically an optimization problem [2, 3, 4], and then formulate DKF problem as a consensus optimization problem, which provides a fresh look at the problem. This results in that DKF problem can be solved by many existing distributed optimization algorithms [16, 17, 18, 19, 20], expecting various DKF algorithms to be derived. As an instance, a new DKF algorithm employing the dual ascent method [20], one of the basic algorithms for distributed optimization problems, is provided in this paper.

This paper is organized as follows. In Section II, we recall CKF problem from the optimization perspective, and connects DKF problem to a distributed optimization problem. A new DKF algorithm based on dual ascent method is proposed in Section III, and numerical experiments evaluating the proposed algorithm is conducted in Section IV.

Notation: For matrices $A_{1}$, …, $A_{n}$, $\operatorname{\text{diag}}(A_{1},\dots,A_{n})$ denotes the block diagonal matrix composed of $A_{1}$ to $A_{n}$. For scalars $a_{1}$,…, $a_{n}$, $[a_{1};\dots;a_{n}]:=[a_{1}^{\top},\dots,a_{n}^{\top}]^{\top}$, and $[A_{1};\dots;A_{n}]$ with matrices $A_{i}$’s is defined similarly. $1_{n}\in\mathbb{R}^{n}$ denotes the vector whose components are all 1, and $I_{n}$ is the identity matrix whose dimension is $n\times n$. The maximum and minimum eigenvalue of a matrix $A$ are denoted by $\sigma_{\max}(A)$ and $\sigma_{\min}(A)$, respectively. For a random variable $x$, $x\sim\mathsf{N}(\mu,\sigma^{2})$ denotes $x$ is normally distributed with the mean $\mu$ and the variance $\sigma^{2}$, and $\mathbb{E}\{x\}$ denotes the expected value of a random variable $x$, i.e., $\mathbb{E}\{x\}=\mu$. The half vectorization of a symmetric matrix $M\in\mathbb{R}^{n\times n}$ is denoted by $\text{vec}_{h}({M})\in\mathbb{R}^{n(n+1)/2}$, whose elements are filled in Column-major order. $i.e.,\text{vec}_{h}({M}):=[M_{1,1};\dots;M_{1,n};M_{2,2};\dots;M_{2,n};\dots;$ $M_{n-1,n-1};M_{n-1,n};M_{n,n}]$ where $M_{i,j}$ is $i,j$ element of $M$, and $\text{vec}_{h}^{-1}({\cdot})$ denotes the inverse function of $\text{vec}_{h}({\cdot})$, $i.e.,\text{vec}_{h}^{-1}({\text{vec}_{h}({M})})=M$. For a function $f(x,y):\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$, $\nabla_{x}f(x,y)$ denotes the gradient vector $\frac{\partial f(x,y)}{\partial x}=[\frac{\partial f(x,y)}{\partial x_{1}};\dots;\frac{\partial f(x,y)}{\partial x_{n}}]$.

Graph theory: For a network consisting of $N$ nodes, the communication among nodes is modeled by a graph $\mathcal{G}$. Let ${\mathcal{A}}=[a_{ij}]\in{\mathbb{R}}^{N\times N}$ be an adjacency matrix associated to ${\mathcal{G}}$ where $a_{ij}$ is a weight of an edge between nodes $i$ and $j$. If node $i$ communicates to node $j$ then, $a_{ij}>0$, or if not $a_{ij}=0$. Assume there is not self edge, i.e., $a_{ii}=0$. The Laplacian matrix associated to the graph $\mathcal{G}$, denoted by $L$ is a $N\times N$ matrix such that $l_{ij,i\neq j}=-a_{ij}$, and $l_{ii}=\sum_{j=1}^{N}a_{ij}$. ${\mathcal{N}}_{i}$ is a set of nodes communicating with node $i$, i.e., ${\mathcal{N}}_{i}=\{j|a_{ij}>0\}$.

II Distributed Kalman-filtering and Its Connection to Consensus Optimization

In this section, we recall CKF problem in terms of optimization, which is the maximum likelihood estimation[2], and establish a connection between DFK and distributed optimization.

Consider a discrete-time linear system with $N$ sensors described by

[TABLE]

where $x_{k}\in{\mathbb{R}}^{n}$ is the state vector of the dynamic system, $y_{k}:=[y_{1,k};\dots;y_{N,k}]\in\mathbb{R}^{m}$ is the output vector, and $y_{i,k}\in\mathbb{R}^{m_{i}}$ is the output associated to sensor $i$. $m_{i}$’s satisfy $\sum^{N}_{i=1}m_{i}=m$. $F$ is the system matrix and $H$ is the output matrix consisting of $H_{i}\in\mathbb{R}^{m_{i}\times n}$ which is the output matrix associated to sensor $i$. $w_{k}\in{\mathbb{R}}^{n}$ with $w_{k}\sim\mathsf{N}(0,Q)$ is the process noise, $v_{i,k}\sim\mathsf{N}(0,R_{i})$ is the measurement noise on sensor $i$, and $v_{k}:=[v_{1,k};\dots;v_{N,k}]\in{\mathbb{R}}^{m}$ with $v_{k}\sim\mathsf{N}(0,\operatorname{\text{diag}}(R_{1},\dots,R_{N}))$. Assume that the pair $(F,H)$ is observable, and each $v_{i,k}$ is uncorrelated to $v_{{j},k}$ for $j\neq i$.

II-A Centralized Kalman-filtering problem from the optimization perspective

If all the measurements from $N$ sensors are collected and processed altogether, the problem can be seen as the one with a imaginary sensor that measures $y_{k}$ with complete knowledge on $H$, thus called centralized Kalman-filtering.The filtering consists of two steps, prediction and correction. In the prediction step, the predicted estimate $\hat{x}_{k|k-1}$ and error covariance matrix $P_{k|k-1}$ are obtained based on the previous estimate, error covariance matrix, and the system dynamics. The update rules are given by

[TABLE]

where $\hat{x}_{k-1}$ and $P_{k-1}$ are estimate and error covariance matrix in previous time, respectively, and $e_{k|k-1}:=x_{k}-\hat{x}_{k|k-1}$, $e_{k}:=x_{k}-\hat{x}_{k}$. Assume that $P_{k}$ is initialized as a positive definite matrix ($P_{0}>0$, usually set as $Q$).

In the correction step, the predicted estimate and the error covariance matrix are updated based on the current measurements containing the measurement noise. The correction step can be regarded as a process to find the optimal parameter (estimate) from the predicted estimate $\hat{x}_{k|k-1}$, error covariance $P_{k|k-1}$, and the observation $y_{k}$. In fact, it is known that this step is an optimization problem (maximum likelihood estimation, MLE[2]) and we recall the details below.

Let ${z}_{k}=[y_{k};\hat{x}_{k|k-1}]\in\mathbb{R}^{m+n}$ and $\bar{H}_{c}=[H;I_{n}]\in\mathbb{R}^{(m+n)\times n}$. Then, ${z}_{k}\sim\mathsf{N}({\bar{H}_{c}}x_{k},S_{k})$ where $S_{k}=\text{diag}\{R,P_{k|k-1}\}$. For the random variable ${z}_{k}$, the likelihood function is given by

[TABLE]

where the right-hand side is nothing but the probability density function of $z_{k}$ with the free variable $\xi_{c}\in\mathbb{R}^{n}$.

Now, the maximum likelihood estimate $\hat{x}_{k}$ is defined as

[TABLE]

Since $\mathfrak{L}(\xi_{c})$ is a monotonically decreasing function with respect to $f_{\text{c}}(\xi_{c}):=\frac{1}{2}(z_{k}-\bar{H}_{c}\xi_{c})^{\top}S^{-1}(z_{k}-\bar{H}_{c}\xi_{c})$, $\hat{x}_{k}$ can also be obtained by

[TABLE]

where $K=(H^{\top}R^{-1}H+P^{-1}_{k|k-1})^{-1}H^{\top}R^{-1}$. With the matrix inversion lemma, the Kalman-gain $K$ can be written as $K=P_{k|k-1}H^{\top}(HP_{k|k-1}H^{\top}+R)^{-1}$, which appears in the standard Kalman-filtering.

On the other hand, by the definition of $P_{k}:={\mathbb{E}}\{(\hat{x}_{k}-x_{k})(\hat{x}_{k}-x_{k})^{\top}\}$, the update rule of the error covariance matrix $P_{k}$ of CKF is given by

[TABLE]

For more details, see [4, 2, 3].

II-B Derivation of distributed Kalman-filtering problem

Now, we consider a sensor network which consists of $N$ sensors and suppose that each sensor runs an estimator without the fusion center. Each estimator in the network tries to find the optimal estimate by processing the local measurement and exchanging information with its neighbors through communication network. The communication network among estimators is modeled by a graph $\mathcal{G}$ and the Laplacian matrix associated with $\mathcal{G}$ is denoted by $L\in\mathbb{R}^{N\times N}$. Under the setting (1), estimator $i$ measures only the local measurement $y_{k,i}$, and the parameters $H_{i}$ and $R_{i}$ are kept private to estimator $i$. It is noted that the pair $(F,H_{i})$ is not necessarily observable. We assume that the graph is connected and undirected i.e., $L=L^{\top}$, and $F$ and $Q$ are open to all estimators.

Similar to CKF, DKF has two steps, local prediction and distributed correction. In the local prediction step, each estimator predicts

[TABLE]

where $\hat{x}_{i,k|k-1}$ and $P_{i,k|k-1}$ are local estimates of $\hat{x}_{k|k-1}$ and $P_{k|k-1}$, respectively, that estimator $i$ holds.

In the distributed correction step, each estimator solves the maximum likelihood estimation in a distributed manner. The objective function of CKF $f_{c}(\xi)$ can be rewritten as

[TABLE]

where $\bar{z}_{i,k}=[y_{i,k};\hat{x}_{i,k|k-1}]$, $\bar{H}_{i}=[H_{i};I_{n}]$, $\bar{S}_{i,k}=\operatorname{\text{diag}}(R_{i},NP_{i,k|k-1})$. We assume that $\hat{x}_{i,k|k-1}=\hat{x}_{k|k-1}$ and $P_{i,k|k-1}=P_{k|k-1}$. This makes sense when the each sensor reached a consensus on $\hat{x}_{i,k-1}$ and $P_{i,k-1}$ in the previous correction step.

Assuming that each estimator holds its own optimization variable $\xi_{i}\in\mathbb{R}^{n}$ for $\xi_{c}$, DKF problem is written as the following consensus optimization problem.

[TABLE]

If there exists a distributed algorithm that finds a minimizer $(\xi_{1}^{*},\dots,\xi_{N}^{*})$, we say that the algorithm solves DKF problem.

Since the kernel of Laplacian $L$ is $\text{span}\{1_{N}\}$, the constraints (4b) can be written with as $(L\otimes I_{n})\xi=0$ where $\xi=[\xi_{1};\dots;\xi_{N}]$. To proceed, we define the Lagrangian to solve the problem (4) as

[TABLE]

where $\lambda\in\mathbb{R}^{Nn}$ is the Lagrange multipliers (dual variable) associated with (4b) and $\bar{L}=(L\otimes I_{n})$. We decompose the Lagrangian into local ones defined by

[TABLE]

For the Lagrangian (5), the partial derivatives over $\xi$ and $\lambda$ are given by

[TABLE]

where $\bar{z}_{k}=[\bar{z}_{1,k};\dots;\bar{z}_{N,k}]$, $\bar{H}=\operatorname{\text{diag}}(\bar{H}_{1},\dots,\bar{H}_{N})$ and $\bar{S}_{k}=\operatorname{\text{diag}}(\bar{S}_{1,k},\dots,\bar{S}_{N,k})$. Then, the optimality condition for ($\xi^{*}$, $\lambda^{*}$) becomes the following saddle point equation (KKT conditions), namely

[TABLE]

where $\xi^{*}:=[\xi_{1}^{*};\dots;\xi_{N}^{*}]$ and $\lambda^{*}:=[\lambda_{1}^{*};\dots;\lambda_{N}^{*}]$.

Lemma 1

The solutions to DKF problem are parameterized as $(\xi^{*},\lambda^{*})=((1_{N}\otimes I_{n})\xi^{\dagger},(1_{N}\otimes I_{n})\tilde{\lambda}+\bar{\lambda})$ where $\xi^{\dagger}\in\mathbb{R}^{n}$ and $\bar{\lambda}\in\mathbb{R}^{Nn}$ are unique vectors and $\tilde{\lambda}\in\mathbb{R}^{n}$ is an arbitrary vector. If $(\xi^{*},\lambda^{*})$ is an optimal solution to DKF problem, then $\xi_{i}^{*}$ is the optimal solution to CFK problem. $\diamond$

Proof:

By multiplying $1_{N}^{\top}\otimes I_{n}$ to the dual feasibility equation in (7), one can obtain

[TABLE]

The primal feasibility equation in (7) implies that $\xi^{*}=(1_{N}\otimes I_{n})\xi^{\dagger}$, hence (8) becomes

[TABLE]

From $\bar{H}^{\top}\bar{S}_{k}^{-1}\bar{H}=\operatorname{\text{diag}}(\frac{1}{N}P^{-1}_{k|k-1}+H_{1}^{\top}R_{1}^{-1}H_{1},\dots,$ $\frac{1}{N}P^{-1}_{k|k-1}+H_{N}^{\top}R_{N}^{-1}H_{N})$, one has

[TABLE]

Since $\sum_{i=1}^{N}H_{i}^{\top}R_{i}^{-1}y_{i,k}=H^{\top}R^{-1}y_{k}$ and $\sum_{i=1}^{N}H_{i}^{\top}R_{i}^{-1}H_{i}=H^{\top}R^{-1}H$, it follows that

[TABLE]

where $K_{k}=(P^{-1}_{k|k-1}+H^{\top}R^{-1}H)^{-1}H^{\top}R^{-1}$ and by the matrix inversion lemma, we have $K_{k}=P_{k|k-1}H^{\top}(R+HP_{k|k-1}H^{\top})^{-1}$. From the fact that the right-hand side of above equation is the same with the update rule (2) of CKF, it follows that $\xi_{i}^{*}=\xi^{\dagger}$ is the optimal estimate of CKF $\hat{x}_{k}$.

On the other hand, one can observe that the optimal dual variable $\lambda^{*}$ is not unique since the dual feasibility equation

[TABLE]

is singular. To find $\eta^{*}$, consider the orthonormal matrix $U=[U_{1}~{}\bar{U}]$ such that $LU=U\Lambda$ where $U_{1}=\frac{1}{\sqrt{N}}1_{N}$, $\bar{U}$ consists of the eigenvectors associated with the non-zero eigenvalues of $L$, denoted by $\sigma_{2},$…$,\sigma_{N}$, and $\Lambda=\operatorname{\text{diag}}({0,\sigma_{2},\dots,\sigma_{N}})$. Left multiplying $U^{\top}\otimes I_{n}$ to the equation (10) yields

[TABLE]

where $\bar{\Lambda}=\operatorname{\text{diag}}(\sigma_{2},\dots,\sigma_{N})$ and $b=\bar{H}^{\top}\bar{S}_{k}^{-1}(\bar{z}_{k}-\bar{H}(1_{N}\otimes I_{n})\xi^{\dagger})$. Hence, the optimal dual variable $\lambda^{*}$ becomes $\lambda^{*}=(U\otimes I_{n})\begin{bmatrix}\tilde{\lambda}^{*};(\bar{\Lambda}^{-1}\bar{U}^{\top}\otimes I_{n})b\end{bmatrix}$ where $\tilde{\lambda}^{*}\in\mathbb{R}^{n}$ is an arbitrary vector. This completes the proof. ∎

II-C Information form of DKF problem

It is well known that the dual of the Kalman-filter is the Information filter which uses the canonical parameterization to represent the normal (Gaussian) distribution [4]. With the canonical parameterization, DKF problem (4) can also be written in information form.

Let $\eta_{i}=(H_{i}^{\top}R^{-1}_{i}H_{i}+\frac{1}{N}\Omega_{i,k|k-1})\xi_{i}$, $\Omega_{i,k|k-1}=P_{i,k|k-1}^{-1}$ and $\tau_{i,k|k-1}=P_{i,k|k-1}^{-1}\hat{x}_{i,k|k-1}$ which are the local decision variable for the information vector of the estimator $i$, the locally predicted information matrix and information vector, respectively. With these transformations, we rewrite the problem (4) as

[TABLE]

where

[TABLE]

and $\Phi_{i}=H_{i}^{\top}R^{-1}_{i}H_{i}+\frac{1}{N}\Omega_{i,k|k-1}$. For the distributed problem (11), the Lagrangian is given by

[TABLE]

where $\eta:=[\eta_{1};\dots;\eta_{N}]$ and $\nu$ is the Lagrange multipliers. The associated saddle point equation becomes

[TABLE]

where $\tilde{z}_{k}=[\tilde{z}_{1,k};\dots;\tilde{z}_{N,k}]$, $\tilde{S}_{k}=\operatorname{\text{diag}}(\tilde{S}_{1,k},\dots,\tilde{S}_{N,k})$, $\tilde{z}_{i,k}=[y_{i,k};\tau_{i,k|k-1}]$ and $\tilde{S}_{i,k}=\operatorname{\text{diag}}(R_{i},N\Omega_{i,k|k-1}^{-1})$.

II-D Interpretations of existing DKF algorithm from the optimization perspective

One of the recent DKF algorithms, Consensus on Information (CI) [14, 15] can be interpreted in the provided framework. CI consists of three steps, prediction, local correction, and consensus. In the prediction step, each estimator predicts the estimate based on the system dynamics and previous estimate similar to the standard information filter algorithm. Each estimator also updates the estimate with local measurements and output matrix in the local correction step. After that, the estimators find the agreed estimate by averaging the local estimates in the consensus step.

In the provided framework, CI can be viewed as the algorithm which solves the problem (11) through the two steps, the local correction step and the consensus step. In the former step, each of estimators finds the local minimizer (estimate) of the local objective function $h_{i}(\cdot)$. Since the partial derivative of $h_{i}(\eta_{i})$ becomes

[TABLE]

and the local minimizer $\eta_{i}^{*}$ can be obtained by $\eta^{*}_{i}=H_{i}^{\top}R_{i}^{-1}y_{i}+\frac{1}{N}\tau_{i,k|k-1}$, which is the local update rule of CI111In the CI, the scalar $\frac{1}{N}$ is neglected [14].. The local minimizer, however, can be different among estimators, since it minimizes only the local objective function $h_{i}(\cdot)$, which violates the constraint (11b).

The consensus step of CI performs a role to find an agreed (average) value of the local estimates, using the doubly stochastic matrix, and the results of the consensus step satisfy the constraint (11b). The agreed estimate, however, may not be the global minimizer of (11), which means that the consensus step cannot guarantee the convergence of the estimates to that of CKF.

III A Solution to DKF Problem

One can observe that (5) is strictly convex, differentiable, and the local objective function $f_{i}(\cdot)$ is a quadratic function, hence strong duality holds. In addition, from the fact $\bar{H}^{\top}\bar{S}_{k}^{-1}\bar{H}$ is a nonsingular and block diagonal matrix, the optimal conditions (7) are already in a distributed form. This implies that the minimizer $\xi^{*}$ can be obtained in a distributed manner as long as $\lambda^{*}$ is given, i.e., $\xi_{i}^{*}=(\bar{H}_{i}^{\top}\bar{S}_{i,k}^{-1}\bar{H}_{i})^{-1}(\bar{H}_{i}^{\top}\bar{S}_{i,k}^{-1}\bar{z}_{i,k}-\sum_{j\in\mathcal{N}_{i}}a_{ij}(\lambda^{*}_{i}-\lambda^{*}_{j}))$.

Based on the above discussion, we see that one possible algorithm solving (4), guaranteeing the asymptotic convergence to the global minimizer $\xi^{*}$, is the dual ascent method [16, 20] which is given by

[TABLE]

where $\alpha_{\lambda}>0$ is a step size. The update rule (12) can be written locally as

[TABLE]

where $K_{i,k}=(H_{i}^{\top}R_{i}^{{-1}}H_{i}+\frac{1}{N}P_{i,k|k-1}^{{-1}})^{-1}H_{i}^{\top}R_{i}^{-1}$, $\psi_{i,l}=(H_{i}^{\top}R_{i}^{-1}H_{i}+\frac{1}{N}P_{i,k|k-1}^{-1})^{-1}\sum_{j\in\mathcal{N}_{i}}a_{ij}(\lambda_{i,l}-\lambda_{j,l})$, and $l$ is the iteration index to find the minimizer.

Regarding the convergence of the update rule (13), we have the following result.

Lemma 2

Assume that the network $\mathcal{G}$ is undirected and connected. Then, the sequence $\{\xi_{i,l}\}$ generated by the dual ascent method (13) converges to $\hat{x}_{k}$ of CKF problem (2), as $l$ goes to infinity, provided that the step size $\alpha_{\lambda}>0$ is chosen such that

[TABLE]

where $\sigma_{N}$ is the maximum eigenvalue of $L$. Moreover, the sequence $\{\lambda_{i,l}\}$ converges to a vector which is uniquely determined by the initial conditions of $\lambda_{i}$’s. $\diamond$

Proof:

Substituting the dual feasibility equation to the primal feasibility equation of (7) yields

[TABLE]

Now let $e^{\lambda}_{l}=\lambda_{l}-\lambda^{*}$. Then, one obtains

[TABLE]

From the identity (15), we have

[TABLE]

Here, $\tilde{A}_{\lambda}$ is a symmetric positive semi-definite matrix which has $n$ simple zero eigenvalues, and it holds that $I-\alpha_{\lambda}\sigma_{\max}(\tilde{A}_{\lambda})I\leq I-\alpha_{\lambda}\tilde{A}_{\lambda}\leq I-\alpha_{\lambda}\sigma_{\min}(\tilde{A}_{\lambda})I$. Since $\sigma_{\min}(\tilde{A}_{\lambda})$ is zero, it follows that if $\alpha_{\lambda}>0$ is chosen such that $\alpha_{\lambda}\sigma_{\max}(\tilde{A}_{\lambda})<2$, all eigenvalues of $I-\alpha_{\lambda}\tilde{A}_{\lambda}$, except $1$, are located inside the unit circle. The bound (14) ensures this.

Regarding the convergence of $\lambda_{l}$, we proceed as follows. With the orthonormal matrix $U$ used in Lemma 1, $\tilde{A}_{\lambda}$ can be written as

[TABLE]

where ${M}_{\text{sub}}\in\mathbb{R}^{(N-1)n\times(N-1)n}$ is a submatrix with the first $n$ rows and first $n$ columns removed. In the new coordinates $\bar{e}^{\lambda}_{l}$, defined by $\bar{e}^{\lambda}_{l}=(U^{\top}\otimes I_{n})e^{\lambda}_{l}$, the error dynamics of the dual variable can be expressed as

[TABLE]

From this equation, we know that the first $n$ components of $\bar{e}^{\lambda}_{l}$, denoted by $\tilde{e}^{\lambda}_{l}$, remains the same for any $l$, i.e., $\tilde{e}_{l}^{\lambda}=\tilde{e}_{0}^{\lambda}$, $\forall l\geq 0$, meaning that $(U_{1}^{\top}\otimes I_{n})e^{\lambda}_{l}=\tilde{e}_{0}^{\lambda},\ \forall l\geq 0$, which means that $\tilde{e}_{0}^{\lambda}=(U_{1}^{\top}\otimes I_{n})e^{\lambda}_{0}$. Moreover, with $\alpha_{\lambda}$ chosen as (14), which guarantees that the matrix $I-\alpha_{\lambda}{M}_{\text{sub}}$ has all its eigenvalues except 1 inside the unit circle, we have $\lim_{l\rightarrow\infty}\bar{e}^{\lambda}_{l}=\begin{bmatrix}\tilde{e}^{\lambda}_{0};0\end{bmatrix}$, from which it follows that

[TABLE]

Recalling that $e^{\lambda}_{l}:=\lambda_{l}-\lambda^{*}$, we have from (17)

[TABLE]

Applying $\lambda^{*}=(U_{1}\otimes I_{n})\tilde{\lambda}^{*}+(\bar{U}\bar{\Lambda}^{-1}\bar{U}^{\top}\otimes I_{n})b$ (for $\tilde{\lambda}^{*}$ and $b$, see the proof of Lemma 1), we have

[TABLE]

where $\text{avg}(\lambda_{i,0})=\frac{1}{N}\sum_{i=1}^{N}\lambda_{i,0}$, and this completes the proof. ∎

Now, we derive an update rule of the error covariance matrix. With the information matrix $\Omega_{k}:=P_{k}^{-1}$, the error covariance update rule (3) can be written as

[TABLE]

Define $\Omega_{i,k}:=H^{\top}_{i}R^{-1}_{i}H_{i}+\frac{1}{N}\Omega_{i,k|k-1}$. Then, the updated information matrix of CKF can be obtained by solving the following distributed optimization problem

[TABLE]

where $\zeta_{i}\in\mathbb{R}^{n(n+1)/2}$ is the decision variable. Note that the minimizer $\zeta^{*}:=[\zeta_{i}^{*};\dots;\zeta_{N}^{*}]\in\mathbb{R}^{Nn(n+1)/2}$ of the above optimization problem is nothing but the average of all $\text{vec}(N\Omega_{i,k})$, which corresponds to $\Omega_{k}$.

Define the Lagrangian for the problem (18) as

[TABLE]

where $\mu\in\mathbb{R}^{Nn(n+1)/2}$ is the dual variable. The saddle point equation for (19) is given by

[TABLE]

where $\bar{z}_{\Omega,k}:=[\text{vec}_{h}({N\Omega_{1,k}});\dots;\text{vec}_{h}({N\Omega_{N,k}})]$, and $\mu^{*}$ is the dual variable of the optimal point. From the similar arguments in the proof of Lemma 1, we have

[TABLE]

This implies that the optimal solution $\zeta^{*}_{i}$ is the half vectorization of the average of the locally predicted information matrix corrected by the global information $H^{\top}R^{-1}H$.

Based on the above arguments, we propose a dual ascent type update rule for the error covariance matrix as

[TABLE]

where $\alpha_{\mu}$ is a step size such that $0<\alpha_{\mu}<2/\sigma_{N}^{2}$, which obtained by the similar arguments in the proof of Lemma 2. Putting all pieces together, we propose a DKF algorithm described in Algorithm 1.

In the structural point of view, the algorithm consists of local prediction step and distributed correction step as in CKF. In the local prediction step, each estimator locally predicts the estimate and the corresponding covariance matrix. In the distributed correction step, each estimator finds the optimal points for the state estimate and its error covariance matrix, iteratively, by using the local measurement information and exchanging information with its neighbors. With sufficiently large $l^{*}$, locally updated $\xi_{i,l^{*}}$ and $P_{i,k}$ converge to those of CKF with tunable size of errors.

IV Numerical Experiments

We have two examples for the developed theory. The first one is a simple academic example, while the second one is more practical one.

Example 1

Consider a system given by

[TABLE]

and $Q=0.1$, $R=\operatorname{\text{diag}}(0.1,0.2,0.3,0.1)$, and suppose that $4$ estimators are connected through a communication network whose Laplacian matrix is given by

[TABLE]

The step sizes for the algorithm are chosen as $\alpha_{\lambda},\alpha_{\mu}=0.01$.

Figure 1(a) shows that the average error norm defined by $\text{avg}(\|e_{i,k}\|)=\frac{1}{N}\sum_{i=1}^{N}\|\hat{x}_{i,k}-x_{k}\|$ decreases more rapidly as $l^{*}$ increases. Figure 1(b) also shows that as $l^{*}$ increases, the average error covariance norm defined by $\text{avg}(\|P_{i,k}\|):=\frac{1}{N}\sum_{i=1}^{N}\|P_{i,k}\|$ approaches $\|P_{k}\|$ which is the norm of the error covariance matrix of CKF. It is seen that, when $l^{*}=50$, there is very little difference between $\text{avg}(\|P_{i,k}\|)$ and $\|P_{k}\|$ of CKF. $\square$

Example 2

In this example, we evaluate DA-DKF with a network consisting of 50 estimators to estimate the state of a target system. The dynamics of the target system is described by

[TABLE]

where $w_{k}\sim\mathsf{N}(0,Q)$ and $Q=0.1$. The first and the third components of $x_{k}$ represent the $x$-axis position and $y$-axis position in the plane, respectively .

The estimator $i$ knows $e^{A}$, $Q$, $H_{i}\in\mathbb{R}^{1\times 4}$ and $R_{i}>0$, and each $H_{i}$ and $R_{i}$ is randomly chosen. The connections among estimators are also randomly selected and the weight is $1$ when connected, and the maximum eigenvalue of $L$ is $18.5$. For all $i$, $P_{i,0}=Q$ and each component of the initial estimate $\hat{x}_{i,0}$ is randomly chosen within $(-15,15)$ as shown in Figure 2(a). The parameters for DA-DKF were chosen as $\alpha,\beta=10^{-5}$, $l^{*}=10$. Figure 2 shows four snapshots of the target system’s position (black cross) and the each estimator’s estimate (red circles). The blue line is the trajectory of the target system. As time goes by (as $k$ increases), the estimates of the distributed Kalman-filters converge to the vicinity of the position of the target system. $\square$

V Conclusions and Future Work

This paper dealt with DKF from the optimization perspective. By observing that the correction step of Kalman-filtering is basically an optimization problem, we formulated DKF problem from the centralized one. The formulated problem is a quadratic consensus optimization problem. One of the recent DKF algorithms, Consensus on Information [14] was reinterpreted from the distributed optimization perspective. In addition, various DKF algorithms can be derived, by employing many existing distributed optimization methods to DKF problem. As an instance, DA-DKF has been presented, employing the distributed dual ascent method, and the algorithm has been validated with numerical experiments.

For the future work, we plan to analyze the effect of the residuals of the previous iteration $k$, especially how the residuals affect the convergence. In addition, researches on more practical obstacles, such as considering the time-varying network topology, reducing communication loads, will be conducted.