Simultaneous Perturbation Stochastic Approximation-based Consensus for Tracking under Unknown-but-Bounded Disturbances

Simultaneous Perturbation Stochastic

Approximation-based Consensus for Tracking under Unknown–but–Bounded Disturbances

Oleg Granichin, Senior Member, IEEE, Victoria Erofeeva, Yury Ivanskiy, and Yuming Jiang,Senior Member, IEEE

Abstract—We consider a setup where a distributed set of sensors working cooperatively can estimate an unknown signal of interest, whereas any individual sensor cannot fulfil the task due to lack of necessary information diversity. This paper deals with these kinds of estimation and tracking problems and focuses on a class of simultaneous perturbation stochastic approximation (SPSA)-based consensus algorithms for the cases when the cor- rupted observations of sensors are transmitted between sensors with communication noise and the communication protocol has to satisfy a prespecified cost constraints on the network topology.

Sufficient conditions are introduced to guarantee the stability of estimates obtained in this way, without resorting to commonly used but stringent conventional statistical assumptions about the observation noise such as randomness, independence, and zero mean. We derive an upper bound of the mean–square error of the estimates in the problem of unknown time-varying parameters tracking under unknown–but–bounded observation errors and noisy communication channels. The result is illustrated by a practical application to the multi-sensor multi-target tracking problem.

Index Terms—Distributed tracking, multi-agent networks, con- sensus algorithm, simultaneous perturbation stochastic approxi- mation, SPSA, randomized algorithm, arbitrary noise, unknown–

but–bounded disturbances, stochastic stability, tracking perfor- mance.

I. INTRODUCTION

Distributed cooperative control of networked systems has been investigated and numerous potential applications to com- plex manufacturing, energy and social systems have been developed [1]–[3] over the past few decades. One of the fundamental concepts in multi-agent cooperative control is consensus. This approach aims to find an agreement between all agents in a network regarding a common value using only local information and communicating among neighboring agents.

The goal of distributed optimization is usually to find the minimum of some loss function F¯(x) = Pn

i=1Fⁱ(x) via interaction between agents. Here, x ∈ R^d and Fⁱ(x) :

The theoretical research in Sections I-VI of this work was supported by the Russian Fund for Basic Research (project no. 20-01-00619). The obtaining of experimental results in Section VII was supported by Russian Science Foundation (project no. 19-71-10012).

O. Granichin, V. Erofeeva and Y. Ivanskiy are with Saint Petersburg State University, 198504, Universitetskii pr. 28, St. Petersburg, Russia. e-mail:

[email protected], [email protected], [email protected].

Y. Jiang is with Department of Information Security and Communication Technology, Norwegian University of Science and Technology (NTNU), NO- 7491, Trondheim, Norway, e-mail: [email protected].

R^d → R is the loss function of agent i, which is typically known only to the agent itself. Studies of consensus and distributed optimization algorithms began from the 1970-80s [4], [5]. Distributed asynchronous stochastic approximation algorithms were studied in [6]. To date, there exist a number of approaches for the case when functionsFⁱ(x)are convex. In particular, the Alternating Direction Method of Multipliers [7], [8], as well as the subgradient method [9], [10] were proposed.

For non-convex tasks, the works [11], [12] develop a large class of distributed algorithms based on various “functional- surrogate units”. The distributed tracking problem is consid- ered when the estimated parameters vary over time.

Recently, for large-scale systems consisting of many in- dividuals (components, targets), a distributed optimization is often used to estimate the unknown parameters which mini- mize a loss function, based on the information obtained by distributed sensors. So-calledmultitarget-multisensor tracking problems have been widely studied in many practical ap- plications such as adaptive mobile networks, cognitive ra- dio systems, target localization in biological networks, fish schooling, bee swarming, and bird flight (see, e.g., [13], [14]). It is well known that distributed tracking algorithms have some significant advantages over the centralized ones or the fusion methods. Centralized algorithms usually require the distributed sensor network to transmit the whole system information into a fusion center to estimate the unknown signals. This leads to the necessity of strong communication capabilities over sensor networks which is hard to provide in many practical situations when sensors may only have the capability to exchange information locally between their neighbors. An alternative approach for multitarget-multisensor tracking problems assumes only local interaction between sensors without the governing data fusion center. A detailed lit- erature overview of the recent advances in the stability analysis of a consensus-based least squares algorithms is performed in [15] for distributed estimation and tracking problems.

The maximum likelihood estimator and stochastic approxi- mation (SA) algorithms with decreasing to zero step-size are actively used to optimize mean-risk functionals. In gradient- free conventional stochastic approximation algorithms, two measurements are used to approximate each component of the vector-gradient of the cost function (implying 2d mea- surements for the d-dimension state space). Simultaneous perturbation stochastic approximation (SPSA) was proposed by Spall [16]. It can be used to solve optimization problems

in the case when it is difficult or impossible to obtain a gradient of the loss function with respect to the parameters being optimized. In any multidimensional case (d1), SPSA requires only two measurements of a loss function at each iteration. The current estimate is updated along a randomly chosen direction∆with±1Bernoulli distributed components.

Traditionally, a stochastic optimization problem under un- certainties focused on finding a set of system parameters that deliver a minimum (or maximum) value to a certain mean-risk functional. In practical applications, these parameters may also vary over time. The problem of tracking changes in system parameters is considered in [17]–[19]. In this paper, such a problem is called the minimum-point tracking of a nonsta- tionary mean-risk functional. In centralized (non-distributed) cases, SPSA-like algorithms for parameter tracking problems were considered in [20]–[22]. The stochastic approximation method with a constant step-size has also been used in [23]

to achieve the approximate mean-squared consensus in multi- agent systems operating under noisy measurement conditions.

