Distributed Aggregate Function Estimation by Biphasically Configured Metropolis-Hasting Weight Model

RADIOENGINEERING, VOL. 26, NO. 2, JUNE 2017 479 DOI: 10.13164/re.2017.0479 APPLICATIONS OF WIRELESS COMMUNICATIONS Distributed Aggregate Function Estimation by Biphasically Configured Metropolis-Hasting Weight Model Martin KENYERES 1, Jozef KENYERES 2, Vladislav SKORPIL1, Radim BURGET1 1 Dept. of Telecommunications, Brno University of Technology, Technická 12, 612 00 Brno, Czech Republic 2 Sipwise GmbH, Europaring F15, 2345 Brunn am Gebirge, Austria [email protected], [email protected]m, [email protected], [email protected] Submitted December 12, 2016 / Accepted April 7, 2017 Abstract. An energy-efficient estimation of an aggregate function can significantly optimize a global event detection or monitoring in wireless sensor networks. This is probably the main reason why an optimization of the complementary consensus algorithms is one of the key challenges of the lifetime extension of the wireless sensor networks on which the attention of many scientists is paid. In this paper, we introduce an optimized weight model for the average consensus algorithm. It is called the Biphasically configured Metropolis-Hasting weight model and is based on a modification of the Metropolis-Hasting weight model by rephrasing the initial configuration into two parts. The first one is the default configuration of the Metropolis-Hasting weight model, while, the other one is based on a recalculation of the weights allocated to the adjacent nodes’ incoming values at the cost of decreasing the value of the weights of the inner states. The whole initial configuration is executed in a fully-distributed manner. In the experimental section, it is proven that our optimized weight model significantly optimizes the Metropolis-Hasting weight model in several aspects and achieves better results compared with other concurrent weight models. Keywords Distributed computing, aggregate function, average consensus algorithm, Metropolis-Hasting weight model, wireless sensor networks 1. Introduction 1.1 Wireless Sensor Networks Wireless sensor networks (WSNs) are systems intended to perform a real-time detection of a stochastic event or to monitor physical quantities [1–2]. They are formed by battery-constrained nodes deployed in a geographical area where a phenomenon of interest is observed. These nodes are equipped with hardware components such as a wireless transceiver, a sensor to sense physical quantities, a central processor unit, a source of energy etc. [3–4] Thus, the nodes are able to obtain necessary information about the observed phenomenon, process it, mutually exchange data and make a meaningful decision on the examined physical quantity [4]. Due to their character, the WSNs find the application in various areas such as military surveillance, a natural disaster detection and its elimination, habitat monitoring, inventory tracking, an acoustic detection, pollution monitoring, medical systems, target tracking, a robotic exploration, a health care (especially, they find the usage in the scenarios considering monitoring elderly patients in a remote area), environment monitoring, a micro surgery, agriculture etc. [5–6]. In many applications, WSNs may be formed by hundreds of nodes potentially situated in inaccessible locations and therefore, a battery recharge or replacement may be complicated [7]. An exhausted battery results in a node death, which can decrease the quality of the final decisions or even prevent a whole system from fulfilling its functionality. As the results, the attention of many scientists has been focused on an optimization of the energy consumption aspect in the last years [8–11]. It is because an effective optimization can significantly increase the network lifetime of a WSN application [12]. In [13], the authors divide the architectures of a global event detection into three categories. The second and the third architecture require a complementary consensus algorithm to estimate aggregate functions in order to ensure a higher credibility of the measured outputs. These architectures do not assume the presence of a fusion center in a network. The implementation of this supplementary algorithm ensures a higher precision of the final decision on the observed phenomenon in many applications [14]. A decision made according to data obtained by independently-measuring nodes secures a more credible output than a decision made in terms of a single measurement and minimizes a change of an incorrect classification [14]. The importance of the consensus algorithm implementation for high-quality monitoring in WSNs is discussed in [15]. 480 M. KENYERES, J. KENYERES, V. SKORPIL, DISTRIBUTED AGGREGATE FUNCTION ESTIMATION BY BIPHASICALLY … 1.2 Average Consensus Algorithm Due to the character of the WSNs, the modern applications are often based on the implementation of distributed mechanisms. The algorithms of distributed computing substitute the older frequently-implemented centralized manner of the computation [14]. Despite its reliability and high precision, the centralized algorithms do not pose the optimal solution for the implementation into the systems formed by battery-constrained devices. One of the most appropriate distributed algorithms for WSN applications is average consensus, which is a fully distributed iterative algorithm primarily for estimation the average from the values of all the nodes present in a network [16]. This algorithm does not require the presence of any fusion center. The nodes are able to estimate the average by a mutual exchange of their inner states with the nodes situated in the adjacent area. The average consensus algorithm is characterized by a high flexibility because its execution is modifiable by the chosen weight model [17]. The weight models differ from each other in several aspects, for example, we can list the convergence rate of the algorithm, the process of the initial configuration, the information that is necessary for its proper functionality, the robustness etc. [17], [18]. In this paper, we focus on the Metropolis-Hasting model, which requires only the locallyavailable information for the initial configuration and therefore, it is one of the most preferred solutions for a real implementation into battery-constrained systems [19]. As mentioned earlier, the average consensus algorithm is primarily proposed for the estimation of the average value. However, tiny modifications can ensure that the algorithm is able to estimate other aggregate functions. One of the other frequently-used applications is the estimation of the network size. The information about the number of the nodes in a network is crucial for a proper functionality of many distributed systems [20]. In this case, the execution of the average consensus algorithm is modified in such a way that one of the nodes has the initial value set to the value equaling 1 (it is called the leader) [21]. The other nodes are set to 0. Subsequently, the nodes converge to the value equaling the reciprocal of the size of a network [21]. However, this modification causes several problems. One of the most significant ones is how to appoint the most suitable node as the leader. It often requires the implementation of other complementary mechanisms to determine this, which is not the optimal solution for battery-constrained devices [21]. As shown later in this paper, a bad choice of the leader can significantly decrease the convergence rate of the algorithm. Thus, an improvement of the leader selection can significantly optimize WSN applications by removing the necessity for other complementary algorithms. 1.3 Motivation All the previously-discussed problems motive us to propose an optimization mechanism of the MetropolisHasting weight model that improves the convergence rates, the number of the necessary messages and minimize the negative effect caused by an inappropriate choice of the leader. An optimization of these mentioned aspects can significantly accelerate and simplify the computation process and therefore optimizes the real-life applications of WSNs. The choice of this weight model is affected by an effort to improve a weight model of the average consensus algorithm that finds the wide usage thanks to its specific character. The Metropolis-Hasting weight model fulfills these criteria because it does not require any global information about the network for its proper initial configuration and so it works in a fully distributed manner. Thus, it is an appropriate solution for an implementation into the WSNs. Thus, this was the main reason that motivated us to focus our research on this weight model. In this paper, we introduce an optimized weight model derived from the Metropolis-Hasting weight model that improves the discussed aspects of this model. The optimized model modifies the weight matrix of the Metropolis-Hasting weight model by an additional step during the initial configuration. Thus, the weight matrix is initially configured twice. The first configuration is the default one defined within the Metropolis-Hasting weight model and the other one poses the novelty proposed by us. The other phase is based on recalculating the weights allocated to the adjacent nodes’ incoming states at cost of decreasing the value of the weight of the inner states in a distributed way. 1.4 Paper Organization In Sec. 2, we turn our attention to the latest papers related to the average consensus optimization. In the next section, we provide the mathematical tool used to model the average consensus algorithm executed in the WSNs, present the main features of this algorithm and adduce important theorems defined within the spectral graph theory. We also introduce concurrent weight models which the optimized weight model is compared with. In Sec. 4, we introduce our optimized weight model, provide mathematical tools to model it and derive the convergence proof. In Sec. 5, we examine the optimization of chosen aspects ensured by the optimized weight model. We focus on an optimization of the average estimation, the network size estimation and the range of the convergence rates caused by the choice of the leader. We compare this optimization in three types of the networks – weakly, averagely and strongly connected. All of them consist of ten randomly generated networks. In Appendix section, we adduce the complete results obtained within our numerical experiments. 2. Related Work This section is devoted to an insight into an optimization mechanism proposed for the average consensus algo- RADIOENGINEERING, VOL. 26, NO. 2, JUNE 2017 481 rithm. We introduce the latest papers dealing with optimization mechanisms. In [19], [22–25], the authors’ attention focuses on the Metropolis-Hasting weight model (also in various applications as a complementary mechanism), which our optimized weight model is derived from. It was developed by Metropolis, Rosenbluth and Teller in 1953 and generalized by Hastings in 1970. It was originally defined within Markov chain Monte Carlo methods and proposed to simulate complex, non-standard multivariate distributions [26]. Its modification for the consensus achievement problem finds the usage in many applications due to its character. For its proper initial configuration, only locally-available information is required, i.e. the number of the neighbors of a particular node and the number of neighbors of the node from its adjacent area. This significantly simplifies the initial configuration phase. Thus, this weight model finds the usage in many applications. Additionally, it also poses a robust solution against a quantization noise [18]. In [27], an average consensus optimization based on the usage of the opportunistic inter-agent communication to achieve the consensus is presented. Each node is endowed by a local criterion determining when to broadcast the inner state to the nodes situated in the adjacent area. In paper [28], a novel consensus protocol is presented that achieves the average state consensus for multi-agent systems in finite time. The protocol contains a linear and a non-linear term. The state consensus is achieved by the non-linear term, while the performance optimization is ensured by the linear term to some degree. The authors of paper [29] present an optimization mechanism based on a division of the computation process into two phases. The first one is the phase of reaching the local consensuses and the second one is the phase of reaching the global consensus. Within the first phase, a network is reorganized into geographically close areas so-called packs. Here, each node converges to the value equaling the average of all the nodes present in a pack. Subsequently, each pack appoints one of the nodes as the head, which communicates with the other heads and converges to the average value. The authors of [30] present a novel continuous-time dynamic average consensus algorithm for networks with the interaction that can be described by weight-balanced directed graphs of a strong connectivity. The nodes are able to track the average of the dynamic inputs with some non-zero steady-error. Its size is controlled by exploiting a design parameter. In [31], a distributed algorithm for average consensus that solves the discrete-time average consensus problem on strongly connected weighted digraphs is presented. Its principle lays in the computation of the average value using the estimation of the left eigenvector associated with the zero eigenvalue of the Laplacian matrix. The authors of [32] built their optimization mechanism on the exploiting of the second-order neighbors. They focus their attention on both the continuous-time case, where the edges are chosen by solving a convex optimization problem formed by utilizing the convex relaxation method, and the discrete-time case, when the edges are chosen using the brute force method. In [33], an optimization mechanism is presented that exploits the prediction of the future value of the inner states. This technique is based on the estimation of the states for the next iterations in terms of the values of the inner states from the previous iterations. The authors of [34] introduce IACA, which is a two-layer improved consensus algorithm of a multi-agent system. The authors propose a new distributed cost optimization method for loading shedding of an islanded microgrid considering cost. The technique solves distributed cost optimization of load shedding by exploiting the synchronization processing of IACA in the layer 2. The authors of [35] present an optimization mechanism minimizing the negative effects caused by a random packet loss. It is based on keeping track of the changes in the state variable, which the neighbors influence causes. The papers [36–38] focus on the Maximum Degree weight model, which is a modification of the Constant weight model, and its applications. The Constant weight model is characterized by the parameter ε, which affects the convergence rate as well as the interval of the convergence. The higher value it takes, the faster the algorithm is. However, a too high value can cause the divergence of the algorithm. The divergence is a type of a failure when the convergence is not reached. Instead of it, the nodes diverge to infinite large values [39]. This error poses a serious problem that stunts a whole network [39]. The Maximum Degree weight model is based on the setting of the parameter ε to the value equaled to the reciprocal of the number of the neighbors of the best-connected node in a network. The initial configuration requires the knowledge about this value and therefore, it is necessary to implement a supplementary mechanism to determine it [19]. In [18], [40], the Best Constant weight model is discussed. Its optimized version is based on the utilization of the knowledge about the second smallest and the largest eigenvalue of the Laplacian matrix [18]. To compose it, it is necessary to know the information about the complete network topology. Thus, this weight model requires a particular centralization for its optimization. 3. Modeling of Average Consensus Algorithm in WSNs In Sec. 3.1, we introduce the used mathematical model of the WSNs executing the average consensus algorithm and the main features of this algorithm. In Sec. 3.2, we discuss and mathematically describe the concurrent weight models of average consensus with which our optimized weight model is compared. 3.1 Used Mathematical Model In order to model the WSNs, a mathematical tool defined within the spectral graph theory is used [41], [42]. A WSN is considered to be an indirect finite graph defined 482 M. KENYERES, J. KENYERES, V. SKORPIL, DISTRIBUTED AGGREGATE FUNCTION ESTIMATION BY BIPHASICALLY … as G = (V,E). The set V is formed by all the vertexes, which are representatives of the particular nodes. Each node is labeled by the unique identity number vi. We assume that the nodes are labeled with the numbers 1, 2,...,N, where N is the size of a network and therefore |V| = N. The mutual connectivity between the nodes is indicated by the existence of an edge. The set E⊂VxV consists of all the edges present in a graph. The edge is labeled as (vi,vj) or eij. We assume the range homogeneity of the nodes and therefore the following statement holds: EE   jiij ee . (1) There are several tools to describe a network topology within the spectral graph theory. One of them is the Laplacian matrix for a description of the mutual connectivity among the nodes. It is a square symmetric matrix defined for all the indirect graphs as follows [43]:         ,otherwise if if ,0 , ,1 ][ ji e dL ij iij E (2) Here, di is the degree of a vertex vi and so, the number of the corresponding node’s neighbors. Except for the mutual connectivity, the Laplacian matrix provides other useful information about the topology. Let us focus on the following sentence [44]: Lemma 1: Let G = ( V , E ) be a graph and let 0 = µ 1 ( L ) ≤ µ 2 ( L ) ≤ µ 3 ( L ) ≤……≤ µ N ( L ) be the ascen - dingly - ordered eigenvalues of the Laplacian matrix of this graph. Then, G is not connected if µ 2 ( L ) = 0 . According to Lemma 1, only the networks with the Laplacian matrix whose second smallest eigenvalue is not equaled to 0 are connected. In the case when µ2(L) equals 0, the average consensus algorithm does not estimate the average from the values of all the nodes but estimates the set of the local averages in each connected subpart of a graph. Therefore, we assume only topologies whose second smallest eigenvalue of the Laplacian matrix is not equaled to 0. The knowledge about the exact value of µ2(L) and µN(L) is necessary for the optimized initial configuration of the Best Constant weight model [26]. A configuration of this model with a smaller positive value ensures the convergence but the execution of the algorithm is slower [18]. As mentioned above, the average consensus algorithm is an iterative distributed algorithm based on a mutual exchange of the current states among the nodes. The algorithm is modeled by the difference equation defined as follows [45]: )()1( kk xWx   . (3) Here, W ⊂ L is a weight matrix of the algorithm and the time-variant vector x(k) ∈ RN  1 gathers all the inner states at kth iteration. We assume that the initial states are labeled as k = 1. The elements of W depend on the used weight model. This matrix also provides useful information about the network topology. The following lemma says about the connectivity of the topology [44]: Lemma 2: Let G = ( V , E ) be a graph and let 1 = λ 1 ( W ) ≥ λ 2 ( W ) ≥ λ 3 ( W ) ≥……≥ λ N ( W ) ≥ –1 be the descendingly - ordered eigenvalues of the weight matrix. Then, G is not connected if max { λ 2 ( W ) , –λ N ( W )} = 1 . Within the spectral theory, it is defined that the value max{λ2(W), –λN(W)} equals the spectral radius ρ of the matrix determined as the difference between the matrix W and matrix defined as 1/N·1×1T [46]. Thus, in terms of the previous statement, we can write as follows: T 2 1 ()max{(),()} N N    W11 W W . (4) Let us focus on the features of the weight matrix W. According to [47], the weight matrix W is required to hold the following conditions: 11W   , (5) TT 1W1  , (6) T 1 ()1 N    W11. (7) Here, the vector 1 is a column vector whose all elements are equaled to 1 (its size is implicitly assessable from the previous context). Fulfilling the formula (7) ensures the convergence of the average consensus algorithm, meanwhile, the formulae (5) and (6) imply that the weight matrix is doubly stochastic (sometimes, labeled as bistochastic) and determines the convergence point [47]. These two formulas also implicate the following statement: T WW . (8) Within our analysis, we use also another descriptive tool defined within the spectral graph theory. It is called the adjacency matrix A ∈ {0, 1}N  N and contains the information about the mutual connectivity between the pairs of the nodes. It is a diagonally symmetric matrix of a square shape for all the indirect graphs. The direct connection (i.e. the existence of an edge) is indicated by the presence of 1 in the corresponding position. Thus, the presence of 0 is an indicator that two nodes are not directly connected to one another. Mathematically, the adjacency matrix is defined as follows [48]: 1, if ( , ) [] 0, otherwise ij ij vv A    E (9) Another useful tool is the identity matrix I ∈ {0,1}N  N defined as (10) [49]. In literature, also the other notation I = diag(1,1,….1) can be found. 1, if [] 0, otherwise ij ij I    (10) Usually, a lower index is allocated to indicate its size (the label IN indicates that the underlying matrix has the size equaled to the number of the size a network). RADIOENGINEERING, VOL. 26, NO. 2, JUNE 2017 483 As mentioned earlier, the average consensus algorithm is an iterative algorithm (regardless of the used weight model) executed in such a way that the nodes converge to the value determined as the average calculated from all the initial values [50]. Therefore, this behavior can be described as follows [45]: T 1 lim ( ) lim (1) (1) k kk kN      11 xWx x . (11) Only the existence of this limit ensures the convergence of the average consensus algorithm. As mentioned, it is achieved by using such a weight matrix that holds the conditions (5-7) [46]. As the algorithm convergences to the value in the infinite [51], it is necessary to implement a mechanism indicating the consensus. We use the mechanism defined as follows:   )}(min{)}(max{ kk xx . (12) Here, the parameter δ determines the precision. Its higher values ensure a higher precision at the cost of a slower convergence rate. In our experiments, we assume that its value is equaled to 0.0001. 3.2 Concurrent Weight Models As mentioned, we compare our optimized weight model with other three concurrent ones. In order to distinguish these models from each other, the weight matrix of a particular model has an upper index with the abbreviated name of the model. We use the following abbreviations:  Metropolis-Hasting weight model – MH  Maximum Degree weight model – MD  Best Constant weight model – BC  Biphasically configured Metropolis-Hasting weight model – BMH The first examined model is the Metropolis-Hasting weight model whose weight matrix is defined as follows [19], [22-25]: 1 MH MH 1, (1 max{ , }) , if ( , ) []1 [],if otherwise 0, ij ij N ij ik kki dd vv WWij              E (13) As mentioned, our contribution optimizes this model. The second model is the Maximum Degree weight model. It is derived from the Constant weight model in such a way that the weighting parameter ε is set to the value equaling the reciprocal of the degree of the best-connected node. Therefore, it is defined as follows [36–38]: max MD max 1/ , if ( , ) []1/,if 0, otherwise ij ij i dvv Wddij         E (14) Within the initial configuration phase, this model requires each node in a network to be aware of the number of the neighbors of the best-connected node. In order to get this information in a distributed manner, it is necessary to implement a complementary algorithm [19]. The last examined model is called the Best Constant weight model. We assume its optimized variant even though it requires the information about the second smallest and the largest eigenvalue of the Laplacian matrix. It is defined as follows [40]: 2 BC 2 2 / ( ( ) ( )), if ( , ) [ ] 12 /( () ()), if 0, otherwise Nij ij i N vv Wd ij            LL E LL (15) 4. Biphasically Configured MetropolisHasting Weight Model In this section, we introduce our Biphasically configured Metropolis-Hasting weight model. In Sec. 4.1, we explain the model mechanics and provide the mathematical description. In Sec. 4.2, the convergence conditions are presented. 4.1 Principle of Biphasically Configured Metropolis-Hasting Weight Model This subsection focuses on an introduction of the main features of our optimized weight model. As mentioned above, it is called the Biphasically configured Metropolis-Hasting weight model and is derived from (as its name implies) the Metropolis-Hasting weight model. As already discussed, the Metropolis-Hasting model is appropriate for the implementation into real-life applications thanks to its simplified demands for the initial configuration. In order to correctly fulfill its functionality, each node has to be aware of the number of its neighbors as well as the number of the neighbors of the adjacent nodes. Thus, only locally available information is necessary for the correct initial configuration. There are several approaches to obtain this information (centralized one, distributed one assuming a phase when this information is distributed in the adjacent area, a manual configuration etc.). The most appropriate solution depends on a particular application. Our optimized weight model is based on rephrasing the initial configuration process of the Metropolis-Hasting weight model into two phases. The first phase is identical to the default configuration of the Metropolis-Hasting weight model. The other phase consists of a recalculation of the weights allocated to the incoming values from the adjacent nodes. Within this phase, each node determines when to do this recalculation according to its unique identity. Thus, each node has to be additionally aware of the diagonal value of the weight matrix WBMH corresponding to all its neighbors – this value is locally available. The 484 M. KENYERES, J. KENYERES, V. SKORPIL, DISTRIBUTED AGGREGATE FUNCTION ESTIMATION BY BIPHASICALLY … recalculation of WBMH has to be executed sequentially, i.e. the node with the identity number equaled to 1 initiates the whole process, updates the matrix WBMH (the active node updates its inner updating rules as well as informs its neighbors about new weights) and only then the node 2 can start the recalculation. Thus, we assume the variability of the weight matrix WBMH(a) during the second phase of the configuration process. The parameter a takes the values from 1, 2,…,N and labels the active node (i.e. the one that is allowed to make the recalculation) as well as the round of the recalculation process (we assume that WBMH (0) = WMH and each label of the round corresponds to the unique number of a node). Thus, we label the node currently allowed to make the update as va. The length of the recalculation process is determined by the size of a network (i.e. the process lasts N rounds). Let us define the set Νa gathering all the nodes from the adjacent area of va. So, we can write as follows:  1][  aiia AvΝ. (16) At the round when a node va is allowed to update WBMH, it has to be aware of the current value in the diagonal corresponding to it and all its neighbors. These values represent the weights of the current inner state. Subsequently, va calculates the growth coefficient defined as follows:       N j ajjj N j ajjjaa a AaW AaWaW v 1 BMH 1 BMHBMH ].[)]1([ }].[)]1([,)]1(min{[ )(  . (17) The choice of the minimal value from the weight of the inner state of the active node and the sum of weights of the inner states of its neighbors ensures that the growth coefficient is never greater than 1. Thus, the convergence of the weight model regardless of the underlying topology (see Sec. 4.2) is guaranteed. Subsequently, the node currently making the recalculation decreases the weight of its inner state to the minimal possible value (the ideal scenario is when this weight is equaled to 0 after finishing this procedure) and distributes this value among the neighbors in terms of the weights of the inner states of these nodes. From the central view, this procedure is described according to the following rules: BMH BMH BMH [()][(1)]()[(1)] for . ia ia a ii ia Wa Wa vWa v   Ν (18) This formula describes an increase of the weights of the incoming values of the adjacent nodes. Their values are increased with the value equaled to the diagonal value of the adjacent nodes (i.e. the weight of their inner states) weighted by the growth coefficient. The Metropolis-Hasting weight model assumes a doubly stochastic matrix (i.e. an edge eij is allocated only one weight – the incoming state of va and vi is weighted with the same value) and therefore, it is necessary to preserve this condition: BMH BMH BMH [()][(1)]()[(1)] for . ai ai a ii ia Wa Wa vWa v   Ν (19) After all the neighbors are allocated a new weight, the nodes with an increased weight have to decrease their diagonal value by the increase of their incoming value weight. Otherwise, the convergence conditions may not be preserved. From the central view, it is possible to describe the previous procedure using tools defined within the spectral graph theory as follows: BMH BMH 1, [ ( )] 1 [ ( )] , for , . N ii ij i a a jji Wa Wa v v     Ν (20) After all the nodes execute the recalculation described above (this procedure is repeated for all N nodes), the weight matrix for Biphasically configured MetropolisHasting weight model is completed and the average consensus algorithm can be executed according to (3). 4.2 Convergence Proof In the following subsection, we provide the sufficient conditions for the convergence of average consensus algorithm. The average consensus algorithm whose weights are symmetric can be described using a weighted graph [52]. The non-zero elements of its adjacency matrix AWG are allocated a strictly positive weight [AWG]ij = wij. Subsequently, it is possible to derive the weighted Laplacian matrix as [52]: WG WG WG diag{ }LdA . (21) Here, the dWG is a weighted degree vector formed by the value of the degrees of the nodes. The weighted degree vector is defined as follows [52]: 1Ad  WGWG . (22) Firstly, we show that the weights of the MetropolisHasting weight model (13) ensure the convergence of the algorithm regardless of the underlying topology. Within the spectral graph theory, it is defined [52] that W = I – LWG and therefore, W is doubly stochastic with an eigenvalue with the magnitude equaled to 1 and associated to the eigenvector with the values N–1/2. Furthermore, the matrix LWG is semidefinite and fulfills the following statement [52]: WG ( ) 1 ( ) for kk k   WL . (23) Thus, the convergence condition (7) can be reformulated as follows [52]:  2)( WG L  ,  Multiplicity of the zero eigenvalue of LWG has multiplicity one. The second statement is satisfied for all the connected graphs. It can be confirmed using the quadratic form defined as follows [52]: 2 TWG () ij ij i j e wvv     E vL v . (24) RADIOENGINEERING, VOL. 26, NO. 2, JUNE 2017 485 The vector v is an eigenvector associated to the weighted Laplacian matrix LWG. As seen, the quadratic form is equaled to the zero value if and only if vi = vj. This statement is valid for the weights of the positive values. This requirement is met by the unique normalized vector v = N–1/2. Let us focus on the first constraint. As (21) and (22) hold, the eigenvalues λWG satisfy the following condition according to Gershgorin circle theorem, which is defined as follows [53]: WG Lii dd  WGWG )(  . (25) The parameter di WG is the weighted degree of the node i defined as:   N j iji wd 1 WG . (26) The value also presents the ith row of the vector di WG. In particular, the following statement is valid [52]: WG WG ()2 i d  L. (27) Thus, di WG ≤ 1 is a sufficient condition for the average consensus algorithm to convergence for all i [52] (when a graph is neither bipartite nor regular). Since it is improbable that the graph describing a WSN is bipartite regular [52], we do not deal with these critical graph topologies. From (13), it is clear that this condition holds for each topology. Thanks to the expression in the numerator in (17), our recalculation always ensures di WG ≤ 1. It secures that the sum in both the rows and the columns does not change despite the recalculation. This guarantees the convergence conditions for our mechanism. Let us analyze the functionality of our optimized weight model. The minimal value from two parameters in the numerator of (17) is chosen to preserve the convergence conditions. The value of [WBMH(a – 1)]aa poses the maximal possible value with which va can decrease the weight of its inner state and distribute it among its adjacent nodes. A decrease with a value greater than [WBMH(a – 1)]aa always causes di WG ≤ 1 not to hold. Thus, in such a scenario, the convergence conditions are not fulfilled. Regarding all the positive values smaller than [WBMH(a – 1)]aa: the convergence conditions are preserved but the optimization is less significant. Now, let us focus on Σj[WBMH(a – 1)jj·[A]aj. There can be a scenario when va can increase the weights of the incoming values of its adjacent nodes with a value that causes some of its neighbors to achieve a negative value in the diagonal after (20). In this scenario, the growth coefficient χ is greater than 1, which results in di WG > 1. Therefore, parameter Σj[WBMH(a – 1)jj·[A]aj ensures that the growth coefficient never exceed the value 1 and so, the convergence conditions are always preserved. The procedure in (20) must be executed in order to keep the weight matrix W doubly-stochastic. This step secures that (5) and (6) always hold. 5. Numerical Experiments and Discussion In this section, we present the results of the numerical experiments executed in Matlab R2015a. All the used software was designed by the authors of this paper. In our experiments, three sets of networks are assumed with randomly generated topologies. We assume weakly, averagely and strongly connected networks. Each set consists of 10 unique topologies with the size of 200 nodes. Due to the limited range of the paper, only one representative of these sets is shown in Fig. 1, Fig. 2 and Fig. 3 respectively. The networks were generated as follows: each free position within the working area of a square shape was allocated the probability equaled to the reciprocal of the number of free positions. Thus, the choice of the position where a node was placed had a uniform distribution. Subsequently, the nodes situated in the transmission range of each node were labeled as its neighbors (i.e. there is an edge between them). In order to ensure a various average connectivity, the transmission range was changed. Fig. 1. Example representative of strongly connected topologies Fig. 2. Example representative of averagely connected topologies 486 M. KENYERES, J. KENYERES, V. SKORPIL, DISTRIBUTED AGGREGATE FUNCTION ESTIMATION BY BIPHASICALLY … Fig. 3. Example representative of weakly connected topologies. In order to evaluate the achieved optimization, we compare our mechanism with other discussed weight models. We compare it with the Metropolis-Hasting weight model, the Maximum-degree weight model, and the Best Constant weight model. These models were chosen because all of them are classified as constant weight models [25], frequently-used and discussed in other papers and optimize average consensus with a similar principle as our mechanism. 5.1 Estimation of Average Value Within the first experiment, we examine the convergence rates and the number of the required messages. In this section, we draw our attention to the average estimation, compare our mechanism with the concurrent weight models and show the difference in the number of the messages that are necessary for the average consensus algorithm to achieve the consensus for each model. As the initial configuration of our optimized weight model depends on the placement of the identity numbers, we use randperm, which is a built-in function in Matlab to generate a vector with a random position of the numbers. Therefore, in order to ensure a higher credibility of our conclusions, our mechanism was executed 100 times (each execution is characterized by a shuffled set of the identity numbers) in each topology. In Tab. 1, our optimized weight model is compared with three other examined weight models. A positive value means an improvement ensured by our mechanism, meanwhile, a negative value indicates that a concurrent model achieves a faster convergence rate. We label these scenarios as positive, respectively, negative optimization. In Tab. 1, we have shown the average (calculated as the average value of the optimization of all ten networks within one set of the networks), the maximal (the optimization achieved in the network where the algorithm is optimized most significantly) and the minimal optimization (the optimization achieved in the network where the algorithm is optimized worse or even negatively) as well as the range of the optimization (all are expressed in %). Within this comparison, the average calculated from 100 convergence rates obtained within these executions is chosen as a representative of 100 executions of the Biphasically configured Metropolis-Hasting weight model. The complete results are shown in Appendix A. The column labeled as CR contains the convergence rates expressed as the number of the iterations. The column OPT [%] is formed by the relative optimizations [%] of our mechanism compared with the concurrent weight models. We can see from the results that our optimized weight model achieves a faster average convergence rate compared with all the concurrent models in all the sets of the networks. However, the Best Constant weight model achieves a faster rate in one of ten strongly connected topologies. The optimization is the most significant compared with the Maximum Degree weight model. Its average value ranges from 40.08% to 54.60%. The optimization of the Metropolis-Hasting model is from the range 8.19% to 15.45%. In this case, the most important fact is that our optimized weight model achieves a faster convergence rate in all the networks. The Best Constant weight model is optimized in the range 14.24% to 23.25%. Furthermore, we can see that (for all the weight models) the less connected the networks are, the higher average optimization our mechanism ensures. Let us focus on the range of the optimization (calculated as the difference between the maximal and the minimal optimization within one set). We can observe that this parameter achieves the highest value for the Best Constant weight model in all the sets. In Appendix A, we have also shown the number of the messages necessary for the average consensus algorithm to be completed. Since this parameter is closely related to the convergence rate, we do not provide a separate analysis. Consequently, we analyze the effect of a random shuffle of the identity numbers on the convergence rates. Convergence rate optimization of average estimation [%] in Weakly connected networks MD MH BC Average 54.60 % 15.45 % 23.25 % Maximum 64.60 % 21.28 % 35.95 % Minimum 47.38 % 12.65 % 8.77 % Range 17.22 % 8.63 % 27.18 % Convergence rate optimization of average estimation [%] in Averagely connected networks MD MH BC Average 46.84 % 11.15 % 15.36 % Maximum 60.04 % 15.64 % 36.21 % Minimum 39.71 % 7.80 % 1.42 % Range 20.33 % 7.84 % 34.79 % Convergence rate optimization of average estimation [%] in Strongly connected networks MD MH BC Average 40.08 % 8.19 % 14.24 % Maximum 69.30 % 16.93 % 49.29 % Minimum 31.47 % 4.29 % -3.42 % Range 37.83 % 12.64 % 52.62 % Tab. 1. Comparison of BMH with others – average est. RADIOENGINEERING, VOL. 26, NO. 2, JUNE 2017 487 Fig. 4. Percentage range of convergence rates caused by shuffle of identity numbers (average est.). The position of the identity numbers can affect only our optimized weight model – the other ones do not assume this value during the configuration process. In Fig. 4, we have depicted the range of the convergence rates calculated from the convergence rates obtained within 100 reparations (we depict the value of the ratio: the range/the average convergence rate expressed in % because the convergence rates differ from each other in different topologies). We can see that the range does not exceed 9 % of the average value in any case. This primarily affects the convergence rate in the weakly connected topologies. In the averagely and the strongly connected networks, its impact is negligible except for one topology of a strong connectivity. In the following part, we examine whether this shuffle can cause that there is an execution whose convergence rate is slower than the convergence rate of one of the concurrent models. Thus, we depict the slowest scenario within the Biphasically configured Metropolis-Hasting weight model with the fastest concurrent model in order to show that a shuffle does not cause a negative optimization when the average optimization is positive. The mentioned comparison for each set of the networks is depicted in Fig. 5 (weakly connected networks), Fig. 6 (averagely connected networks) and Fig. 7 (strongly connected networks). The white column represents the slowest convergence rate within the Biphasically configured Metropolis-Hasting weight model, while, the black one represents the fastest concurrent weight model. We can see that a positive optimization is preserved in all the cases when a positive average optimization is. Thus, a random allocation of the identity numbers has only a minimal impact on the convergence rate. Convergence rate Fig. 5. Comparison of slowest scenario of BMH with fastest concurrent model – weak connectivity – average est. Fig. 6. Comparison of slowest scenario of BMH with fastest concurrent model – average connectivity – average est. Networ k 1122334455667788991010 0 100 200 300 400 500 600 700 800 BMH - slowest scenario Fastest concurrent model Fig. 7. Comparison of slowest scenario of BMH with fastest concurrent model – strong connectivity – average est. 494 M. KENYERES, J. KENYERES, V. SKORPIL, DISTRIBUTED AGGREGATE FUNCTION ESTIMATION BY BIPHASICALLY … Appendix B The convergence rates of the network size estimation Maximum Degree Metropolis-Hasting Best Constant BMH CR OPT [%] CR OPT [%] CR OPT [%] CR Network #1 2930.50 56.19 1553.11 17.35 1733.61 25.95 1283.72 Network #2 1733.00 54.28 928.87 14.70 955.4 17.08 792.30 Network #3 2399.74 56.36 1250.58 16.26 1334.13 21.50 1047.27 Network #4 3435.97 56.15 1775.21 15.13 2159.40 30.23 1506.69 Weakly connected Network #5 3104.41 56.44 1605.87 15.79 1812.49 25.39 1352.31 Network #6 2278.28 47.75 1402.14 15.10 1383.30 13.94 1190.48 Network #7 4828.82 58.99 2502.52 20.86 3117.06 36.47 1980.42 Network #8 1206.68 49.27 740.11 17.29 726.44 15.73 612.15 Network #9 2788.35 64.46 1234.56 19.74 1516.88 34.67 990.91 Network #10 1498.38 50.58 860.51 13.95 861.25 14.03 740.43 Network #1 486.48 43.85 304.95 10.43 290.33 5.92 273.14 Network #2 563.86 44.62 350.46 10.90 314.70 0.78 312.26 Network #3 591.57 47.23 342.56 8.87 321.78 2.98 312.18 Network #4 815.70 51.73 439.52 10.41 475.28 17.15 393.76 Averagely connected Network #5 481.54 45.82 288.13 9.44 273.42 4.57 260.92 Network #6 1151.97 60.00 530.98 13.21 679.12 32.14 460.84 Network #7 562.69 42.26 365.90 11.20 341.12 4.75 324.91 Network #8 682.59 45.68 426.31 13.03 413.60 10.36 370.76 Network #9 444.69 36.90 311.02 9.78 286.20 1.96 280.60 Network #10 421.22 41.85 267.08 8.29 248.19 1.31 244.94 Network #1 339.23 53.51 175.28 10.02 231.73 31.94 157.72 Network #2 643.01 68.16 241.43 15.21 389.26 47.41 204.72 Network #3 208.51 36.87 138.93 5.25 136.07 3.26 131.64 Network #4 295.72 44.93 175.35 7.13 159.99 -1.78 162.84 Strongly connected Network #5 240.14 32.02 172.34 5.27 167.13 2.32 163.25 Network #6 284.81 39.61 183.94 6.49 174.17 1.25 172.00 Network #7 248.03 34.88 173.06 6.67 182.39 11.44 161.52 Network #8 279.66 40.23 182.12 8.21 177.24 5.69 167.16 Network #9 251.15 30.76 187.15 7.09 156.52 -11.1 173.89 Network #10 268.61 36.99 182.65 7.34 172.51 1.90 169.24 The number of the messages necessary for the network size estimation Maximum Degree Metropolis-Hasting Best Constant BMH MN MN MN MN Network #1 583169.50 309068.89 344988.39 255460.28 Network #2 344872.97 184845.13 190140.52 157667.70 Network #3 477548.26 248865.42 265491.87 208406.73 Network #4 683758.03 353266.79 429720.60 299831.31 Weakly connected Network #5 617777.59 319568.13 360685.51 269109.69 Network #6 453377.72 279025.86 275276.70 236905.52 Network #7 960935.18 498001.48 620294.94 394103.58 Network #8 240129.32 147281.89 144561.56 121817.85 Network #9 554881.65 245677.44 301859.12 197191.09 Network #10 298177.62 171241.49 171388.75 147345.57 Network #1 96809.52 60685.05 57775.67 54354.86 Network #2 112208.14 69741.54 62625.30 62139.74 Network #3 117722.43 68169.44 64034.22 62123.82 Network #4 162324.30 87464.48 94580.72 78358.24 Averagely connected Network #5 95826.46 57337.87 54410.58 51923.08 Network #6 229242.03 105665,02 135144.88 91707.16 Network #7 111975.31 72814.10 67882.88 64657.09 Network #8 135835.41 84835.69 82306.40 73781.24 Network #9 88493.31 61892.98 56953.80 55839.40 Network #10 83822.78 53148.92 49389.81 48743.06 Network #1 67506.77 34880.72 46114.27 31386.28 Network #2 127958.99 48044.57 77462.74 40739.28 Network #3 41493.49 27647.07 27077.93 26196.36 Network #4 58848.28 34894.65 31838.01 32405.16 Strongly connected Network #5 47787.86 34295.66 33258.87 32486.75 Network #6 56677.19 36604.06 34659.83 34228.00 Network #7 49357.97 34438.94 36295.61 32142.48 Network #8 55652.34 36241.88 35270.76 33264.84 Network #9 49978.85 37242.85 31147.48 34604.11 Network #10 53453.39 36347.35 34329.49 33678.76 RADIOENGINEERING, VOL. 26, NO. 2, JUNE 2017 495 Appendix C The range of the convergence rates caused by the change of the leader Maximum Degree Metropolis-Hasting Best Constant BMH CR OPT [%] CR OPT [%] CR OPT [%] CR Network #1 2700 50.77 1574 15.54 2518 47.21 1329.34 Network #2 1441 54.44 767 14.40 1343 51.11 656.53 Network #3 1560 60.75 714 14.24 1638 62.62 612.36 Network #4 3637 58.69 1822 17.54 3255 53.84 1502.48 Weakly connected Network #5 2814 49.48 1612 11.81 2625 45.84 1421.58 Network #6 1949 40.85 1267 9.01 1874 38.48 1152.87 Network #7 4996 65.26 2429 28.55 4795 63.80 1735.61 Network #8 512 49.63 271 4.84 630 59.07 257.88 Network #9 2158 62.78 1028 21.86 2036 60.55 803.26 Network #10 937 38.19 692 16.31 929 37.66 579.16 Network #1 345 43.09 205 4.22 267 26.46 196.34 Network #2 261 34.31 181 5.28 348 50.74 171.44 Network #3 335 48.96 185 7.58 386 55.70 170.98 Network #4 522 56.52 253 10.28 486 53.29 226.99 Averagely connected Network #5 379 57.87 176 9.28 273 41.51 159.67 Network #6 863 68.08 329 16.28 824 66.57 275.44 Network #7 248 50.77 138 11.54 259 52.86 122.08 Network #8 239 25.21 199 10.17 444 59.74 178.76 Network #9 207 31.83 147 4.00 258 45.30 141.12 Network #10 289 39.45 187 6.42 208 15.87 175.00 Network #1 560 83.19 136 30.79 283 66.74 94.13 Network #2 834 74.68 293 27.93 615 65.67 211.16 Network #3 119 52.16 60 5.12 92 38.12 56.93 Network #4 141 51.09 75 8.05 158 56.35 68.96 Strongly connected Network #5 863 68.08 329 16.28 824 66.57 275.44 Network #6 245 49.68 130 5.17 199 38.05 123.28 Network #7 196 41.24 126 8.60 143 19.47 115.16 Network #8 130 47.85 78 13.09 168 59.65 67.79 Network #9 111 22.72 88 2.52 144 40.43 85.78 Network #10 92 34.78 61 1.64 160 62.50 60.00