# Distributed Function Calculation over Noisy Networks.

1. Introduction

Consider a network where each node i has the knowledge of its initial value [x.sub.0][i] [member of] R. The objective for each node is to compute some function relating to the initial values of all nodes, that is, [x.sub.0] , [x.sub.0] , ..., [x.sub.0] [n], under the constraint that each node transmits its value to its neighbors and updates its own value with some protocol at each time-step. Such problem is called distributed function calculation (DFC) and it has received considerable attention in recent years [1,2] due to the development of distributed networks and systems.

Various protocols are designed to achieve such a goal, for example, the nearest-neighbor rule . Moreover, it leads to a well-studied field in control theory, distributed consensus [3-5] which can be viewed as a special case of this DFC problem; some model predictive control (MPC) methods are proposed to solve the consensus problem [6-8]. When the observations of each node are noise-free,  shows that, instead of requiring asymptotic convergence, finite observations are enough to compute a predefined function. Reference  shows that consensus function calculation can be performed in a minimal number of time-steps, and, furthermore, the minimal number of time-steps can be fully characterized algebraically and graphically .

When it comes to the noisy network where each node only obtains a noisy (or uncertain) measurement of its neighbors' values (e.g., the transmission channel between nodes is noisy (bounded noise) and quantized (given a certain quantization scheme)), the estimation on initial value becomes a very challenging problem. It has been shown in  that each node can obtain an unbiased estimate of any desired linear function as a linear combination of the noisy values received from its neighbors, along with its own values, using only the first-order statistics of the noise. If the second order statistics of the noise are also known, the variance of the estimation error can be minimized distributedly by refining its estimate of the linear function. Quantization and other communication issues are addressed and analyzed in .

Instead of stimulating the network with different initial conditions  to characterize statistical properties, the focus of this paper is on what the best estimation is on the initial condition of the network given a number of successive observations. The main contribution (Section 3) is to propose an algorithm to solve initial values of all nodes with uncertain output measurements, for example, outputs with bounded channel noises and quantization, over a finite (minimal) estimation horizon. Minimal number of time-steps is characterized to fully recover the initial states of all nodes from merely the information of a randomly chosen node and its neighbors.

Furthermore, when minimal number of time-steps is not satisfied, we propose an approach to solve the upper and lower bounds on the initial values of all the nodes. These bounds are the tightest possible given the system model, output measurements, and the rough bounds on the unknown initial state. After obtaining the exact value or the tightest bounds, any exact value or bounds of objective function related to [x.sub.0] can be therefore computed. Simulations about undirected and directed networks are given in Section 4 to verify the theoretical results. Finally, Section 5 summarizes this paper.

The notation is fairly standard. For a matrix A [member of] [R.sup.MxN], A[i, j] [member of] R denotes the element in the ith row and jth column. For a column vector [alpha] [member of] [R.sup.N], [alpha][i] denotes its ith element. Similarly for a row vector [beta] [member of] [R.sup.1xN], [beta][i] denotes its ith element. We denote [e.sup.T.sub.r] = [0, ..., 0, [1.sub.rth], 0, ..., 0] [member of] [R.sup.1xN]. For vectors x, y [member of] [R.sup.n], x < y and x [less than or equal to] y denote element-wise inequalities. [I.sub.N] denotes the identity matrix with dimension N.

2. Problem Formulation

The underlying weighted direct graph of the considered network is denoted by G = (V, E, W), where V = {[v.sub.1], ..., [v.sub.n]} is the set of nodes, E [subset] V x V is the set of edges, and the adjacency matrix W = [{W[i, j]}.sub.i,j=1, ..., n] [member of] [N.sup.nxn.sub.[greater than or equal to]0], with nonnegative element W[i, j] = 1 when there is a link from i to j and equaling zero when there is no link from i to j. Let x[i] [member of] R denote the state of node i, which might represent a physical quantity such as attitude, position, temperature, voltage, and so forth.

Considering the classical protocol, the dynamics of a network of discrete-time integrator agents is defined by

[x.sub.i] [k + 1] = [n.summation over (j=1)] [d.sub.ij] [k] [x.sub.j] [k], i = 1, ..., n, (1)

where k [member of] 0, 1, ... is the discrete-time index; [d.sub.ij] [k] is the (i, j) entry of a row-stochastic matrix D = [[d.sub.ij]] [member of] [R.sup.nxn] at the discrete-time index k with the additional assumption that [d.sub.ij] [k] > 0 for all i = 1, ..., n and [d.sub.ij] [k] > 0, [for all]i [not equal to] j, if (j, i) [member of] [[epsilon].sub.n] and [d.sub.ij] [k] = 0 otherwise. Intuitively, the information state of each vehicle is updated as the weighted average of its current state and the current states of its neighbors. Equation (1) can be written in matrix form as

