Stable Desynchronization for Wireless Sensor Networks: (III) Stability Analysis

TL;DR

This paper analyzes the stability of desynchronization algorithms in wireless sensor networks using dynamical systems, eigenvalue bounds, and stability theorems, providing conditions for system stability at equilibrium.

Contribution

It introduces a formal stability analysis framework for desynchronization algorithms, extending previous work with eigenvalue bounds and stability criteria.

Findings

Eigenvalues determine stability conditions.

Stability depends on the number of nodes within certain bounds.

The analysis confirms stability at equilibrium under specified conditions.

Abstract

In this paper, we use dynamical systems to analyze stability of desynchronization algorithms at equilibrium. We start by illustrating the equilibrium of a dynamic systems and formalizing force components and time phases. Then, we use Linear Approximation to obtain Jaconian (J) matrixes which are used to find the eigenvalues. Next, we employ the Hirst and Macey theorem and Gershgorins theorem to find the bounds of those eigenvalues. Finally, if the number of nodes (n) is within such bounds, the systems are stable at equilibrium. (This paper is the last part of the series Stable Desynchronization for Wireless Sensor Networks - (I) Concepts and Algorithms (II) Performance Evaluation (III) Stability Analysis)

Equations

$${\Delta }_1=$$

$${\Delta }_2=$$

$$⋮$$

$${\Delta }_{n-2}=$$

$${\Delta }_{n-1}=$$

$$\displaystyle+\frac{1}{\Delta_{n}}+\frac{1}{\Delta_{n}+\Delta_{n-1}}+\ldots+\frac{1}{\Delta_{n}+\Delta_{n-1}+\ldots+\Delta_{\frac{n}{2}+2}}\Bigg{)}$$

$${\Delta }_n=$$

$$\displaystyle+\frac{1}{\Delta_{n}}+\frac{1}{\Delta_{n}+\Delta_{n-1}}+\ldots+\frac{1}{\Delta_{n}+\Delta_{n-1}+\ldots+\Delta_{\frac{n}{2}+2}}\Bigg{)}$$

$${\Delta }_1=$$

$${\Delta }_2=$$

$$⋮$$

$${\Delta }_{n-2}=$$

$${\Delta }_{n-1}=$$

$$\displaystyle+\frac{1}{\Delta_{n}}+\frac{1}{\Delta_{n}+\Delta_{n-1}}+\ldots+\frac{1}{\Delta_{n}+\Delta_{n-1}+\ldots+\Delta_{\lceil\frac{n}{2}\rceil+1}}\Bigg{)}$$

$${\Delta }_n=$$

$$\displaystyle+\frac{1}{\Delta_{n}}+\frac{1}{\Delta_{n}+\Delta_{n-1}}+\ldots+\frac{1}{\Delta_{n}+\Delta_{n-1}+\ldots+\Delta_{\lceil\frac{n}{2}\rceil+1}}\Bigg{)}$$

$$J=\left(\begin{array}{ccccc}\frac{\partial {\Delta }_1}{\partial {\Delta }_1}& \frac{\partial {\Delta }_1}{\partial {\Delta }_2}& \cdots & \frac{\partial {\Delta }_1}{\partial {\Delta }_{n-1}}& \frac{\partial {\Delta }_1}{\partial {\Delta }_n}\\ \frac{\partial {\Delta }_2}{\partial {\Delta }_1}& \frac{\partial {\Delta }_2}{\partial {\Delta }_2}& \cdots & \frac{\partial {\Delta }_2}{\partial {\Delta }_{n-1}}& \frac{\partial {\Delta }_2}{\partial {\Delta }_n}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ \frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_1}& \frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_2}& \cdots & \frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_{n-1}}& \frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_n}\\ \frac{\partial {\Delta }_n}{\partial {\Delta }_1}& \frac{\partial {\Delta }_n}{\partial {\Delta }_2}& \cdots & \frac{\partial {\Delta }_n}{\partial {\Delta }_{n-1}}& \frac{\partial {\Delta }_n}{\partial {\Delta }_n}\end{array}\right)$$

$$\displaystyle{\frac{\partial\Delta_{i}}{\partial\Delta_{j}}}=\left\{\begin{array}[]{l l}1&\quad\text{if }j=i+1\\ 0&\quad\text{otherwise}\\ \end{array}\right.$$

$$\frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_1}=$$

$$=$$

$$\frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_2}=$$

$$=$$

$$⋮$$

$$\frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_{\frac{n}{2}-1}}=$$