Contributions. In the case of differentiable time-varying loss functions and almost arbitrary external bounded noise, an upper bound of the mean square estimation error was derived in [20] for estimates of the SPSA type algorithms with constant step-size. This upper bound may be sufficiently small compared to the significant level of noise when the rate of change of parameters is low enough. One of the main conditions is a strong convexity property of the minimized mean-risk functional. In this paper we weaken this assumption by combining SPSA with the consensus algorithm from [23].

We propose a new SPSA-based consensus algorithm for dis- tributed tracking under unknown–but–bounded disturbances.

The preliminary concept of this paper is presented in [24].

In many practical applications, the network processes the data under certain constraints, and the data transmission is accompanied by noise. In this paper, compared with [24], we consider such noisy data transmission and a communication protocol with prespecified cost constraints on the network topology. Also, we study a more general type of simultaneous perturbation and we choose the current points of observations in a more general manner. We obtain an upper bound of the mean square error of estimates of unknown time-varing param- eters tracking. Communication cost constraints are satisfied by exploiting a specific intentionally randomized topology of the network communication graph.

The paper is organized as follows. The preliminary infor- mation regarding concepts of the graph theory and network topology is given in Section II. A formal problem setting of a distributed non-constrained time-varying mean-risk optimiza- tion with noisy local communications is given in Section III.

The main result including assumptions and the SPSA-based consensus algorithm for tracking is presented in Section IV.

In Section V, the efficiency of the proposed algorithm is illustrated through the numerical simulation.

II. PRELIMINARIES

Let (Ω,F, P) be the underlying probability space cor- responding to sample space Ω, set of all events F, and probability measureP.Edenotes mathematical expectation.

A. Concepts of Graph Theory

Given a network consisting of n sensors. Let the inter- action between sensors be described by the directed graph G = (N,E), where N ={1, . . . , n} is a set of vertices and E ⊆ N × N is a set of edges. A subgraph of G is a graph G¯= (NG¯,EG¯), whereNG¯⊆ N andEG¯⊆ E. Denote byi∈ N an identifier ofi-th sensor and(j, i)∈ E if there is a directed edge from sensorj to sensori. The latter means that sensorj is able to transmit data to sensor i. For a sensor i∈ N, the set of neighbors is defined as Nⁱ ={j ∈ N : (j, i) ∈ E}.

Thein-degreeofi∈ N equals|Nⁱ|. Here and after,| · |is the cardinality of a set, and the identifier ofi-th sensor is used as a superscript and not as an exponent.

Let c^i,j > 0 be the weight associated with the edge (j, i)∈ E, and c^i,j = 0 whenever (j, i)∈ E/ . Let C = [c^i,j], be the weighted adjacency matrix, or simply connectivity matrix. Denote by G_C = (N_C,E_C) the weighted directed graph, where NC ≡ N and EC ≡ E. We assume that weightc^i,j is the cost of data transmission through the edge (j, i) ∈ EC. The weighted in-degree of i ∈ NC is defined as deg⁺_i(C) =Pn

j=1c^i,j, the maximum in-degree among all nodes contained in the graphGC asdeg⁺_max(C), and the diag- onal matrix as D(C) = diag_n(col(deg⁺₁(C), . . . ,deg⁺_n(C))).

Then, L(C) =D(C)−C is theLaplacianof graphGC. Definition 1. A directed graph G is said to be strongly connected if for every pair of nodes j, i ∈ N, there exists a path of directed edges that goes from j toi.

Denote the eigenvalues of Laplacian L(C) by λ₁, . . . , λ_n and arrange them in ascending order of real parts: 0 ≤ Re(λ₁) ≤ Re(λ₂) ≤ . . . ≤ Re(λ_n). It is known, that if the graph is strongly connected then λ1 = 0 and all other eigenvalues ofLare in the open right half of the complex plane (see, e.g., [3]). The eigenvalue of matrix C with maximum absolute magnitude is defined asλmax(C).

B. Network Topology Constraints

In practice, we have constraints on the bandwidth of com- munication channels, network response time, hardware and financial requirements, etc. In this paper, we associate these constraints with matrixC, which characterizes the cost of data transmission in the network. In many practical applications, we may represent cost constraints of sensor i ∈ N as a predefined upper boundQ:deg⁺_i(C)≤Q. This bound may be thought of as the total cost of communication with neighbors of sensor i. To satisfy this constraint, we may generate at each time instant t subgraph GB_t ⊂ GC associated with the weighted connectivity matrix Bt such that deg⁺_i (Bt) ≤ Q.

Obviously, the cost constraint deg⁺_i (Bt) ≤ Q may not be satisfied for given Bt = C andQ, e.g. when n= 6,GC is the complete graph withc^i,j = 1,i6=j,c^i,i= 0, and Q <5.

One possible solution is to use a randomized topology, when we drop5−Qedges randomly. Such randomized strategy for Q= 1is similar to the scheme used ingossipalgorithms [25].

Moreover, random subgraphs naturally arise in many practical applications.

Next, we consider a communication protocol needed to satisfy a predefined averaged cost constraint.

Definition 2. Random subgraphG_Bt satisfies the averaged cost constraints with levelQif

E deg⁺_max(B_t)≤Q. (1) In the example considered above we are able to satisfy averaged cost constraints if each sensorirandomly selects its neighbors out of allj∈ Nⁱwith probability ^Q

deg⁺_i(Bt)= 0.2Q at each time instant t.

III. DISTRIBUTEDTRACKING

A. Non-stationary Mean-risk Functional

Let Ξ be a set, {f_ξⁱ(θ)}_ξ∈Ξ be a family of differentiable functions: ∀i ∈ N f_ξⁱ(θ) : R^d → R. We assume that parameterθcannot be directly measured. Hence, we introduce a sequence of measurement points xⁱ₁,xⁱ₂, . . . , i∈ N chosen according to an observation plan. The values yⁱ₁, y₂ⁱ, . . . of the functions f_ξⁱ