[x.sub.k+1] = ([D.sub.n] [k] [cross product] [I.sub.m]) x [k]. (2)

Definition 1 (DFC for the ith node ). Let [f.sub.i] be a function of the initial values; that is, [f.sub.i] ([x.sub.0] , [x.sub.0], ..., [x.sub.0] [n]) is calculable by node i if it can be calculated by node i after running iteration (2) for a sufficiently large number of time-steps.

Remark 2. In particular, we are interested in distributed linear function calculation; that is, f is linear. The distributed average consensus problem  can be categorized into this.

Definition 3 (distributed asymptotic consensus ). System (2) is said to asymptotically achieve distributed consensus if [lim.sub.k[right arrow][infinity]] [x.sub.k] = 1[c.sup.T] [x.sub.0], where 1 is the column vector with all components equal to 1 and [c.sup.T] is some constant row vector. In other words, the values of all nodes converge to the same linear combination of the initial node values [c.sup.T] [x.sub.0].

The set of all computable linear [f.sub.i] is characterized and the algorithm for each node i to compute its correspondent [f.sub.i] is proposed in  for noiseless case. However, there could be noise, uncertainty caused by quantizations, both at the process and at the observation. By taking uncertainty into account, the system can be generalized and modeled as a discrete-time network system of the form

[x.sub.k+1] = A[x.sub.k] + [B.sub.d] [d.sub.k], [y.sub.k] = C[x.sub.k] + [D.sub.d] [d.sub.k], (3)

for k [member of] N, where N := {0, 1, ..., N - 1} is the estimation horizon and [mathematical expression not reproducible] are the state, output, and noise vectors (the word "noise" here has broader meanings, e.g., noise, quantization error, and uncertainty), respectively. [mathematical expression not reproducible] are the interconnection matrix and the transmission channel noise and quantization distribution matrix, respectively, while [mathematical expression not reproducible] are the output channel disturbance distribution matrices. C is ([deg.sub.i] + 1) x n matrix (where deg; is the degree of node i) with a single 1 at each row denoting the information for the node i: its own and neighbors' values at each time-step.

3. Solution to Robust DFC Problem

3.1. Output Model. Assuming that each node in the network has some memory and computational ability, this section proposes a new scheme for each node to solve the DFC problem subject to bounded channel noises by using a system model in (3) together with output measurements over a finite estimation horizon. Firstly, a necessary and sufficient condition for each node obtaining the exact initial state value is derived. It is shown that when this condition is satisfied, the initial state value can be easily solved by doing matrix calculations using the local observations of each node.

Then, for the case in which the exact calculation of initial state value is not possible (i.e., the necessary and sufficient condition does not hold), an approach to estimate the initial condition of network is proposed. This approach computes upper and lower bounds on the initial state value by using Linear Programming (LP) optimization techniques and these bounds are the tightest possible given the network model, available output measurements, and the a priori rough bounds on the initial state and noises.

The first problem of exact initial state calculation is formulated as follows.

Problem 1. With the system dynamics defined as in (3), assume that [y.sub.0], [y.sub.1], ..., are given and

(A1) (A, C) is observable;

(A2) [n.sub.y] > [n.sub.d].

Find the sufficient and necessary condition such that the exact value of the initial condition [x.sub.0] can be solved for any d.

Here, assume that the successive observations [y.sub.k] are available for, say node i (we ignore the superscript i for notational convenience), all k [member of] N. For the estimation horizon N, let

[mathematical expression not reproducible], (4)

where [mathematical expression not reproducible]. Using an iterative computation, the system algebraic formulation is described as

[mathematical expression not reproducible]. (5)

For the reason for simplicity, the output over the estimation horizon N is represented by

y = [C.sub.y0] [x.sub.0] + [D.sub.yd] d. (6)

3.2. Exact Solution to DFC. We will have Theorem 4 to solve this problem.

Theorem 4. Let all data and assumptions be as given in Problem 1; then the exact value of the initial state [x.sub.0] can be solved from (6) for any disturbance d if and only if there exists [mathematical expression not reproducible] which spans the null space of 21yd such that

rank [[U.sup.T] [C.sub.y0]] = n. (7)

Proof. It is seen from (6) that [x.sub.0] can be solved for all d if and only if there always exists [mathematical expression not reproducible] such that

[mathematical expression not reproducible]. (8)

Note that all M in (8) satisfying [mathematical expression not reproducible] are given by M = Z[U.sup.T], where [mathematical expression not reproducible] is a free variable and UT spans the null space of [D.sub.yd]. Then condition (8) is satisfied if and only if there always exists Z such that