$$=$$

$$\frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_{\frac{n}{2}}}=$$

$$\frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_{\frac{n}{2}+1}}=$$

$$\frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_{\frac{n}{2}+2}}=$$

$$=$$

$$⋮$$

$$\frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_{n-1}}=$$

$$=$$

$$\frac{\partial {\Delta }_{n-1}}{\partial {\Delta }_n}=$$

$$=$$

$$\frac{\partial {\Delta }_n}{\partial {\Delta }_1}=$$

$$\frac{\partial {\Delta }_n}{\partial {\Delta }_2}=$$

$$⋮$$

$$\frac{\partial {\Delta }_n}{\partial {\Delta }_{\frac{n}{2}-1}}=$$

$$\frac{\partial {\Delta }_n}{\partial {\Delta }_{\frac{n}{2}}}=$$

$$\frac{\partial {\Delta }_n}{\partial {\Delta }_{\frac{n}{2}+1}}=$$

$$\frac{\partial {\Delta }_n}{\partial {\Delta }_{\frac{n}{2}+2}}=$$

$$⋮$$

$$\frac{\partial {\Delta }_n}{\partial {\Delta }_{n-1}}=$$

$$\frac{\partial {\Delta }_n}{\partial {\Delta }_n}=$$

$$J=\left(\begin{array}{cccccccccc}0& 1& \cdots & 0& 0& 0& 0& \cdots & 0& 0\\ 0& 0& \ddots & 0& 0& 0& 0& \cdots & 0& 0\\ ⋮& ⋮& \ddots & \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 0& 1& 0& 0& \cdots & 0& 0\\ 0& 0& \cdots & 0& 0& 1& 0& \cdots & 0& 0\\ 0& 0& \cdots & 0& 0& 0& 1& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 0& 0& 0& 0& \cdots & 1& 0\\ A\sum _1& A\sum _2& \cdots & A\sum _{\frac{n}{2}-1}& 0& 0& -A\sum _{\frac{n}{2}-1}& \cdots & -A\sum _2& 1-A\sum _1\\ 1-A\sum _1& -A\sum _2& \cdots & -A\sum _{\frac{n}{2}-1}& 0& 0& A\sum _{\frac{n}{2}-1}& \cdots & A\sum _2& A\sum _1\end{array}\right).$$

$$J=\left(\begin{array}{ccccccccc}0& 1& \cdots & 0& 0& 0& \cdots & 0& 0\\ 0& 0& \ddots & 0& 0& 0& \cdots & 0& 0\\ ⋮& ⋮& \ddots & \ddots & ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 0& 1& 0& \cdots & 0& 0\\ 0& 0& \cdots & 0& 0& 1& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 0& 0& 0& \cdots & 1& 0\\ A\sum _1& A\sum _2& \cdots & A\sum _{\frac{n}{2}-1}& 0& -A\sum _{\frac{n}{2}-1}& \cdots & -A\sum _2& 1-A\sum _1\\ 1-A\sum _1& -A\sum _2& \cdots & -A\sum _{\frac{n}{2}-1}& 0& A\sum _{\frac{n}{2}-1}& \cdots & A\sum _2& A\sum _1\end{array}\right).$$

$$|J-\lambda I|$$

$$=\left|\begin{array}{cccccccccc}-\lambda & 1& \cdots & 0& 0& 0& 0& \cdots & 0& 0\\ 0& -\lambda & \ddots & 0& 0& 0& 0& \cdots & 0& 0\\ ⋮& ⋮& \ddots & \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & -\lambda & 1& 0& 0& \cdots & 0& 0\\ 0& 0& \cdots & 0& -\lambda & 1& 0& \cdots & 0& 0\\ 0& 0& \cdots & 0& 0& -\lambda & 1& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & \ddots & ⋮& ⋮\\ 0& 0& \cdots & 0& 0& 0& 0& \cdots & 1& 0\\ 1& 0& \cdots & 0& 0& 0& 0& \cdots & -\lambda & 1-\lambda \\ 1-A\sum _1& -A\sum _2& \cdots & -A\sum _{\frac{n}{2}-1}& 0& 0& A\sum _{\frac{n}{2}-1}& \cdots & A\sum _2& A\sum _1-\lambda \end{array}\right|$$