t(·)are observable at every time instant t= 1,2, . . .with additive externalunknown-but-bounded noisevⁱ_t yⁱ_t=f_ξⁱ_t(xⁱ_t) +v_tⁱ, (2) where {ξt}, ξt ∈Ξ, is a non-controllable deterministic (e.g., Ξ =Nandξt=t) or random sequence. In the latter case we assume that a probability distribution ofξtexists and may be known or unknown.

LetFt−1be theσ-algebra of all probabilistic events which happened up to time instant t,EFt−1 denotes the conditional mathematical expectation with respect to the σ-algebraFt−1. We consider an optimization problem in which the cost function F¯_t(θ)is expressed as the sum of local contributions F_tⁱ(θ) = EFt−1f_ξⁱ

t(θ) and all of them depend on a common optimization variable θ. Moreover, minimizer θ of F¯t(θ) may vary over time. Formally, the non-stationary mean-risk optimization problem is as follows: estimate the time-varying minimum point θtof the distributed function

F¯t(θ) =X

i∈N

F_tⁱ(θ) =EFt−1

X

i∈N

f_ξⁱ_t(θ)→min

θ . (3)

More precisely, the problem is to obtain an estimate θbt of an unknown vectorθtbased on the observationsyⁱ₁, y₂ⁱ, . . . , yⁱ_t and measurement points xⁱ₁,xⁱ₂, . . . ,xⁱ_t, i ∈ N, through the minimization of time-varying mean-risk functional(3).

Remark. There exist two special cases of measurement model (2) related to the different types of noise vⁱ_tand distur- bance ξt: (i) If drift θt is deterministic thenF_tⁱ(θ) = f_ξⁱ_t(θ), ξ_t=t, and the measurement model may be defined in more conventional way as y_tⁱ = F_tⁱ(xⁱ_t) +vⁱ_t; (ii) If noise v_tⁱ has a probability distribution then we may consider it as a part of disturbance ξ_t. The measurement model for this case is yⁱ_t=f_ξⁱ

t(xⁱ_t).

B. Communication Network

In centralized networks, it is required to transmit all needed information such as y₁ⁱ, yⁱ₂, . . . , y_tⁱ,xⁱ₁,xⁱ₂, . . . ,xⁱ_t, i ∈ N, to a fusion center in order to estimate the unknown vectorθt. In such networks, robustness of the fusion center and quality of the communication channels become a critical factor. In many

situations, sensors may only have the capability to exchange information with their neighbors. The communication with neighbors may be much cheaper or much faster then trans- mission to fusion center as well. Moreover, the information may be transmitted over the noisy communication channels and with delays, and the network topology may vary over time. Also, in practice, the network cost constraints naturally arise. For example, we don’t have communication channels with infinite bandwidth and the response time of the network should be practically reasonable. All these factors motivate the development of distributed decentralized algorithms.

To formalize the distributed setting, we assume that at time instanttsensors are able to communicate with their neighbors through the network defined by graphGBt= (NBt,EBt). The corresponding connectivity matrixB_tsatisfies some averaged cost constraints (1) with levelQ.

We also assume that sensoriobtains its current estimateθb_tⁱ based on its noisy observation (2) and, if the set N_tⁱ ={j ∈ N_Bt : (j, i)∈ E_Bt}is not empty, also on the current estimates transmitted by its neighbors through the noisy channels

θ˜_t^i,j=θb^j_t+w^i,j_t , j∈ N_tⁱ, (4) wherew^i,j_t is communication noise. Ifj6∈ N_tⁱwe setθ˜^i,j_t = 0.

C. Example

In this subsection, we present an example illustrating the considered problem statement. Given a distributed network consisting of n = 6 planar sensors identified by i ∈ N = {1,2, . . . ,6}. The state of sensor i is sⁱ ∈ R². We assume that the states are known and doesn’t depend on time, i.e.

the sensors are stationary. In the sensing range of the sen- sors, there are m = 10 moving planar targets identified by l ∈ M = {1,2, . . . ,10}. The goal of each sensor i is to estimate the states of all targets r^l_t∈R² at time instantt.

Letθt= col(r¹_t, . . . ,r¹⁰_t )∈R²⁰be the common state vector of all targets,θˆ_tⁱ= col(ˆr^i,1_t , . . . ,ˆr^i,10_t )be a common vector of i-th sensor current estimates. Each target l∈ M changes the position according to the following dynamics:

r^l_t=r^l_t−1+ζ_t−1^l , l∈ M, (5) whereζ_t−1^l are random vectors uniformly distributed in a ball.

We assume that at time instantt sensor i is able to measure the squared distance ρ^i,l_t = ρ(sⁱ,r^l_t) = kr^l_t−sⁱk² to some moving target r^l_t.

The network is modelled by complete graphGC, for which we have the followingtopology constraints: each sensori∈ N at each time instant t is able to measure the noisy squared distance to only one target l ∈ M and to receive estimates θb^j_t and measurements ρ^j,l_t only from one randomly chosen neighbor j ∈ N_tⁱ. This leads to the communication protocol satisfying averaged cost constraints with level Q= 1 consid- ered as example at the end of Section II.

Let sensori receive the current estimate and measurement from some neighbor with identifier j ∈ N. Denote by u = col(i, j, l) the vector, where the first element is the identifier of a sensor, the intermediate element is the identifier of a neighbor, which shares its information with sensor i,