M[C.sub.y0] = Z[U.sup.T] [C.sub.y0] = In (9)

is satisfied. Note that there always exists Z which is the left inverse of [U.sup.T] [C.sub.y0] and can be solved from (9) if and only if [U.sup.T] [C.sub.y0] has full column rank, so that [x.sub.0] = Z[U.sup.T] y, where Z and [U.sup.T] can be obtained by using singular value decomposition .

Remark 5. The idea of Theorem 4 is that even if the noise/ quantization is unknown, we can still cancel the noise/quantization effect by using simple matrix calculations such that finding M satisfies (8).

Remark 6. It is noted that [mathematical expression not reproducible] has rank m which increases as the horizon length N increases and that [C.sub.y0] has full column rank n as long as N is larger than the observability index of the pair (A, C). Therefore, condition (7) is satisfied if and only if N is large enough such that both [U.sup.T] and [C.sub.y0] have rank larger than or equal to n.

Corollary 7. Let all data and assumptions be as given in Theorem 4. Assume that [D.sub.d] has full column rank and take v as the observability index of the pair (A, C); then the exact value of the initial state [x.sub.0] can be solved from (6) for any disturbance d if and only if the measurement horizon length is

N [greater than or equal to] max {v, n/[n.sub.y] - [n.sub.d]}. (10)

Proof. [D.sub.d] has full column rank, which implies [D.sub.yd] always has full column rank. So there exists [mathematical expression not reproducible] with rank [N.sub.y] - [N.sub.d] which spans the null space of [D.sup.T.sub.yd]. Note that rank [[U.sup.T] [C.sub.y0]] = min {([N.sub.y] - [N.sub.d]), rank [(.sub.Cy0])}; hence (7) is satisfied if and only if [N.sub.y] - [N.sub.d] = N x ([n.sub.y] - [n.sub.d]) [greater than or equal to] n and N [greater than or equal to] v.

Theorem 4 gives the sufficient and necessary condition which ensures the exact value of [x.sub.0] can be calculated. However, it may not be always implementable due to the limitation of the number of successive observations in practical network systems. Hence, it is useful to investigate the techniques for estimating initial state when output measurements are not sufficient. The rest of this section studies the case when condition (7) is not satisfied. A state estimation approach is proposed by which the initial state [x.sub.0] is estimated by solving its tightest upper and lower bounds. LP optimization techniques are used to obtain the bounds and the cases for upper and lower bounds which are dealt with separately.

3.3. Optimal Bounds for DFC. Here, assume that the rough upper and lower bounds [[x.bar].sub.r], [[bar.x].sub.r] and [[d.bar].sub.k], [[bar.d].sub.k] on the initial state and disturbances, respectively, are available and have the form

[[x.bar].sub.r] [less than or equal to] [x.sub.0] [less than or equal to] [[bar.x].sub.r], [[d.bar].sub.k] [less than or equal to] [d.sub.k] [less than or equal to] [[bar.d].sub.k].

Let [x.sub.0] [i] represent the ith element of the initial state values [x.sub.0] and [[bar.x].sub.0i]; [[x.bar].sub.0i] stand for the tightest upper and lower bounds on [x.sub.0i], where i = 1, 2, ..., n, given the system model, output measurements over the estimation horizon, and the a priori rough bounds on [x.sub.0] and d; then the problem formulation for obtaining tightest upper and lower bounds on [x.sub.0] is described as follows.

Problem 2. Let the system dynamics be as defined in (3). Assume that the rough upper and lower bounds on [x.sub.0] and d are known and that [y.sub.0], [y.sub.1], ..., [y.sub.N-1] are given. Find the upper and lower bounds [[bar.x].sub.0i] and [[x.bar].sub.0i], for i = 1, 2, ..., n, defined by

[mathematical expression not reproducible]. (12)

The following theorem shows that evaluating the upper bound on each [x.sub.0] [i] can be obtained by solving an LP problem. Note that we use boldface to denote variables in the optimization.

Theorem 8. Let all data be as defined previously and [R.sup.m.sub.+] := {v [member of] [R.sup.m] : v [greater than or equal to] 0}. Then [x.sub.0] [i] [less than or equal to] [[bar.x].sub.0i] for i = 1, 2, ..., n, if and only if there exist [mathematical expression not reproducible] such that

[mathematical expression not reproducible]. (13)

Proof. As it can be seen that [x.sub.0] [i] = [e.sup.T.sub.i] [x.sub.0], then a manipulation verifies that the following identity

[mathematical expression not reproducible]. (14)