$$=\left|\begin{array}{cccccccccc}-\lambda +\frac{1}{{\lambda }^{n-2}}& 0& \cdots & 0& 0& 0& 0& \cdots & 0& \frac{1-\lambda }{{\lambda }^{n-2}}\\ \frac{1}{{\lambda }^{n-3}}& -\lambda & \cdots & 0& 0& 0& 0& \cdots & 0& \frac{1-\lambda }{{\lambda }^{n-3}}\\ \frac{1}{{\lambda }^{n-4}}& 0& \ddots & 0& 0& 0& 0& \cdots & 0& \frac{1-\lambda }{{\lambda }^{n-4}}\\ ⋮& ⋮& \ddots & \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ \frac{1}{{\lambda }^{\frac{n}{2}-2}}& 0& \cdots & 0& -\lambda & 0& 0& \cdots & 0& \frac{1-\lambda }{{\lambda }^{\frac{n}{2}-2}}\\ \frac{1}{{\lambda }^{\frac{n}{2}-3}}& 0& \cdots & 0& 0& -\lambda & 0& \cdots & 0& \frac{1-\lambda }{{\lambda }^{\frac{n}{2}-3}}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & \ddots & ⋮& ⋮\\ \frac{1}{\lambda }& 0& \cdots & 0& 0& 0& 0& \ddots & 0& \frac{1-\lambda }{\lambda }\\ 1& 0& \cdots & 0& 0& 0& 0& \cdots & -\lambda & 1-\lambda \\ 1-A\sum _1& -A\sum _2& \cdots & -A\sum _{\frac{n}{2}-1}& 0& 0& A\sum _{\frac{n}{2}-1}& \cdots & A\sum _2& A\sum _1-\lambda \end{array}\right|$$

$$=\left|\begin{array}{cccccccccc}-\lambda +\frac{1}{{\lambda }^{n-2}}& 0& \cdots & 0& 0& 0& 0& \cdots & 0& \frac{1-\lambda }{{\lambda }^{n-2}}\\ {\lambda }^2& -\lambda & \cdots & 0& 0& 0& 0& \cdots & 0& 0\\ {\lambda }^3& 0& \ddots & 0& 0& 0& 0& \cdots & 0& 0\\ ⋮& ⋮& \ddots & \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ {\lambda }^{\frac{n}{2}}& 0& \cdots & 0& -\lambda & 0& 0& \cdots & 0& 0\\ {\lambda }^{\frac{n}{2}+1}& 0& \cdots & 0& 0& -\lambda & 0& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & \ddots & ⋮& ⋮\\ {\lambda }^{n-2}& 0& \cdots & 0& 0& 0& 0& \ddots & 0& 0\\ {\lambda }^{n-1}& 0& \cdots & 0& 0& 0& 0& \cdots & -\lambda & 0\\ 1-A\sum _1& -A\sum _2& \cdots & -A\sum _{\frac{n}{2}-1}& 0& 0& A\sum _{\frac{n}{2}-1}& \cdots & A\sum _2& A\sum _1-\lambda \end{array}\right|$$

$$=\left|\begin{array}{cccccccccc}-\lambda +\frac{1}{{\lambda }^{n-2}}& 0& \cdots & 0& 0& 0& 0& \cdots & 0& \frac{1-\lambda }{{\lambda }^{n-2}}\\ {\lambda }^2& -\lambda & \cdots & 0& 0& 0& 0& \cdots & 0& 0\\ {\lambda }^3& 0& \ddots & 0& 0& 0& 0& \cdots & 0& 0\\ ⋮& ⋮& \ddots & \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ {\lambda }^{\frac{n}{2}}& 0& \cdots & 0& -\lambda & 0& 0& \cdots & 0& 0\\ {\lambda }^{\frac{n}{2}+1}& 0& \cdots & 0& 0& -\lambda & 0& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & \ddots & ⋮& ⋮\\ {\lambda }^{n-2}& 0& \cdots & 0& 0& 0& 0& \ddots & 0& 0\\ {\lambda }^{n-1}& 0& \cdots & 0& 0& 0& 0& \cdots & -\lambda & 0\\ B& 0& \cdots & 0& 0& 0& 0& \cdots & 0& A\sum _1-\lambda \end{array}\right|$$

Full text

I Stability Analysis of DWARF (the Single-Hop Desynchronization Algorithm)

To prove that the system is stable, we begin by transforming the system into a non-linear dynamic system. The evenness of the number of nodes affects the analysis. Therefore, we divide the non-linear dynamic system into two cases: when $n$ is even and when $n$ is odd where $n$ is the number of nodes.