and the last element is the identifier of a target, which this sensor observe at time instant t. Note that in general there may be several intermediate elements. Also, denote by ρ¯t(u) = ρ(sⁱ,r^l(u)_t )−ρ(s^j,r^l(u)_t ) the residual between measurements of sensor i and its neighborj. Here and after, l(u) : R^|u| → R gives the last element of u. In this case, using the square difference formula we derive

C^ur^l_t=D^u_t ⇒C^uTC^ur^l_t=C^uTD^u_t ⇒I^ur^l_t=H_t^u, (6) where I^u = [C^uTC^u]⁰C^uTC^u, H_t^u = [C^uTC^u]⁰C^uTD^u_t, C^u = 2(s^j −sⁱ)^T, D_t^u = ¯ρt(u) +ks^jk²− ksⁱk², and [·]⁰ denotes a vector or matrix Moore–Penrose inverse.

Denote byUⁱthe set of all vectorsuwith the first elementi.

Letuⁱ_t∈Uⁱ be a random variable and inputxⁱ= ˆθⁱ_tbe fixed.

We consider observation model (2) as follows y_tⁱ=f_ξⁱ

t(xⁱ) =kI^uitˆr^i,l(uit)−H_t^uitk², (7) where ξt consists of all uⁱ_t generated at time instant t, i.e.

ξt= col(θt,u¹_t,u²_t, . . . ,u⁶_t).

This leads us to following individual mean-risk sub- functionals F_tⁱ(xⁱ) = EFt−1f_ξⁱ

t(xⁱ), which are equal to

1

|Uⁱ|

P

uⁱ∈UⁱkI^uiˆr^i,l(ui)−H_t^uik²when positions of all targets do not evolve over time.

IV. MAINRESULT

In this section, we present the main result of this paper. All proofs are relegated to Appendix.

A. SPSA-based Consensus Algorithm

Let∆ⁱ_k, k= 1,2, . . . , i∈ N,be an observed sequence of independent random vectors inR^d, called thesimultaneous test perturbation, with symmetrical distribution functions Pⁱ_k(·), and let Kⁱ_k(·) :R^d → R^d, k = 1,2, . . . , be a set of vector functions (kernels).

Let us take fixed nonrandom initial vectorsθb₀ⁱ ∈R^d, positive step-size α, gain coefficientγ, and choose sequences of such nonnegative numbers{β_k⁺}and{β⁻_k}thatβk =β_k⁺+β_k⁻>0.

We consider the algorithm with two observations of distributed sub-functions f_ξⁱ

t(θ) for each agent i ∈ N for constructing sequences of measurement points {xⁱ_t} and estimates {bθ_tⁱ}:















xⁱ_2k=θbⁱ_2k−2+β_k⁺∆ⁱ_k,xⁱ_2k−1=θb_2k−2ⁱ −β_k⁻∆ⁱ_k, bθⁱ_2k−1=θb_2k−2ⁱ ,

bθⁱ_2k=θbⁱ_2k−1−α_yⁱ

2k−yⁱ_2k−1

βk Kⁱ_k(∆ⁱ_k)+

γP

j∈N_2k−1ⁱ b^i,j_2k−1(˜θ_2k−1^i,j −θbⁱ_2k−1) .

(8)

Algorithm (8) consists of two parts: (i) The first one is similar to SPSA-like algorithm from [20]. This part represents a pseudo-gradient of sub-functions f_ξⁱ

t(θ); (ii) The second one coincides with Local Voting Protocol (LVP) from [23], where it is used to solve load balancing problem in stochastic networks. This part is determined for each sensor i by the weighted sum of differences between the information about the current estimate bθⁱ_2k−1 of sensor iand noisy information about the estimates of its neighbors.

Further, we denote by θ¯_t = col(bθ¹_t, . . . ,θb_tⁿ) the com- mon vector of estimates of all sensors at time instant t and by θ˜¯_t = col(˜θ^1,1_t ,θ˜^2,1_t , . . . ,θ˜^n,1_t ,θ˜^1,2_t , . . . ,θ˜_t^n,n) the cor- responding vector of data transmitted over the noisy chan- nels. Also we introduce the following notations: y¯t = diag_n(col(y¹_t, . . . , yⁿ_t)),∆¯k = col(K¹_k(∆¹_k), . . . ,Kⁿ_k(∆ⁿ_k)).

Using the notations introduced above, the algorithm (8) can be rewritten in the following form

θ¯2k= ¯θ_2k−1−α

y¯2k−y¯_2k−1 β_k ⊗Id

∆¯k+ γ

L¯2k−1⊗I_d θ˜¯_2k−1

(9) where(n×n²) matrix L¯_2k−1 is defined in such a way that its i-th row consists of zeros except the elements from the position (j−1)n+ 1tojn which coincide withi-th row of L(B_2k−1).

B. Main Assumptions

For anyi∈ N let us formulate assumptions about functions F_tⁱ(x), f_ξⁱ(x), disturbances, network topology, randomized perturbations∆ⁱ_k, and noises.

1: The functions F_tⁱ(·) are convex, they have a common minimum pointθt and

∀x∈R^d hx−θt,EFt−1∇f_ξⁱ_t(x)i ≥0.

2:∀ξ∈Ξ,∀i∈ N the gradient∇f_ξⁱ(x)satisfies the Lipschitz condition:∀x₁,x₂∈R^d

k∇f_ξⁱ(x₁)− ∇f_ξⁱ(x₂)k ≤Mkx1−x₂k with the same constantM >0.

3: The gradient ∇f_ξⁱ_t is uniformly bounded in the mean- square sense at the minimum pointθt:∀tEk∇f_ξⁱ

t(θt)k²≤g²₂, Eh∇f_ξⁱ

t(θt),∇f_ξⁱ_t−1(θ_t−1)i ≤g₂²(g2= 0ifξtis not a random parameter, i.e. f_ξⁱ

t(x) =F_tⁱ(x)).