is valid for all [mathematical expression not reproducible], where

[mathematical expression not reproducible], (15)

where * denotes terms readily inferred from symmetry.

It follows that a sufficient condition for [x.sub.0] [i] [less than or equal to] [[bar.x].sub.0i] subject to [[x.bar].sub.r] [less than or equal to] [x.sub.0] [less than or equal to] [[bar.x].sub.r], [d.bar] [less than or equal to] d [less than or equal to] [bar.d], and y = [C.sub.y0] [x.sub.0] + [D.sub.yd] d is

[mathematical expression not reproducible]. (16)

That the condition is necessary follows from Farkas' Lemma . Note that evaluating [[bar.x].sub.0i] subject to (16) is equivalent to solving the LP problem as follows:

[mathematical expression not reproducible]. (17)

Therefore, the upper bound on each [x.sub.0] [i] for i = 1, 2, ..., n can be computed by solving (17).

Similarly, the next theorem is proposed to solve the optimal lower bound on each [x.sub.0] [i] for i = 1, 2, ..., n. Since the proof can be obtained by following that of Theorem 8, it is omitted.

Theorem 9. Let all data be as defined above and [R.sup.m.sub.-] := {v [member of] [R.sup.m] : v [less than or equal to] 0}. Then [[x.bar].sub.0i] [less than or equal to] [x.sub.0] [i] for i = 1, 2, ..., n if and only if there exist [mathematical expression not reproducible] such that L([[x.bar].sub.0i], [[bar.[[mu]]].sub.i], [[[mu].bar].sub.i], [[bar.[tau]].sub.i], [[[tau].bar].sub.i], [[lambda].sub.i]) [??] 0, where L(x, x, x, x, x, x) is defined in (16). Similarly to (17), evaluating the tightest lower bound on [x.sub.0] [i] is equivalent to solving the following LP problem:

[mathematical expression not reproducible]. (18)

Remark 10. From Theorems 8 and 9, it can be seen that sufficient and necessary conditions for [[bar.x].sub.0i] and [[x.bar].sub.0i] to be upper and lower bounds on [x.sub.0] [i] are given in the form of the LP conditions (17) and (18), respectively. Then the tightest bounds on the initial state are obtained by solving 2 x n LP problems.

Remark 11. Since LP problems can be easily solved, another benefit of the proposed method is that we can use it as a recursive online implementation; that is, at each time-step k, node i uses Theorems 8 and 9 to compute the upper and lower bounds on the initial state [x.sub.0] until condition (7) is satisfied, in which case the initial conditions can be fully recovered.

4. Simulation

4.1. Directed Network Example. Next, we consider a directed network with 14 nodes in Figure 1. Its corresponding adjacency matrix W and the underlying graphical representation are shown as follows:

[mathematical expression not reproducible]. (19)

Take the sampling time to be e = 1/6; then the system matrices are obtained as follows:

[mathematical expression not reproducible]. (20)

We randomly choose an initial state

[mathematical expression not reproducible]. (21)

to stimulate the system. It is assumed that there are 6 transmission noises between nodes which are zero-mean bounded from [-1, 1]. Also, assume that the rough upper and lower bounds on the initial state [[bar.x].sub.r] and [[x.bar].sub.r] are 1 and -1, respectively. From (7), it can be seen that the minimal number of steps N for direct calculation on [x.sub.0] is 7.

Therefore, in the scenario that there are enough output measurements over the horizon length N = 7, condition (7) is satisfied. So there exists M satisfying (8) such that [x.sub.0] = My. By taking the matrix calculation, the initial state [x.sub.0] is recovered.

For the scenario that there are only output measurements over the horizon length N = 6, the full recovery of [x.sub.0] is not possible. By using the method given in Theorems 8 and 9, the tightest upper and lower bounds on [x.sub.0] are solved and given as follows:

[mathematical expression not reproducible]. (22)

As it can be seen that the difference between the upper and lower bounds is small, [x.sub.0] is well-estimated.

By using the exact value or the tightest bounds for [x.sub.0], any exact value or bounds of objective function related to [x.sub.0] can be therefore obtained. For example, here, the true average consensus value of this system, that is, 1/14 ([[summation].sup.14.sub.i=1] [x.sub.0] [i]), is -0.0933, while the bounds on it obtained from our method are [-0.1768,-0.0083] when N = 6 which are very tight with respect to the noise level [+ or -] 1.

Figure 2 illustrates the bounds on the average consensus value as the horizon of output measurements reaches 7 in which the exact average consensus value is solved. It is shown that the bounds become tighter and tighter as more and more measurement information becomes available.

5. Conclusion