1. when $n$ is even: Figure 1a illustrates the equilibrium of a dynamic system when the number of nodes is even. Noticeably, node 0 and $n/2$ are exactly at the opposite side of each other. At the first snapshot of the system, node 0 adjusts its phase based on the force function of DWARF. After adjustment, we re-label node 1 to 0, node 2 to 1, …, node $n-1$ to $n-2$, and node 0 to $n-1$ for analysis at the next snapshot (see Figure 1b). Therefore, to transform into a difference equation, $\Delta_{1}$ in the next snapshot is $\Delta_{2}$ in the previous snapshot, $\Delta_{2}$ in the next snapshot is $\Delta_{3}$ in the previous snapshot, and so on. However, $\Delta_{n-1}$ and $\Delta_{n}$ in the next snapshot are $\Delta_{n}$ and $\Delta_{1}$ in the previous snapshot adjusted by the force function, respectively.

Therefore, the non-linear dynamic system when $n$ is even can be expressed as the following difference equations:

[TABLE]

We note that the force from node $n/2$ in the previous snapshot is already balanced when we consider node 0. Therefore, $\Delta_{\frac{n}{2}}$ and $\Delta_{\frac{n}{2}+1}$ do not appear in the force equation for adjusting $\Delta_{n-1}$ and $\Delta_{n}$ in the next snapshot.

2. when $n$ is odd: Similarly, Figure 2 illustrates the non-linear dynamic system when $n$ is odd. The only difference from when $n$ is even is that there is no node that is opposite to node 0. Therefore, in the difference equations, only $\Delta_{\lceil\frac{n}{2}\rceil}$ from the previous snapshot does not appear in the force equation for adjusting $\Delta_{n-1}$ and $\Delta_{n}$ in the next snapshot.

[TABLE]

Due to the non-linear adaptation function, the standard linear dynamic system analysis does not suffice. Therefore, we locally analyse stability of the system around a fixed point which is the equilibrium point. We begin by linear approximation to find the Jacobian at the equilibrium. Then, from the Jacobian, we find the bound of eigenvalues which is the crucial part of stability analysis.

I-A Linear Approximation

The Jacobian ($J$) of a difference equations system is defined as follows:

[TABLE]

We consider the Jacobian when $n$ is even and odd separately. In both cases, after finding the Jacobian, we substitute each $\Delta_{i}$ with $T/n$ which is the phase interval between each node at the equilibrium.

I-A1 $n$ is even

For $\Delta_{i},1\leq i\leq n-2$, partial derivatives of ${\partial\Delta_{i}/\partial\Delta_{j}},1\leq j\leq n$ are the followings:

[TABLE]

For $\Delta_{n-1}$, we find its partial derivatives as follows:

[TABLE]

For $\Delta_{n}$, we find its partial derivatives which is similar to $\Delta_{n-1}$ as follows.

[TABLE]

Let $\sum\limits_{s}$ stands for $\sum\limits_{i=s}^{\frac{n}{2}-1}{\frac{1}{i^{2}}}$, where $s\in[1,\frac{n}{2}-1]$ and let $A=Kn^{2}/T$. The Jacobian matrix of the system, when $n$ is even, is

[TABLE]

I-A2 $n$ is odd

Similarly, when $n$ is odd, we can find its Jacobian with the same procedure as when $n$ is even.We obtain the similar Jacobian matrix as follows:

[TABLE]

I-B Finding Eigenvalues

After finding the Jacobian of the linear approximation at the equilibrium, we find the eigenvalues by solving the equation $|J-\lambda I|=0$. If all eigenvalues lay on a unit circle, the system is stable at the equilibrium. Therefore, we begin by finding the determinant of the matrix $J-\lambda I$. We use row operations to transform the determinant of $|J-\lambda I|$ into the determinant of a triangular matrix. Then, the determinant of the triangular matrix is the multiplication of diagonal entries.

The following procedure is to transform the determinant of $|J-\lambda I|$ into the determinant of triangular matrix when $n$ is even:

[TABLE]

where

[TABLE]

Similarly, we use the same procedure to find the determinant of $|J-\lambda I|$ when $n$ is odd. The result matrix is the same as above except that

[TABLE]

We note that the difference between Equation 14 and 15 is the exponents of $\lambda$ at two middle terms.

The determinant of a triangular matrix is the product of diagonal terms. That is

[TABLE]

For terms $-\lambda$, it is obvious that we have repeated eigenvalues $\lambda$ that are equal to [math] and lay on a unit circle.

For the last term $A\sum\limits_{1}-\lambda$, the value of $\lambda$ depends on $n$ as follows:

[TABLE]

In this case, $\lambda\geq 0$, however, we have to find the least upper limit of $n$ that leads to $\lambda\leq 1$.

[TABLE]

Therefore, if the number of nodes $n$ is less than $3.178\times 10^{9}$ nodes, the eigenvalue $\lambda$ lies between 0 and 1 in this case.

For the first term, $-\lambda+\frac{1}{\lambda^{n-2}}+\frac{B}{A\sum\limits_{1}-\lambda}\left(\frac{\lambda-1}{\lambda^{n-2}}\right)$, the $\lambda$ must not be [math]. If $n$ is even, we derive the polynomial of $\lambda$ as follows:

[TABLE]

To get the eigenvalues, we have to solve the polynomial

[TABLE]

Similarly, if $n$ is odd, we have to solve the polynomial

[TABLE]

However, according to the Abel’s Impossibility theorem (or Abel-Ruffini theorem), we cannot find a general algebraic solution to polynomial equations of degree five or higher ([3]). Therefore, instead of finding the exact values, we find the upper and lower bound of the eigenvalues to ensure that they lay on a unit circle.

I-C The Bound of Eigenvalues

To find the bound of polynomial roots (i.e. our eigenvalues), we follow the following theorem of Hirst and Macey ([1]).

Theorem I.1** (Hirst and Macey Bound).**

Given $f:\mathbb{C}\rightarrow\mathbb{C}$ defined by $f(z)=z^{n}+a_{n-1}z^{n-1}+\ldots+a_{1}z+a_{0}$, where $a_{0},a_{1},\ldots,a_{n}\in\mathcal{C}$, and $n$ a positive integer. If $z$ is a zero of $f$, then

[TABLE]

Therefore, to bound all eigenvalues within a unit circle, $\sum_{i=0}^{n-1}|a_{i}|$ must be less than or equal 1. From both Equation 17 and 18,

[TABLE]

where $A=Kn^{2}/T=38.597n^{-1.874}(T/1000)(n^{2}/T)=0.038597n^{0.126}$.

For the first case, $-1+2A\sum_{i=1}^{\frac{n}{2}-1}\frac{1}{i^{2}}\leq 2A\sum_{i=1}^{\frac{n}{2}-1}\frac{1}{i^{2}}-1$ is always true regardless of the number of nodes $n$. Therefore, we consider the second case:

[TABLE]

When $n$ is large, $\sum_{i=1}^{\frac{n}{2}-1}\frac{1}{i^{2}}$ is approximately equal to the Reimann zeta function $\zeta(2)=\sum_{i=1}^{\infty}\frac{1}{i^{2}}=\pi^{2}/6\approx 1.645$. Therefore, we get the following:

[TABLE]

II Stability Analysis of M-DWARF (the Multi-Hop Desynchronization Algorithm)

The multi-hop stability analysis in this chapter is more complicated than the analysis of single-hop networks due to the connectivity and topology which affect the analysis. We begin by formalizing components of forces to a general form. Then, we transform the system into a dynamic system equations and locally analyse at the equilibrium.

II-A Definition

We first define several variables and notations to be used throughout the analysis in this section.

We define $n$ to be the number of nodes in the system and $T$ be the time period. We assume $n$ is even. However, in the case of odd $n$, the system can be analysed with the same procedure.

We define $c_{i,j}$ to be a connectivity status in two-hop communication from node $j$ to node $i$. If node $i$ perceives the presence of node $j$ (as a one-hop neighbor or by relaying relative phase), $c_{i,j}$ is 1. Otherwise, $c_{i,j}$ is 0. For example, in a simple 3-node chain network with node labelled 0, 1, and 2 respectively, if communication links are symmetry, the values of $c_{0,1},c_{1,0},c_{1,2}$, and $c_{2,1}$ are 1 whereas the values of $c_{0,2}$ and $c_{2,0}$ are 0. We assume that nodes are labelled sequentially in increasing order from [math] to $n-1$ according the the increasing phase in the global period ring.

We define the modulo notation $[x]_{n}$ to stand for $x\mod n$ for brevity.

We define $\Delta_{i}$ to be a phase interval between node $i$ and node $[i+1]_{n}$, where $\Delta_{i}\in(0,T]$.

II-B Force Component