4: The drift is bounded: a) kθt −θ_t−1k ≤ δθ < ∞, or Ekθt−θ_t−1k² ≤ δ_θ² and Ekθt−θ_t−1kkθ_t−1−θ_t−2k ≤ δ_θ² if a sequence{ξt} is random;

b) EF2k−2|f_ξⁱ

2k(x)−f_ξⁱ_2k−1(x)|^q ≤δ_θ^q(g^q₀+g₁^qkx−θ_2k−2k^q) for q= 1,2 and for anyi∈ N.

5: a) Graphs GB_t, t= 0, . . .are i.i.d., i.e. the random events of appearance of “time-varying” edge (j, i) in graphGB_t are independent and identically distributed for the fixed pair(j, i), i∈ N, j∈ N_maxⁱ =∪tN_tⁱ.

b) For all i∈ N, j ∈ N_tⁱ weights b^i,j_t are independent ran- dom variables with mean values (mathematical expectations):

Eb^i,j_t = b^i,j_av, and bounded variances: EkB_t−B_avk² ≤ σ_B² whereB_av = [b^i,j_av].

c)EP

j∈N_tⁱ(b^i,j_t )²≤ _n−1^Q2 .

d) GraphGB_av is strongly connected.

6: For k = 1,2, . . . , the successive differences ˜v_kⁱ =vⁱ_2k − v_2k−1ⁱ of observation noise are bounded:|˜vⁱ_k| ≤ cv <∞, or E(˜v_kⁱ)²≤c²_v if a sequence{˜vⁱ_k} is random.

7:Fort= 1,2, . . . ,∀i∈ N,∀j∈ N the communication noise w_t^i,jis random i.i.d. (independent identically distributed) with zero-meanEw^i,j_t = 0and bounded disturbances:Ekw_t^i,jk²≤

σ_w². All random vectors and valuesw_t^i,j,b^i,j_t ,ξ_t, andξ_t+1 are mutually independent (if they are random).

8: For anyi, j∈ N, k= 1,2, . . . , a) Vectors∆ⁱ_k are mutually independent.

b) ∆ⁱ_k andξ_2k−1, ξ2k (if they are random) do not depend on the σ-algebraF_2k−2.

c) If ξ_2k−1, ξ2k,¯v_kⁱ,w^i,j_2k−1, b^i,j_2k−1 are random, then random vectors ∆ⁱ_k and elements ξ_2k−1, ξ2k,˜v_kⁱ, w_2k−1^i,j , b^i,j_2k−1 are independent.

d) For k = 1,2, . . . , vectors ∆ⁱ_k and vector functions Kⁱ_k(·) along with simultaneous perturbation symmetrical distribution functionsP_k(·)satisfy the conditions

Z

xP_k(dx) = Z

xkKⁱ_k(x)k²P_k(dx) = Z

Kⁱ_k(x)P_k(dx) = 0, Z

he,xiKⁱ_k(x)Pk(dx) =he,1di1d, Z

kxk²Pk(dx)≤c²_∆,(10) Z

kKⁱ_k(x)k²Pk(dx)≤c²_∆, Z

kKⁱ_k(x)k²kxk²Pk(dx)≤c⁴_∆. Note that all Assumptions 1–8 are standard for the consid- ered problem.

Remark.Usually, it is practically reasonable to define{∆ⁱ_k} as a sequence of independent Bernoulli random vectors from R^d with each component independently taking values ±^√1₂ with probabilities ¹₂ andKⁱ_k(x)≡xas kernel functions. For this case, we have c_∆ = 1. The case, when β_k⁺ = β_k⁻ and decreasing to zero sequence α_k is used instead of constant step-size α, corresponds to the SPSA algorithm in [16].

The similar algorithm with randomly varying truncations and randomized difference was studied in [26] where the case β_k⁻= 0was additionally considered.

Example. Return back to the example from Section III-C and check Assumptions 1–5.

1. Using (6) and (7), we obtain for gradient hx−θ_t,EFt−1∇f_ξⁱ

t(x)i=EFt−1(x^i,l(uit)−r^l(uit))^T[I^ui]^T I^ui(x^i,l(uit)−r^l(uit))≥0.

2. Using (7), we obtain k∇f_ξⁱ(x₁) − ∇f_ξⁱ(x₂)k = k2[I^ui]^TI^ui(x^l(u₁ ⁱ⁾−x^l(u₂ ⁱ⁾)k ≤ Mkx2−x1k, where M = maxik2[I^ui]^TI^uik.

3. ∇f_ξⁱ_t(θt) = 0. Hence,g2= 0.

4. Assumption about the drift holds for δ_θ = nδ_ζ and by virtue of drift model (5) when ζ_tⁱ are random i.i.d. vectors with Eζ_tⁱ = 0, and Ekζtk² ≤δ_ζ²,g₀ = 4√

2¯s²,g₁= 8√ 2¯s², wheres¯= maxi,jksⁱ−s^jk.

5. a), c), d) hold by the construction; in b) b^i,j_av = 0.2, i6=j, b^i,i_av = 0,σ_B² = 4.8.

C. Analysis of the Estimation Accuracy

To analyze the quality of estimates we apply the following definition for the problem of minimum tracking for mean-risk functional (3):

Definition 2.A sequence of estimates{θ¯_2kⁱ } hasan asymp- totically efficient upper boundL >¯ 0of residuals of estimation if∀ε >0∃k¯ such that∀k >¯k

q

Ekθ¯ⁱ_2k−1_n⊗θ_2kk²≤L¯+ε.

Denote λ¯2 = Re(λ2(L(Bav))), λ¯m = λ

1

max2 (L(Bav)^TL(Bav)), c+ = maxk β_k⁺

β_k, β˜ = maxk 1 β_k,