We have studied the problem of distributed function calculation in a connected network operating with noise. In particular, this is not focused on the nearest-neighbor protocol, where each node updates its own value at each step as a linear combination of its own previous values and those of its neighbors. Each node tries to compute a linear function of all the initial states of all nodes. Algorithm to compute the function in question is proposed when the number of successive observations is large enough; we also propose an approach to obtain tightest bounds on such a function based on the available observations while the number of observations is not large enough.

http://dx.doi.org/10.1155/2016/6093293

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

This work was supported by National Natural Science Foundation of China (NNSFC) under Grant no. 61322304.

References

 A. Giridhar and P. R. Kumar, "Computing and communicating functions over sensor networks," IEEE Journal on Selected Areas in Communications, vol. 23, no. 4, pp. 755-764, 2005.

 S. Sundaram and C. N. Hadjicostis, "Distributed function calculation and consensus using linear iterative strategies," IEEE Journal on Selected Areas in Communications, vol. 26, no. 4, pp. 650-660, 2008.

 A. Jadbabaie, J. Lin, and A. S. Morse, "Coordination of groups of mobile autonomous agents using nearest neighbor rules," IEEE Transactions on Automatic Control, vol. 48, no. 6, pp. 988-1001, 2003.

 W. Ren and R. W. Beard, "Consensus seeking in multiagent systems under dynamically changing interaction topologies," IEEE Transactions on Automatic Control, vol. 50, no. 5, pp. 655-661, 2005.

 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.

 Z. Cheng, H.-T. Zhang, M.-C. Fan, and G. Chen, "Distributed consensus of multi-agent systems with input constraints: a model predictive control approach," IEEE Transactions on Circuits and Systems. I. Regular Papers, vol. 62, no. 3, pp. 825-834, 2015.

 Z. M. Cheng, M.-C. Fan, and H.-T. Zhang, "Distributed MPC based consensus for single-integrator multi-agent systems," ISA Transactions, vol. 58, pp. 112-120, 2015.

 H.-T. Zhang, M. Z. Q. Chen, and T. Zhou, "Improve consensus via decentralized predictive mechanisms," EPL, vol. 86, no. 4, Article ID 40011, 2009.

 Y. Yuan, G.-B. Stan, L. Shi, and J. Goncalves, "Decentralised final value theorem for discrete-time LTI systems with application to minimal-time distributed consensus," in Proceedings of the 48th IEEE Conference on Decision and Control held jointly with 2009 28th Chinese Control Conference (CDC/CCC '09), pp. 2664-2669, Shanghai, China, December 2009.

 Y. Yuan, G.-B. Stan, M. Barahona, L. Shi, and J. Goncalves, "Decentralised minimum-time consensus," Automatica, vol. 49, no. 5, pp. 1227-1235, 2013.

 S. Sundaram and C. N. Hadjicostis, "Distributed calculation of linear functions in noisy networks via linear iterations," in Proceedings of the 47th IEEE Conference on Decision and Control (CDC '08), pp. 5462-5467, Cancun, Mexico, December 2008.

 R. Carli, F. Fagnani, A. Speranzon, and S. Zampieri, "Communication constraints in the average consensus problem," Automatica, vol. 44, no. 3, pp. 671-684, 2008.

 L. Xiao and S. Boyd, "Fast linear iterations for distributed averaging," Systems & Control Letters, vol. 53, no. 1, pp. 65-78, 2004.

 Z. Zhang and I. M. Jaimoukha, "An optimal solution to an H_/[H.sub.[infinity]] fault detection problem," in Proceedings of the 50th IEEE Conference on Decision and Control and European Control Conference (CDC-ECC '11), pp. 903-908, Orlando, Fla, USA, December 2011.

 R. Jagannathan and S. Schaible, "Duality in generalized fractional programming via Farkas' lemma," Journal of Optimization Theory and Applications, vol. 41, no. 3, pp. 417-424,1983.

Zhidun Zeng, (1) Xin Yang, (2) Ze Zhang, (3) Xiaoyu Mo, (4) and Zhiqiang Long (2)

(1) College of Physics and Optoelectronics, Taiyuan University of Technology, Taiyuan 030600, China

(2) College of Mechatronics Engineering and Automation, National University of Defense Technology, Changsha 410000, China

(3) Expedia, Angel Building, St. John Street, London EC1V4EX, UK

(4) School of Automation, Key Laboratory of Image Processing and Intelligent Control, Huazhong University of Science and Technology, Wuhan 430074, China

Correspondence should be addressed to Ze Zhang; zzhang1@expedia.com

Received 21 December 2015; Revised 23 March 2016; Accepted 23 March 2016