To construct a dynamic system, we first analyse how force with the absorption mechanism affects the system equations. As described in Section 4.2.3, a force can be absorbed if it is not originated from the closest phase neighbor. Therefore, respect to the node $i$’s point of view, some forces have no effect to node $i$ whereas some forces do. We classify forces that affect node $i$ into three components: closest component, resistance component, and absorption component. The final form of forces that affect node $i$ will be in the following form:

[TABLE]

Closest component: Respect to the node $i$’s point of view, the closest component is the force from the closest phase neighbor of node $i$ and node $i$ perceives its presence. This force is not absorbed by any node. This force component from node $j$ to node $i$ will appear in the equation if, between node $j$ and node $i$, node $i$ does not perceive the presence of any node. In other words, node $j$ is the closest perceived phase neighbor of node $i$. We firstly define the closest component for node $i$ by using the combination of logic and algebraic expression as follows (we use it only for clarification purpose and we will change it to pure algebraic expression later),

[TABLE]

where $f_{i,[j]_{n}}^{+}$ is a positive (clockwise) force and $f_{i,[j]_{n}}^{-}$ is a negative (counter-clockwise) force. We note that $f_{i,[j]_{n}}^{+}$ appears only when $c_{i,[j]_{n}}=1$ and all $c_{i,[k]_{n}}=0$, where $j<k<i+n$. Similary, $f_{i,[j]_{n}}^{-}$ appears only when $c_{i,[j]_{n}}=1$ and all $c_{i,[k]_{n}}=0$, where $i<k<j$.

Then, we convert the logic expression into the algebraic expression. For logical negation, $\neg c_{i,[j]_{n}}$ can be algebraically expressed as $(1-c_{i,[j]_{n}})$. For logical and ($\land$), we can express algebraically by using multiplication instead. For example, $c_{i,[j]_{n}}\land c_{i,[k]_{n}}$ can be expressed as $c_{i,[j]_{n}}c_{i,[k]_{n}}$.

Therefore, Equation LABEL:eq:closest can be expressed algebraically as the following,

[TABLE]

Let $R_{i,j}^{+}$ be $\prod_{k=1}^{i+n-(j+1)}(1-c_{i,[i+(n-k)]_{n}})\in{0,1}$ and $R_{i,j}^{-}$ be $\prod_{k=1}^{j-1}(1-c_{i,[i+k]_{n}})\in{0,1}$. We derive the following form of the closet component to be used in our analysis,

[TABLE]

Resistance component: Respect to the node $i$’s point of view, there is a resistance component originated from node $j$ if the following criteria are satisfied:

• •

Node $i$ perceives the presence of node $j$.

• •

There is at least one node $k$ following node $j$ in the time period ring (in each force direction).

• •

Node $i$ perceives the presence of node $k$.

For example, with respect to the node [math]’s point of view, if node $1$ is the closest phase neighbor of node [math] and there is node $2$ following node $1$ in clockwise direction, and node [math] perceives the presence of both node $1$ and node $2$, then, there is a force difference between node $1$ and node $2$ in the form of $f_{0,1}-f_{0,2}$ (see Section 4.2.3). The part $f_{0,1}$ is called resistance component and the part $f_{0,2}$ is called absorption component which we will describe later. We note that, if there is no node following the closest phase neighbor in each direction, there is no resistance and absorption components.

Therefore, we define the resistance component for node $i$ by using the combination of logic and algebraic expression as follows,

[TABLE]

Then, we convert the logical expression into the algebraic expression.

Let $s^{+}(v_{i,k}^{+})$ be an algebraic or function to represent a logical or expression of $(c_{i,[i+(n-k)]_{n}}\lor c_{i,[i+(n-k-1)]_{n}}\lor\cdots\lor c_{i,[i+(\frac{n}{2}+1)]_{n}})$, where $v_{i,k}^{+}=c_{i,[i+(\frac{n}{2}+1)]_{n}}2^{\frac{n}{2}-1-k}+\cdots+c_{i,[i+(n-k-1)]_{n}}2^{1}+c_{i,[i+(n-k)]_{n}}2^{0}$.

Similarly, let $s^{-}(v_{i,k}^{-})$ be an algebraic function to represent a logical or expression of $(c_{i,[i+k]_{n}}\lor c_{i,[i+(k+1)]_{n}}\lor\cdots\lor c_{i,[i+(\frac{n}{2}-1)]_{n}})$, where $v_{i,k}^{-}=c_{i,[i+(\frac{n}{2}-1)]_{n}}2^{\frac{n}{2}-1-k}+\cdots+c_{i,[i+(k+1)]_{n}}2^{1}+c_{i,[i+k]_{n}}2^{0}$.