¯

c = maxk

_β+ k

β_k

² + _β−

k

β_k

²

, β¯ = maxk (β_k⁺)²

β_k + ^(β

− k)² β_k , cm = λ¯²_m + σ²_B, c1 = c∆λ¯mM(δθg1β˜ + c∆), c2 = 2c²_∆(δ_θ²g₁²β˜² + c²_∆M²), cµ = (¯λ2 − αc1)/cm, cd =p

1−α²c2cm/(¯λ2−αc1)².

The following theorem shows the asymptotically efficient upper bound of estimation residuals provided by algorithm (8).

Theorem 1: If Assumptions 1–8 hold, β¯ = minkβk > 0, positive constantαis sufficiently small:

α <

λ¯₂ c1+√

c2cm

(11) and

c_µ(1−c_d)< αγ < c_µ(1 +c_d) (12) thenthe averaged cost constraint (1) holds and the sequence of estimates provided by algorithm (8) has an asymptotically efficient upper bound which equals to

L¯= 1 µ

h+p

h²+lµ

, (13)

whereµ= 2γ¯λ₂−α(c_mγ²+α(2γc₁+c₂)), h=γc₃+c₄, l=αγ²Q²σ²_w+c₅,

c₃= 2√

nλ¯_mδ_θ+αλ¯_mc_∆M(δ_θg₀β˜+ ¯βc²_∆), c₄=M c₊+c⁴_∆g₁¯c+ 2c²_∆(1 +c₊)

δ_θ+c⁵_∆M²β,¯ c5= ^4nδ_α^2θ + 2αc²_∆

β˜²n(c²_v+δ²_θg₀²) +c²_∆cnM¯ (cv+δθg0)+

c³_∆nMβ(M δ¯ θ+δθc++g2) +

2M n(δ_θ²c++c³_∆β) +¯ c²_∆n δ_θ²(1 +c+)²+g₂²+M²β¯²c⁴_∆ . See the proof of Theorem 1 in Appendix.

Remarks. 1. The boundLin the Theorem 1 is tight, so there exists no L⁰ < L such that the statement of the Theorem 1 still holds if all inequalities from the Assumptions 1–8 are replaced by equation.

2. The observation noisev_tⁱin Theorem 1 can be said to be almost arbitrary since it may either be nonrandom but bounded or it may also be a realization of some stochastic process with arbitrary internal dependencies. In particular, to prove the results of Theorem 1, there is no need to assume that v_tⁱ andF_t−1 are not dependent.

3. The proof of Theorem 1 allows for consideration of random sequences{β_k⁺}and{β⁻_k}whose values at iterationk are measurable under the correspondingσ-algebraF2k−2. This fact is sometimes useful from a practical point of view.

4. The result of the Theorem 1 shows that for the case without drift (δθ = 0) and g2 = 0 under any noise level cv

the asymptotic upper bound can be made infinitely small by choosing sufficiently smallαandβ_k^±. At the same time, in the case of drift, the bigger drift normδθ can be compensated by choosing a bigger step-sizeαandβ_n^±. This leads to a tradeoff between makingαsmaller because of noisy observations and makingαbigger due to the drift of optimal points.

V. SIMULATION

In this section, we present the numerical experiments, which illustrate the performance of the suggested algorithm. We apply the algorithm to the problem described in Section III-C.

The starting positions of the targets are chosen randomly from the interval[0; 100]. The states of the targets evolve over time as follows: r^l_t = r^l_t−1 +χ^l_t−1. Let χ^l_t−1 be a random vector uniformly distributed on the ball of radius equal to 0.2. The sensors don’t move and their coordinates are random values uniformly distributed in interval[100; 120]. We consider observation model (2) defined asy_tⁱ=kI^uitˆr^i,l(uit)−H_t^uitk²+ vⁱ_t,wherev_tⁱis modelled as unknown-but-bounded disturbance falling withing interval [0.6; -0.6].

Algorithm (8) working on each node has the following parameters: α= 0.03,β = 1.5,γ = 1.5. The initial estimate on each sensor for each target coordinate was chosen randomly from the interval [50; 100]. Fig. 1 shows how the residuals evolve over time. Figures show that there exists time instant t starting with which the estimations converge to the actual value and move next to it.

Fig. 1. ResidualskI^uitˆr^i,l(uit)−H^u

i t

t k²obtained by nodes.

VI. CONCLUSION

In this paper, we proposed a new SPSA-based consen- sus algorithm for distributed tracking under unknown–but–

bounded disturbances. Compared to the SPSA algorithm, this algorithm is suitable for distributed problems due to the relaxed assumption regarding the strong convexity of the min- imized mean-risk functional. In many practical applications, the network processes the data under certain constraints, and the data transmission is accompanied by noise. In this paper, we consider such noisy data transmission and the case where a communication protocol has to satisfy prespecified cost constraints. Communication cost constrains are satisfied by exploiting a specific intentionally randomized topology of the network communication graph. We obtain an upper bound on the mean square error of estimates of tracking unknown time- varing parameters under unknown–but–bounded observation errors and noisy communication channels.

-20 0 20 40 60 80 100 120

x -40

-20 0 20 40 60 80 100 120

y

Fig. 2. The estimatesˆr^i,l_t obtained by nodes and actual targets positions r^i,l_t . (Empty circles denote sensor positions, targets movement is depicted as a series of shaded circles and plus signs show the estimated target positions.) The figure shows sparse data for clarity: each 50th position of targets and the estimates.

REFERENCES

[1] R. Olfati-Saber and R. M. Murray, “Consensus problems in networks of agents with switching topology and time-delays,”IEEE Transactions on Automatic Control, vol. 49, no. 9, pp. 1520–1533, 2004.