We note that, if we write a binary string $c_{i,[i+(\frac{n}{2}+1)]_{n}}\cdots c_{i,[i+(n-k-1)]_{n}}c_{i,[i+(n-k)]_{n}}$, $v_{k}^{+}$ is an integer value in base 10 of this binary string where $v_{i,k}^{+}\in\{0,1,\cdots,2^{\frac{n}{2}-k}-1\}$. If all $c_{i,[j]_{n}}=0$, then $v_{i,k}^{+}=0$. The result is similar for $v_{i,k}^{-}$.

Let $I_{A}(x)$ be an indicator function as follows,

[TABLE]

Therefore, $s^{+}(v_{i,k}^{+})=I_{\mathbb{N}-\{0\}}(v_{i,k}^{+})$ and $s^{-}(v_{i,k}^{-})=I_{\mathbb{N}-\{0\}}(v_{i,k}^{-})$.

From Equation 26, we derive the equivalent algebraic expression as follows,

[TABLE]

Let $S_{i,j}^{+}=s^{+}(v_{i,i+n-j+1}^{+})\in\{0,1\}$ and $S_{i,j}^{-}=s^{-}(v_{i,j-i+1}^{-})\in\{0,1\}$. We derive the following form for the resistance component to be used in our analysis,

[TABLE]

Absorption component: Respect to the node $i$’s point of view, there is an absorption component originated from node $j$ if the following criteria are satisfied:

• •

Node $i$ perceives the presence of node $j$.

• •

There is at least one node $k$ stays between node $i$ and node $j$ in the time period ring (in each force direction).

• •

Node $i$ perceives the presence of node $k$.

Therefore, we define the absorption component for node $i$ by using the combination of logic and algebraic expression as follows,

[TABLE]

Based on the same procedure deriving the resistance component, we derive the following form for the absorption component to be used in our analysis,

[TABLE]

where $T_{i,j}^{+}$ and $T_{i,j}^{-}$ are the algebraic expressions of the logical or terms in Equation 30.

Therefore, respect to the node $i$’s point of view, the total force at node $i$ is as follows,

[TABLE]

In the next section, we analyse the stability of the M-DWARF algorithm based on the derived total force.

II-C Stability Analysis

As same as the analysis of single-hop networks, to prove that the system is stable, we begin by transforming the system into a non-linear dynamic system.

Let $f_{i,j}^{+}$ and $f_{i,j}^{-}$ be the positive and negative forces from node $j$ to node $i$ respectively,

[TABLE]

where $\Delta_{k}\in(0,T]$ is the phase difference between node $[k+1]_{n}$ and node $k$.

From Equation 32, the total force at node $i$ is

[TABLE]

Let $\Delta_{i}$ be the current phase difference between node $[i+1]_{n}$ and node $i$ and $\Delta_{i}^{\prime}$ be the phase difference between node $[i+1]_{n}$ and node $i$ in the next time period. The dynamic system of a multi-hop network running the M-DWARF algorithm is as follows:

[TABLE]

We write the transition of $\Delta_{i}$ in a general form as the following,

[TABLE]

Then, we linearly approximate the dynamic system at the equilibrium.

II-C1 Linear Approximation

The Jacobian ($J$) of a difference equations system is defined as follows:

[TABLE]

Therefore, from Equation 36, each element in a row of Jacobian is

[TABLE]

For the first term, $\frac{\partial\Delta_{i}}{\partial\Delta_{p}}$ is 1 if $p=i$. . If $p\neq i$, this term is zero. Formally,

[TABLE]

Then, we find $\frac{\partial F_{i}}{\partial\Delta_{p}}$ as follows:

[TABLE]

For $\frac{\partial F_{[i+1]_{n}}}{\partial\Delta_{p}}$, we can find by replacing $i$ with $i+1$ in Equation 40. The result is as follows:

[TABLE]

Consequently, from Equation 38, 39, 40, and 41, we find $\frac{\partial\Delta_{i}^{\prime}}{\partial\Delta_{p}}$ as the following:

[TABLE]

At this point, we already get the general form of each row of the Jacobian matrix. To locally approximate at the equilibrium, we have to substitute $c_{i,j},R_{i,j},S_{i,j},T_{i,j}$, and $\Delta_{i}$ for all $i,j\in{0,1,\cdots,n-1}$ where $j\neq i$ with the value at the equilibrium state. Then, we have find the maximum value of $n$ that bounds the eigenvalues within a unit circle to guarantee the stability. However, such values depend on the topology. For example, the connectivity $c_{i,j}$ depends on whether node $i$ can perceive the presence of node $j$ or not. If node $i$ and $j$ are within two-hop communication, this value is 1, otherwise, this value is 0.