[2] W. Ren and R. W. Beard,Distributed consensus in multi-vehicle coop- erative control. Springer, 2008.

[3] F. L. Lewis, H. Zhang, K. Hengster-Movric, and A. Das,Cooperative control of multi-agent systems: optimal and adaptive design approaches.

Springer Science & Business Media, 2013.

[4] M. DeGroot, “Reaching a consensus,”J. Am. Stat. Assoc., vol. 69, pp.

118–121, 1974.

[5] V. Borkar and P. Varaiya, “Asymptotic agreement in distributed estima- tion,”IEEE Trans. Autom. Control, vol. 27, pp. 650–655, 1982.

[6] J. N. Tsitsiklis, D. P. Bertsekas, and M. Athans, “Distributed asyn- chronous deterministic and stochastic gradient optimization algorithms,”

in1984 American Control Conference, 1984, pp. 484–489.

[7] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Ecksteinet al., “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends in Machine Learning, vol. 3, no. 1, pp. 1–122, 2011.

[8] A. Falsone, I. Notarnicola, G. Notarstefano, and M. Prandini, “Tracking- admm for distributed constraint-coupled optimization,”Automatica, vol.

117, p. 108962, 2020.

[9] M. Rabbat and R. Nowak, “Distributed optimization in sensor networks,”

in Proceedings of the 3rd international symposium on Information processing in sensor networks. ACM, 2004, pp. 20–27.

[10] A. Nedic and A. Ozdaglar, “Distributed subgradient methods for multi- agent optimization,”IEEE Transactions on Automatic Control, vol. 54, no. 1, p. 48, 2009.

[11] M. Zhu and S. Mart´ınez, “Discrete-time dynamic average consensus,”

Automatica, vol. 46, no. 2, pp. 322–329, 2010.

[12] P. Di Lorenzo and G. Scutari, “Next: In-network nonconvex optimiza- tion,” IEEE Transactions on Signal and Information Processing over Networks, vol. 2, no. 2, pp. 120–136, 2016.

[13] S. Aeron, V. Saligrama, and D. A. Castanon, “Efficient sensor man- agement policies for distributed target tracking in multihop sensor networks,”IEEE Transactions on Signal Processing, vol. 56, no. 6, pp.

2562–2574, 2008.

[14] A. H. Sayed, S.-Y. Tu, J. Chen, X. Zhao, and Z. J. Towfic, “Diffusion strategies for adaptation and learning over networks: An examination of distributed strategies and network behavior,”IEEE Signal Process.

Mag., vol. 30, pp. 155–171, 2013.

[15] S. Xie and L. Guo, “Analysis of normalized least mean squares-based consensus adaptive filters under a general information condition,”SIAM J. Control Optim., vol. 56, pp. 3404–3431, 2018.

[16] J. C. Spall, “Multivariate stochastic approximation using a simultaneous perturbation gradient approximation,”IEEE Transactions on Automatic Control, vol. 37, no. 3, pp. 332–341, 1992.

[17] B. Polyak, Introduction to optimization. Optimization Software, Publications Division (New York), 1987.

[18] H. J. Kushner and G. Yin, Stochastic Approximation and Recursive Algorithms and Applications. New York, Springer–Verlag, 2003.

[19] V. S. Borkar,Stochastic approximation: a dynamical systems viewpoint.

Springer, 2009.

[20] O. Granichin and N. Amelina, “Simultaneous perturbation stochastic approximation for tracking under unknown but bounded disturbances,”

IEEE Transactions on Automatic Control, vol. 60, no. 6, pp. 1653–1658, 2015.

[21] J. Zhu and J. C. Spall, “Tracking capability of stochastic gradient algorithm with constant gain,” inDecision and Control (CDC), 2016 IEEE 55th Conference on. IEEE, 2016, pp. 4522–4527.

[22] ——, “Probabilistic bounds in tracking a discrete-time varying process,”

in2018 IEEE Conference on Decision and Control (CDC). IEEE, 2018, pp. 4849–4854.

[23] N. Amelina, A. Fradkov, Y. Jiang, and D. J. Vergados, “Approximate consensus in stochastic networks with application to load balancing,”

IEEE Transactions on Information Theory, vol. 61, no. 4, pp. 1739–

1752, 2015.

[24] V. Erofeeva, O. Granichin, N. Amelina, Y. Ivanskiy, and Y. Jiang,

“Distributed tracking via simultaneous perturbation stochastic approximation-based consensus algorithm,” in 2019 IEEE 58th Conference on Decision and Control (CDC). IEEE, 2019, pp.

6050–6055.

[25] S. Boyd, A. Ghosh, B. Prabhakar, and D. Shah, “Randomized gossip algorithms,”IEEE/ACM Transactions on Networking (TON), vol. 14, no. SI, pp. 2508–2530, 2006.

[26] H.-F. Chen, T. E. Duncan, and B. Pasik-Duncan, “A kiefer-wolfowitz al- gorithm with randomized differences,”IEEE Transactions on Automatic Control, vol. 44, no. 3, pp. 442–453, 1999.

APPENDIX

Proof of Theorem 1. At first, we prove that the averaged graph GB_av corresponds the average cost constraint (1).

By virtue of Assumption 5c and Cauchy-Bunyakovsky- Schwarz inequality, we get

ECost(Bt) =E X

i∈N

b^i,j_t ≤ s

nE X

i∈N

(b^i,j_t )²≤ r

nQ² n =Q.

Hence, the average cost constraint (1) holds.

At second, we study the asymptotic mean square properties of residualsη_k=kθ¯_2k−1_n⊗θ_2kk.

Denote¯s_k= _β^α

k((¯y_2k−y¯_2k−1)⊗I_d) ¯∆_k,dⁱ_t=θbⁱ