In this paper, we prove the stability of the multi-hop star topology as an example because this topology is less complicated to prove and may be the easiest example for the reader to follow the proof. For other topologies, we conjecture that the local stability can also be assured as well.

For the multi-hop star topology, the time interval between arbitrary two consecutive nodes is $T/n$ at the equilibrium. Therefore, each row of the Jacobian matrix from Equation 42 is changed by substituting all $\Delta_{k}$ with $T/n$ and re-indexing the summation as follows,

[TABLE]

Then, for the star topology, nodes perceive the presence of all two-hop neighbors. Therefore, all terms $c_{i,j},R_{i,j},S_{i,j}$, and $T_{i,j}$ are 1. Consequently, the row of the Jacobian matrix at the equilibrium is the following:

[TABLE]

The result of the Jacobian matrix at the equilibrium is the circulant matrix as shown below:

[TABLE]

where the diagonal entry D_{0}=1-2A\Bigg{(}1+\sum_{j=1}^{n-2}1/j^{2}\Bigg{)} and $A=Kn^{2}/T$.

We can find each eigenvalue of a circulant matrix with a general solution presented in [4]. However, we only need to guarantee that all eigenvalues lay on a unit circle. Therefore, we find the bound of eigenvalues instead.

II-C2 The Bound of Eigenvalues

To find the bound of an $n\times n$ matrix, we use the Gershgorin’s Circle Theorem ([2, 5])

Theorem II.1** (Gershgorin’s Theorem Round 1).**

Every eigenvalue $\lambda$ of $n\times n$ matrix $A$ satisfies:

[TABLE]

In other words, every eigenvalue lies within at least one of Gershgorin discs centered at $A_{i,i}$ with radius $\sum_{j\neq i}|A_{i,j}|$, where $A_{i,i}$ is the diagonal entry of a matrix.

In our circulant matrix, all diagonal entries and the sums of elements in each row and each column are the same. Therefore, in our matrix, all Gershgorin’s discs are centered at $D_{0}=1-2A(1+\sum_{j=1}^{n-2}1/j^{2})$ with radius $r=2A(1+\sum_{j=1}^{\frac{n}{2}-2}1/j^{2})$, where $A=Kn^{2}/T$.

Then, we find the maximum number of nodes $n$ that guarantees the Gershgorin’ discs are in a unit circle.

Let $\vec{D}$ be a vector drawn from the origin $(0,0)$ to the center of the Gershgorin’s disc $(1-2A(1+\sum_{j=1}^{n-2}1/j^{2}),0)$. Due to the imaginary part of $\vec{D}$ is zero, the magnitude $|\vec{D}|$ is $|1-2A(1+\sum_{j=1}^{n-2}1/j^{2}|$. Therefore, we derive the following to guarantee the Gershgorin’s discs are in a unit circle:

[TABLE]

From the M-DWARF algorithm, we substitute $K$ with $38.597\times n^{-1.874}\times T/1000$. Additionally, when $n$ is large, the value of $\sum_{j=1}^{n-2}1/j^{2}$ and $\sum_{j=1}^{\frac{n}{2}-2}1/j^{2}$ converge to the Reimann zeta function $\zeta(2)=\sum_{i=1}^{\infty}1/i^{2}=\pi^{2}/6\approx 1.645$. From Equation 45, we get the following:

[TABLE]

The condition $1-0.20417813n^{0.126}\leq\text{ }1-0.20417813n^{0.126}$ is always true regardless of the number of nodes $n$. Then, we consider the following condition:

[TABLE]

Therefore, if the number of nodes is less than $2.99\times 10^{5}$ nodes, every eigenvalue is guaranteed to lay in a unit circle. In other words, the non-linear dynamic system for the multi-hop star topology is locally stable at the equilibrium. If there is a small perturbation around the equilibrium, the system is able to converge back to the equilibrium.

To prove the stability of other topologies, we can substitute $\Delta_{i}$ and $c_{i,j}$ in the Jacobian matrix with the value at the equilibrium of each topology. Then, finding the eigenvalues of the substituted Jacobian matrix. If we can bound that every eigenvalue lies in a unit circle, the algorithm is locally stable for such topologies. We conjecture that the proof of other topologies is similar to the proof of the star topology with the similar procedure.