2d^t−1₂ e−θ_t, d¯t = col{d¹_t, . . . ,dⁿ_t}, where d·e is a ceiling function,

¯

wt = col{w^1,1_t ,w^2,1_t , . . . ,w^n,1_t , w^1,2_t , . . . ,w^n,n_t }, v¯t = col{˜v_t¹, . . . ,˜v_tⁿ},

Let F˜¯_k−1 = σ{F_k−1,¯v_2k−1,v¯2k, ξ_2k−1, ξ2k,∆¯k, B_2k−1} be the σ-algebra of probabilistic events generated by Fk−1,¯v2k−1,v¯2k, ξ2k−1, ξ2k,∆¯k, B2k−1, and F¯k−1 = σ{Fk−1,v¯2k−1,v¯2k, ξ2k−1, ξ2k,∆¯k}, F˜k−1=σ{Fk−1,v¯2k−1,v¯2k, ξ2k−1, ξ2k},

Fk−1⊂F˜k−1⊂F¯k−1⊂F˜¯k−1⊂ Fk.

By virtue of communication model (4), we obtain θ˜¯_t = 1_n⊗θ¯_t+ ¯w_t and, according to the algorithm (9), we have ηk=

kθ¯2k−2−1n⊗θ2k−¯sk−αγL¯2k−1(1n⊗θ¯2k−2+ ¯wt)k=

=k¯gk−¯sk+αγL¯_2k−1w¯_2k−1k

whereg¯_k= (I_nd−αγL(B_2k−1)⊗I_d)¯d_2k−2+1_n⊗(θ_2k−2− θ_2k) since it is not so hard to prove that (L(B_2k−1) ⊗ I_d)1_n⊗θ_2k−2= 0based on the properties of Laplasian matrix L(B_2k−1). Taking the conditional expectation over σ-algebra

˜¯

Fk−1, by virtue of Assumption6, we derive EF˜¯k−1η²_k=k¯gk−¯skk²+α²γ²kB2k−1k²σ²_w sinceEF˜¯k−1w¯2k−1= 0.

Assumption 5c gives the bound: EF^¯k−1kB_2k−1k² ≤ Q². Taking the conditional expectation over σ-algebra F¯_k−1, by virtue of Assumption5b, we get

EF^¯k−1η_k²=kg¯k−¯skk²+α²γ²(Q²σ_w² +σ²_Bη_k−1² ) sinceEF^¯k−1(L(B2k−1)− L(Bav))¯d_2k−2= 0.

So, we obtainEF^¯k−1η_k²=

k¯gkk²+k¯skk²−2h¯gk,¯ski+α²γ²(Q²σ²_w+σ_B²η²_k−1). (14) By virtue of Assumption 8c,d we have EF˜k−1˜vkKⁱ_k(∆ⁱ_k) = EF^˜k−1v˜_kEF^˜k−1Kⁱ_k(∆ⁱ_k) =EF^˜k−1˜v_k·0 = 0.Hence, taking the conditional expectation over σ-algebraF˜k−1 of both sides of the (14) and using observation model (2), we can assert the bound forEF^˜k−1η²_k as follows:

EF^˜k−1η_k²≤EF^˜k−1k¯g_kk²−2 α β_k

X

i∈N

hdⁱ_2k,EF^˜k−1

f˜_kⁱKⁱ_k(∆ⁱ_k)i+

+2α βk

X

i∈N

hαγ(L(Bav)dⁱ_2k−2,EF^˜k−1

f˜_kⁱKⁱ_k(∆ⁱ_k)i+

+α² β_k²

X

i∈N

EF^˜k−1

˜

v_kⁱ + ˜f_kⁱ2

kKⁱ_k(∆ⁱ_k)k²+

α²γ²(Q²σ_w² +σ²_Bη_k−1² ) (15) wheref˜_kⁱ =f_ξⁱ

2k(x2k)−f_ξⁱ

2k−1(x_2k−1).

Under fulfilment of Assumption 5d, we haveλ¯2 >0 (see [1]). Hence, for the first term in (14) we derive

EF˜k−1k¯gkk²≤d^T_2k−2(Ind−αγ(L(Bav)⊗Id))^T× (I_nd−αγ(L(Bav)⊗I_d))d_2k−2+EF^˜k−12αγ× d^T_2k−2(Ind−αγ(L(Bav)⊗Id))^T1n⊗(θ2k−2−θ2k)+

k1n⊗(θ_2k−2−θ2k)k²≤η²_k−1−d^T_2k−2αγ×

(L(Bav)⊗Id)^Td_2k−2−d^T_2k−2αγ(L(Bav)⊗Id)d_2k−2+ α²γ²d^T_2k−2(Ind− L(Bav)⊗Id)^T(L(Bav)⊗Id)d_2k−2+ EF˜k−12αγη_k−1√

nk(Id−αγ(L(Bav)kkθ_2k−2−θ2kk+ 4nδ_θ²

≤(1−2αγλ¯2+α²γ²¯λ²_m)η²_k−1+ 4αγ√

nλ¯mδθηk−1+ 4nδ²_θ. (16) For any x,z ∈ R^d, by virtue of Taylor representation of f_ξⁱ

t(x)for t^± = 2k−¹₂±¹₂, we have f_ξⁱ

t±(x) =f_ξⁱ

t±(z) +h∇f_ξⁱ

t±(z+ρ^±_ξ

t±(x−z)),x−zi, (17) whereρ^±_ξ

t± ∈(0,1).

For differencef˜_kⁱ, adding and subtractingh∇f_ξⁱ

t±(z),xⁱ_t±− zi, we derive:

f˜_kⁱ =X

t^±

±f_tⁱ±(z)±h∇f_ξⁱ

t±(z),xⁱ_t±−zi ±M¯_tⁱ±(z) (18)