Modeling of Multihop Wireless Sensor Networks with MAC,Queuing,and Cooperation

Abstract

We present a Markovian decision process (MDP) framework for multihop wireless sensor networks (MHWSNs) to bound the network performance of both energy constrained (EC) networks and energy harvesting (EH) networks, both with and without relay cooperation. The model provides the fundamental performance limit that a MHWSN can theoretically achieve, under the general constraints from medium access control, routing, and energy harvesting. We observe that the analyses for EC and EH networks fall into two branches of MDP theory, which are finite-horizon process and infinite-horizon process, respectively. The performance metrics for EC and EH networks are different. For EC networks, an appropriate metric is the network lifetime; for EH networks, an appropriate metric is, for example, the network throughput. To efficiently solve the models with high dimension, for the EC networks, we propose a novel computational algorithm by taking advantage of the stochastic shortest path structure of the problem; for the EH networks, we propose a dual linear programming based algorithm by considering the sparsity of the transition matrix. Under the unified MDP framework, numerical results demonstrate the advantages of cooperation for the optimal performance, in both EC and EH networks.

1. Introduction

This paper considers modeling the optimal performance of multihop wireless sensor networks (MHWSNs), which can be energy constrained (EC) networks or energy harvesting (EH) networks, both with and without relay cooperation. The network performance of a MHWSN is a complex function of sensors' harvested energy, traffic volume, routing protocol, and medium access (MAC) technique. An analytical approach that is suitable for WSN traffic and that derives the optimum would serve as a valuable benchmark for heuristics MAC and routing protocols. However, such a multihop framework has not been presented. In this paper, we explore the optimum by presenting a unified Markov Decision Process (MDP) analysis that can analyze both EC and EH networks. We observe that the treatments for both networks fall into two branches of MDP theory, the finite-horizon process and the infinite-horizon process, respectively.

Existing analyses have mainly focused on a single node [1], a single link [2], or a single-hop network [3]. Departing from the literature, this paper addresses the general multihop network with neither assumptions on node distribution nor assumptions on total traffic pattern, motivated by the following reasons. First, since the communication range of a single node is limited, single-hop deployment [3] cannot always meet the requirements of many applications, such as environmental and structural monitoring, security and defense, and wireless body area networks [4]. Second, the traffic pattern of nodes in a MHWSN is unbalanced, invalidating the homogeneous or deterministic traffic assumptions, which are widely assumed for wireless ad hoc or mesh networks. In particular, some nodes near the Sink consume their energy faster because they must use their energy to relay other nodes' data. In EC networks, this unbalance creates the well-known “energy hole” problem [5]. In EH networks, the impact on energy consumption is still nontrivial. Moreover, the extent to which the traffic is unbalanced is not known a priori and is usually determined by a routing protocol. Therefore, towards quantifying the network performance, careful consideration should be taken to capture the peculiarity of the traffic pattern. Third, because the network's energy consumption rate is closely related to nodes' interaction in various aspects including channel access, routing, and cooperation, an optimal solution obtained by focusing only on a point-to-point link does not imply the global optimum for the network-wide objective. In fact, quantifying the latter problem is fairly challenging.

On another track, cooperative transmission (CT) is a mixture of communication protocol and physical (PHY) layer combining scheme that allows spatially separated nodes to collaborate to transmit a single message by forming a virtual multiple-input-single-output (VMISO) array. A receiver, gaining a signal-to-noise ratio (SNR) advantage of 10 dB to 20 dB through PHY combining [6], can decode the message at an extended range [7]. CT range extension has been experimentally demonstrated in [8], and it shows significant impact on Layer Two and Layer Three of a WSN by eliminating the “energy hole” in EC network [9, 10] and by providing better Quality of Service in EH networks [11]. Time-division CT (TDCT) is one type of cooperation, in which the source node multicasts the packet, and the cooperators that decode the packet forward it in orthogonal time slot [12]. Our previous work OSC-MAC [13] shows notable lifetime improvement over state-of-the-art MAC protocols by using TDCT to balance energy. Yet, it is noted that the performance gap between the CT MAC protocol [9, 10] with fixed routes and the theoretically optimal value is unknown. Also note that the analysis in [11], among others, assumes perfect link scheduling, which contrasts with the practical situation wherein link activities are subject to MAC constraints and packet collisions occur. All the above discussions have further motivated our analysis to comprehensively consider energy harvesting, routing, MAC, and cooperation, from the network perspective. Our previous work [14] models the lifetime of battery based multihop sensor networks, which, however, does not apply to energy unconstrained networks (e.g., through energy harvesting). This paper extends the model in [14] to provide an analysis framework that apply to both energy constrained (EC) networks and energy harvesting (EH) networks.

The rest of the paper is organized as follows. Section 2 presents the related work. Section 3 describes the system model and assumptions. An MDP formulation is detailed in Section 4. The computational method and numerical evaluations for non-CT and CT networks are presented in Sections 5 and 6, respectively. Finally, Section 7 concludes the paper.

2. Related Work

In this section, we describe existing works related to modeling for wireless sensor networks in terms of the scope of modeling, the assumptions, and the technical tools. Then we highlight the challenges solved and contributions made in this work in comparison.

Existing analyses are mainly limited to a single node, a single link, or a single-hop network [15]. For example, [1] focuses on energy management of a single energy harvesting (EH) node for achieving maximum throughput and minimizing delay, using an MDP [16] model. While it is assumed in [1] that the characteristic of traffic and energy is stationary (same as in this work), it also assumes infinite packet queue and energy queue, which is not practical. In [17], a method based on stochastic timed automata and statistical model checking is proposed to model a WSN protocol called timing-sync protocol, which however focuses only on node behavior. In [18], an MDP model is used to model the sensor processor (CPU) energy and Petri nets are used to model an energy constrained (EC) sensor node. In [2], an MDP model is presented for EH tags in a point-to-point link, assuming a time-slotted system and a stationary energy harvesting model. Transmission of EH nodes in a single link and in the context of a fading channel is tackled in [19]. In [3], an MDP model is used to formulate the lifetime of a single-hop EC network, which considers sensor scheduling in a scenario where only a fraction of sensors collect information and communicate directly with a one-hop Sink. In [20], considering nodes that are equipped with a hybrid energy storage system, the authors provide an MDP model for a single-hop EH network. The major limitations of these works are their inability to analyze the performance of a network that extends beyond one-hop deployment and the lack of analysis for cooperative transmission.

For multihop networks, the past analyses have exploited various approaches. The optimal lifetime of a multihop WSN is derived in [21] using linear programming (LP) formulation based on traffic balance conditions. However, the analyses in [21] are limited to only the routing layer. In [22], a model based on mixed integer linear programming is proposed to identify energy-efficient route from the source sensor to the Sink; yet, it considers only energy minimization and does not produce network-wide objective optimization. In [23], the authors provide a Lyapunov analysis for multihop EH WSNs, which considers the dynamics of energy and queue, to maximize the utility. However, [23] is based on queue stability, which is not suitable for low-traffic WSNs. Reference [24] studies the upper bound on the lifetime of data gathering sensor networks, while [25] derives the upper bound on the lifetime of large-scale networks using the theory of coverage process. Both [24, 25] are based on the assumption that the data sources are deployed with a particular probability density function or process; thus they are unable to capture the inhomogeneous traffic characteristic in a multihop network. We note that all the aforementioned works fail to consider MAC constraints, as they assume perfect link scheduling and thus ignore packet losses and retransmissions, which contrasts the practical situation wherein link activities are subject to MAC constraints and packet collisions occur. Moreover, none of these works, except for [21], considers cooperative transmission.

For cooperative networks, the existing analyses are typically limited to single-hop networks. In [26], the authors consider the lifetime maximization in an amplify-and-forward cooperative network and model the energy dissipation of nodes as a Markov chain. In [27], the authors develop a general probability model to study the outage performance of the best-relay selection with adaptive decode-and-forward cooperative network. Relatively less works have focused on analysis of multihop cooperative networks. The LP model in [28] has been recently extended to multihop CT networks [29] by considering all single-input-single-output (SISO) and virtual multiple-input-single-output (VMISO) links. Again, none of these works has considered MAC layer link constraints, and therefore the bounds that they provide are optimistic. Moreover, because the interference range has significant effect on the performance of the MAC, and the interference range of CT could differ from that of a SISO transmission, a complete evaluation of multihop cooperative networks is hard to obtain without the inclusion of MAC constraints.

The challenges in modeling a multihop WSN reside mainly in the scope of modeling and the computational complexity. In this paper, we report a unified modeling framework with reasonable complexity that provides the fundamental performance limit that a MHWSN can theoretically achieve. Our model applies to both energy constrained networks and energy harvesting networks, and it captures the general constraints from MAC, routing, and energy dynamics. We further present two computational methods taking the advantage of the characteristics of EC networks and EH networks, respectively.

3. System Model and Assumptions

Each sensor node has two queues, the packet queue where data packets are stored and the energy queue where energy is stored. Both the packet queue and the energy queue have limited capacity. We assume that no dynamic transmit power control technique is used; that is, all the nodes use the same transmit power level, which means CT is used exclusively for range extension purpose instead of power reduction.

3.1. Time Slots in the System

Unlike [30] where the time slot is varied and is defined as the interval between two random events on a node, that is, between instances of energy arrivals, we consider a mini-time-slotted system as shown in Figure 2, where the slot duration is constant and slots are normalized to integral units $t\in \left\{\mathrm{1,2},3,\dots \right\}$ . By time t, we refer to the duration within $[t,t+1]$ . The randomness of events occurs during the slot. The reason to use constant slot interval is twofold: first, this provides the time from the perspective of a network that comprises many nodes; second, this allows us to model the packet transmission duration, against the instantaneous information capacity, because packetizing is very common in wireless standards, including IEEE 802.15.4, in which rate changing is infeasible in the duration of one-packet transmission.

In addition, we remark that the time slots in this paper are different and more general than those in the link scheduling literature. In the latter case, the time is divided into slots equal to the length of packet transmission time, and all nodes are synchronized to slot boundaries. Our model does not have such a constraint, so the packet length can span multiple slots and nodes are not necessarily synchronized to slot boundaries.

3.2. Topology and Link Models

3.2.1. Topology Model

The multihop network topology is modeled as a directed graph $G=(V,E)$ . $V=\{1,\dots ,N\}$ is the set of nodes except the Sink, whose ID is $0$ , and E is the set of links comprising both single-input-single-output (SISO) links and VMISO links; that is, $E=E^{\mathrm{s}\mathrm{o}}\cup E^{\mathrm{v}\mathrm{o}}$ . A SISO link l exists if its source and the destination are within the maximum transmission range. Given the transmission power $P_{tx}$ and the required receiving power $P_{rx,\min }$ , the maximum transmission range $d_{\max }$ is defined as $d_{\max }={\left(k{P}_{tx}/P_{rx,\mathrm{m}\mathrm{i}\mathrm{n}}\right)}^{1/\rho }$ , where ρ is the path loss exponent and k is the constant of proportionality [31]. A VMISO link exists if the source and the destination are within maximum transmission of CT with the cooperators. A VMISO link can be formed between the cooperating nodes (the initiator and cooperators) and the VMISO receiver. In the case of VMISO link, we assume the destination is the Sink; for example, Nodes A and B form a VMISO link to the Sink in Figure 1. Therefore, a link is represented by the 3-tuple, $〈s(l),d(l),\mathcal{C}\left(s(l)\right)〉\in E$ , where $s(l)$ and $d(l)$ are the source and the destination, respectively; $\mathcal{C}(s(l))$ is the cooperators of $s(l)$ . This link is a valid VMISO link if and only if

$$\begin{array}{c}{\left(10^{G_{\left|\mathcal{C}\left(s\left(l\right)\right)\right|}/10}\sum ^{i\in s\left(l\right)\cup \mathcal{C}\left(s\left(l\right)\right)}{D}_{i,d\left(l\right)}^{-\rho }\right)}^{-1/\rho }\le d_{\max },\end{array}$$

(1)

the same as in [21], where $D_{i,d(l)}$ is the distance between Node i and the VMISO receiver and $G_n$ is the SNR gain as a function of the number of cooperating nodes; for example, $G_2=10\hspace{0.17em}$ dB, $G_3=13.5\hspace{0.17em}$ dB according to [32] for BPSK modulation at a bit error rate of $10^{-3}$ . Note that, in the case of SISO link, $\mathcal{C}(s(l))=⌀$ and, in the case of VMISO link, $d(l)=0$ .

Figure 1

An illustration of interference model and VMISO link. Node v is the receiver of Node u. Node u's hidden nodes in the gray area will interfere with v's reception.

Figure 2

A timeline illustration of the decision process model for the network. The decision epochs are at the beginning of each minislot.

3.2.2. Interference Model

A node has three ranges associated with it: the transmission range ( $\mathrm{T}\mathrm{X}$ ), the interference range ( $\mathrm{I}\mathrm{F}$ ), and the carrier sensing range ( $\mathrm{C}\mathrm{S}$ ), as shown in Figure 1. Similar to [33], we assume that link $(i,j,·)$ interferes with link $(u,v,·)$ if either the distance $\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(i,v)\le \mathrm{I}\mathrm{F}(i)$ or $i=v$ . For instance, in Figure 1, Node u and any node in the gray area can start transmission at the same time because they are out of carrier sensing (CS) range of (i.e., hidden from) each other. Also, located in Node v's interference (IF) range, the hidden node's transmission interferes with Node v's reception. This phenomenon directly results from the CSMA mechanism of the MAC layer and the relationship between the three ranges (TX, IF, and CS) [34].

3.2.3. Medium Access Control Model

We assume that carrier-sensing-medium-access (CSMA) is performed before a node attempts to transmit. Further, similar to [35], it is assumed that (i) a node cannot transmit and receive at the same time; (ii) a node can transmit if none of its neighbors (in $\mathrm{C}\mathrm{S}$ range) is transmitting; (iii) link errors result only from collisions due to hidden terminals; (iv) nodes receive with perfect capture; that is, a packet is successfully decoded if the receiver and none of its neighbors are transmitting at the start of packet; and (v) the propagation time is zero. Our model can be extended to incorporate fading and shadowing effects of the wireless channel in computing packet errors; however, in this paper we assume (iii) to simplify the presentation.

3.3. Traffic and Energy Models

3.3.1. Traffic Model

The transmission duration of a link $l\in E$ is assumed exponentially distributed with the expectation ${\bar{T}}_l$ , the same as in [35, 36]. Note that this assumption is inaccurate; however, it is necessary to make the MDP problem tractable [36]. The memoryless property of the exponential distribution indicates that, given a packet is being transmitted at the beginning of a time slot, it completes the transmission within the slot of length $\Delta T$ with probability ${\beta }_l=1-e^{-\Delta T/{\bar{T}}_l}$ , regardless of the number of slots it has been transmitting. Making $\Delta T$ small enough allows us to assume that at most one link can finish transmission during a time slot; that is, the probability that more than one link finishes transmission is $o(\Delta t)$ .

The number of exogenous (self-generated) packets at Node i of commodity d (destined to Sink d), $f_i^d(t)$ , during any slot t is modeled as a binary random variable (RV), with probability mass functions (PMF) $\mathrm{Pr}[f_i^d(t)=1]={\alpha }^d(i)$ and $\mathrm{Pr}[f_i^d(t)=0]=1-{\alpha }^d(i)$ , where ${\alpha }^d(i)$ is assumed i.i.d. over time slots.

3.3.2. Energy Harvesting Model

A realistic power consumption model for a sensor node has been studied in [37], which decouples the power usage into baseband DSP circuit, the RF front-end circuit for transmitting and receiving, the power amplifier (PA), and low noise amplifier (LNA). A similar circuit-level analysis is also performed in [38] for energy modeling in cooperative multiple-input-multiple-output (MIMO) sensor networks. In [39], six different statistical models have been used to fit empirical datasets of a solar-powered energy harvesting sensor node. Our recent work [40] provides an energy model for energy harvesting node considering harvesting and leakage in a supercapacitor.

While we are aware of more complex energy models, for simplicity, in this paper we assume a simpler model. However, we note that it should be straightforward to incorporate complex energy consumption and harvesting models into our framework. In this study, the energy harvested at Node i during any slot t is modeled as a binary RV, which is denoted by $h_i(t)\in \{\mathrm{0,1}\}$ . This RV has PMF $\mathrm{Pr}[h_i(t)=1]=\gamma (i)$ and $\mathrm{Pr}[h_i(t)=0]=1-\gamma (i)$ . Representing the energy harvesting rate, $\gamma (i)$ is assumed i.i.d. over time slots. Thus, similar to [1, 41], the energy of Node i evolves as

$$\begin{array}{c}{\mathbf{e}}_i\left(t+1\right)=\min \left\{{\left({\mathbf{e}}_i\left(t\right)-e_{\text{com}}^{\left(i,l\right)}\right)}^{+}+h_i\left(t\right),e_{\max }\right\},\end{array}$$

(2)

where $e_{\text{com}}^{(i,l)}$ is the energy consumption of Node i if it participates in a finished link transmission l and $e_{\max }$ is the battery capacity of a node. We note that other EH models can also be applied, for example, the correlated EH model [2] where the energy harvested in a slot is correlated to the previous slot. Further, we assume that the energy harvested in a particular slot is available at the end of the slot.

4. Markov Decision Process Formulation

In this section, we describe the Markov Decision Process formulation for the EC and EH networks. We model the network state space by a tuple consisting of the transmission set in the network considering the MAC constraints, the queuing level, and the energy level for each node. Next we specify the action space at each time slot, which affects the network state. Then we analyze the system transition dynamics by decoupling the system into the tuple components with condition and then combining them based on the chain rule.

4.1. Energy Constrained (EC) Networks

4.1.1. Network State Space

The network state is defined as $\mathbf{s}\triangleq \{\mathcal{L},\mathbf{q},\mathbf{e}\}$ to include the transmission set $\mathcal{L}$ , the queuing level q of each node, and the energy e of each node. Note that the state space is not the Cartesian product of each component's state space, because they are interactive (e.g., an active link cannot have an empty queue at its source node).

The transmission set $\mathcal{L}(t)$ includes the links that are active (in transmission) at time t, which must not violate the carrier sensing constraints of the MAC. Note that it does not mean that all the receivers of links in the transmission set will receive successfully. Denote the state space of $\mathcal{L}$ as $\mathbb{L}$ . This component includes the collisions in the MAC layer. We denote the links that are free of collision as $\Phi (\mathcal{L})\subset \mathcal{L}$ . $\Phi (\mathcal{L})$ is deterministic given $\mathcal{L}$ , and $\Phi ={\Phi }^{\mathrm{s}\mathrm{o}}\cup {\Phi }^{\mathrm{v}\mathrm{o}}$ .

To find $\mathbb{L}$ , it is equivalent to find the matchings of a graph. Given a graph $G=(V,E)$ , a matching $M(G)$ in G is a set of nonadjacent edges (i.e., no two edges share a common vertex), implying $\mathcal{L}\in \{M(G)\}$ . Further, the MAC layer CSMA adds constraints on $\mathcal{L}$ , where for any two links $l_1,l_2\in \mathcal{L}$ , the sources of the two links ( $s(l_1),s(l_2)$ ) must be out of carrier sensing (CS) range of each other. Therefore, $\mathbb{L}$ is determined by graph G and the binary hearing matrix $\mathbf{H}=[h_{ij}{]}_{i,j\in V}$ of the network, where the element $h_{i,j}=1$ if and only if (iff) Node i and Node j are within CS range of each other. Though enumerating matchings of a graph is NP-complete, Bron-Kerbosch algorithm has been shown to be one of the fastest [42] and has been used in this paper.

4.1.2. Decision Epochs and Action Space

A decision epoch corresponds to the beginning of a time slot. The set of decision epochs are denoted by $\{\mathrm{1,2},\dots ,T\}$ . When T is finite (infinite), the decision problem is referred to as a finite-horizon (infinite-horizon) problem. An action, $a(t)$ , is to admit a new link to the “remaining” transmission set $\mathcal{L}(t)$ , which comprises those links that did not complete transmission during slot $t-1$ . Note that $\mathcal{L}(t)\subset \mathcal{L}(t-1)\hspace{0.17em}\cup \hspace{0.17em}a(t-1)$ . The action space at time t when the system state is $\mathbf{s}(t)$ is represented by

$$\begin{array}{c}\mathcal{A}\left(\mathbf{s}\left(t\right)\right)=\left\{l\in E:{\mathbf{q}}_{s\left(l\right)}\left(t\right)>0,\hspace{0.17em}\hspace{0.17em}l\notin \mathcal{L}\left(t\right),\hspace{0.17em}\hspace{0.17em}\mathcal{L}\left(t\right)\cup l\subset \mathbb{L},E_{\mathrm{s}\mathrm{u}\mathrm{f}}\right\},\end{array}$$

(3)

where $E_{\mathrm{s}\mathrm{u}\mathrm{f}}$ indicates ${\mathbf{e}}_i\ge e_{\text{com}}^{(i,l)}$ , $\forall i\in l$ . Therefore, a link l can join $\mathcal{L}(t)$ iff its source has a nonempty queue, $\mathcal{L}(t)\cup l$ is also a transmission set, and the participating nodes have enough energy. Note that it is assumed that at most one link can be admitted to $\mathcal{L}(t)$ , because $\Delta T$ is small. Also, note that the null set $⌀\in \mathcal{A}(\mathbf{s}(t))$ . As a result, an action can be either a “CT” link, a “non-CT” link, or null (no new transmission).

4.1.3. State Transition Dynamics

We first characterize the state transition of each component.

( 1) Transmission Set Dynamics. Let action $a(t)\in \mathcal{A}(\mathbf{s}(t))$ denote the link admitted to the transmission set. Denote $p\left\{z^l(t)\mid \mathcal{L}(t),a(t)\right\}$ as the probability that only link $l\in \mathcal{L}(t)\cup a(t)$ finishes transmission during slot t and $p\left\{z^{⌀}(t)\mid \mathcal{L}(t),a(t)\right\}$ as the probability that no link finishes transmission during slot t. The time index is dropped when there is no ambiguity. For abbreviation, we denote the transition kernel as $p\left\{{\mathcal{L}}^{\prime }\mid \mathcal{L},a\right\}:=p\left\{\mathcal{L}\left(t+1\right)\mid \mathcal{L}\left(t\right),a\left(t\right)\right\}$ ,

$$\begin{array}{c}p\left\{z^l\mid \mathcal{L},a\right\}={\beta }_l{\displaystyle \prod }_{\begin{array}{c}\underset{k\ne l}{k\in \mathcal{L}\cup a}\end{array}}\left(1-{\beta }_k\right),\hspace{1em}\forall l\in \mathcal{L}\cup a,\\ p\left\{z^{⌀}\mid \mathcal{L},a\right\}=\left\{\begin{array}{cc}{\displaystyle \prod }_{l\in \mathcal{L}\cup a}\left(1-{\beta }_l\right)\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}\mathcal{L}\cup a\ne ⌀,\hfill \\ 1\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}\mathcal{L}\cup a=⌀.\hfill \end{array}\right.\end{array}$$

(4)

Define ${\mathcal{D}}^{(L)}\triangleq \left\{\mathcal{L}\cup a\right\}\setminus {\mathcal{L}}^{\prime }$ ; then we can get

$$\begin{array}{c}p\left\{{\mathcal{L}}^{\prime }\mid \mathcal{L},a\right\}=\left\{\begin{array}{cc}p\left\{z^l\mid \mathcal{L},a\right\}\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{\mathcal{D}}^{\left(L\right)}=l,\hfill \\ p\left\{z^{⌀}\mid \mathcal{L},a\right\}\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{\mathcal{D}}^{\left(L\right)}=⌀,\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}\right.\end{array}$$

(5)

( 2) Queue Length Dynamics. Considering that EC networks are most suitable for light traffic applications, such as monitoring, we assume that at most one node self-generates a packet during a slot. Let ${\mathbf{I}}_i$ represent a vector with its ith element being $1$ and other elements being $0$ , $1\le i\le N$ , and let ${\mathbf{I}}_0$ be the zero vector. The probability that only Node i or no node ( $i=0$ ) self-generates a packet is given by

$$\begin{array}{c}P\left\{{\mathbf{I}}_i\right\}=\left\{\begin{array}{cc}\alpha \left(i\right){\displaystyle \prod }_{\begin{array}{c}j\in V,j\ne i\end{array}}\left(1-\alpha \left(j\right)\right)\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}i\in V,\hfill \\ {\displaystyle \prod }_{\begin{array}{c}j\in V\end{array}}\left(1-\alpha \left(j\right)\right)\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}i=0,\hfill \end{array}\right.\end{array}$$

(6)

where $\alpha (i)$ is the probability of a self-generated packet arrival for Node i in a time slot and is assumed i.i.d. over time slots.

Then the queue evolution can be described by the following cases.

Case 1 (( $c_1$ ): ${\mathbf{q}}^{\prime }=\mathbf{q}$ ).

This can happen for two reasons: ( $c_{1a}$ ) no link transmission affects the queue state and no new exogenous arrival affects the queues and ( $c_{1b}$ ) some link $(i,0,·)$ destined to the Sink successfully transmits and the source self-generates a new packet:

$$\begin{array}{c}{P}_{c_{1a}}=\left(\mathbb{1}\left\{\Delta \mathcal{L}=⌀\right\}+\mathbb{1}\left\{\Delta \mathcal{L}\in \bar{\Psi }\left(\mathcal{L}\cup a\right)\right\}+\mathbb{1}\left\{\Delta \mathcal{L}\in \Psi \left(\mathcal{L}\cup a\right),q_{d\left(\Delta \mathcal{L}\right)}=q^{\max }\right\}\right)·\left(P\left\{{\mathbf{I}}_0\right\}+{\displaystyle \sum }_{q_k=q^{\max },k\in V}p\left\{{\mathbf{I}}_k\right\}\right),\\ P_{c_{1b}}=\mathbb{1}\left\{\Delta \mathcal{L}\in \Psi \left(\mathcal{L}\cup a\right),s\left(\Delta \mathcal{L}\right)=i,d\left(\Delta \mathcal{L}\right)=0\right\}·P\left\{{\mathbf{I}}_i\right\}.\end{array}$$

(7)

Case 2 (( $c_2$ ): ${\mathbf{q}}^{\prime }=\mathbf{q}+{\mathbf{I}}_j$ , $j\in V$ ).

This case includes two possibilities: ( $c_{2a}$ ) some SISO link $(·,j,⌀)$ that is directed to Node j successfully transmits and the source self-generates a new packet and ( $c_{2b}$ ) no link transmission affects the queue state and there is a new exogenous packet arrival to Node j:

$$\begin{array}{c}{P}_{c_{2a}}\left(j\right)=\mathbb{1}\left\{\Delta \mathcal{L}\in \Psi \left(\mathcal{L}\cup a\right),d\left(\Delta \mathcal{L}\right)=j\right\}·P\left\{{\mathbf{I}}_{s\left(\Delta \mathcal{L}\right)}\right\},\\ P_{c_{2b}}\left(j\right)=\left(\mathbb{1}\left\{\Delta \mathcal{L}=⌀\right\}+\mathbb{1}\left\{\Delta \mathcal{L}\in \bar{\Psi }\left(\mathcal{L}\cup a\right)\right\}+\mathbb{1}\left\{\Delta \mathcal{L}\in \Psi \left(\mathcal{L}\cup a\right),q_{d\left(\Delta \mathcal{L}\right)}=q^{\mathrm{m}\mathrm{a}\mathrm{x}}\right\}\right)·P\left\{{\mathbf{I}}_j\right\}.\end{array}$$

(8)

Case 3 (( $c_3$ ): ${\mathbf{q}}^{\prime }=\mathbf{q}-{\mathbf{I}}_i$ , $\hspace{0.17em}\hspace{0.17em}i\in V$ ).

This happens because some link $(i,0,·)$ destined to the Sink successfully transmits and no new exogenous arrival affects the queues:

$$\begin{array}{c}{P}_{c_3}\left(i\right)=\mathbb{1}\left\{\Delta \mathcal{L}\in \Psi \left(\mathcal{L}\cup a\right),\hspace{0.17em}\hspace{0.17em}\Delta \mathcal{L}=\left(i,0,·\right)\right\}·\left(P\left\{{\mathbf{I}}_0\right\}+{\displaystyle \sum }_{q_k^{\prime }=q^{\max },k\in V}P\left\{{\mathbf{I}}_k\right\}\right).\end{array}$$

(9)

Case 4 (( $c_4$ ): ${\mathbf{q}}^{\prime }=\mathbf{q}-{\mathbf{I}}_i+{\mathbf{I}}_j$ , $\hspace{0.17em}\hspace{0.17em}i,j\in V$ , $\hspace{0.17em}\hspace{0.17em}i\ne j$ ).

This case includes two possibilities: ( $c_{4a}$ ) some SISO link $(i,j,⌀)$ successfully transmits and there was no change in queue due to new exogenous arrival and ( $c_{4b}$ ) link $(i,0,·)$ successfully transmits and Node j self-generates a new packet:

$$\begin{array}{c}{P}_{c_{4a}}\left(i,j\right)=\mathbb{1}\left\{\Delta \mathcal{L}\in \Psi \left(\mathcal{L}\cup a\right),\hspace{0.17em}\hspace{0.17em}\Delta \mathcal{L}=\left(i,j,⌀\right)\right\}·\left(P\left\{{\mathbf{I}}_0\right\}+{\displaystyle \sum }_{q_k^{\prime }=q^{\max },k\in V}P\left\{{\mathbf{I}}_k\right\}\right),\\ P_{c_{4b}}\left(i,j\right)=\mathbb{1}\left\{\Delta \mathcal{L}\in \Psi \left(\mathcal{L}\cup a\right),\hspace{0.17em}\hspace{0.17em}\Delta \mathcal{L}=\left(i,0,·\right)\right\}·P\left\{{\mathbf{I}}_j\right\}.\end{array}$$

(10)

Case 5 (( $c_5$ ): ${\mathbf{q}}^{\prime }=\mathbf{q}-{\mathbf{I}}_i+2{\mathbf{I}}_j$ , $\hspace{0.17em}\hspace{0.17em}i,j\in V$ , $\hspace{0.17em}\hspace{0.17em}i\ne j$ ).

This occurs because the SISO link $(i,j,⌀)$ successfully transmits and Node j self-generates a new packet:

$$\begin{array}{c}{P}_{c_5}\left(i,j\right)=\mathbb{1}\left\{\Delta \mathcal{L}\in \Psi \left(\mathcal{L}\cup a\right),\hspace{0.17em}\hspace{0.17em}\Delta \mathcal{L}=\left(i,j,⌀\right)\right\}·P\left\{{\mathbf{I}}_j\right\}.\end{array}$$

(11)

Case 6 (( $c_6$ ): ${\mathbf{q}}^{\prime }=\mathbf{q}-{\mathbf{I}}_i+{\mathbf{I}}_j+{\mathbf{I}}_k$ , $\hspace{0.17em}\hspace{0.17em}i,j,k\in V$ , $\hspace{0.17em}\hspace{0.17em}i\ne j\ne k$ ).

This occurs because the SISO link $(i,j,⌀)$ successfully transmits and a different Node k self-generates a new packet:

$$\begin{array}{c}{P}_{c_6}\left(i,j,k\right)=\mathbb{1}\left\{\Delta \mathcal{L}\in \Psi \left(\mathcal{L}\cup a\right),\hspace{0.17em}\hspace{0.17em}\Delta \mathcal{L}=\left(i,j,⌀\right)\right\}·P\left\{{\mathbf{I}}_k\right\}+\mathbb{1}\left\{\Delta \mathcal{L}\in \Psi \left(\mathcal{L}\cup a\right),\hspace{0.17em}\hspace{0.17em}\Delta \mathcal{L}=\left(i,k,⌀\right)\right\}·P\left\{{\mathbf{I}}_j\right\}.\end{array}$$

(12)

Therefore, the transition kernel of $\mathbf{q}(t)$ is expressed as (23).

( 3) Energy Evolution Dynamics. The process $\mathbf{e}(t)$ is dictated by transmission energy and receiving energy consumption incurred during a finished link transmission, regardless of success or failure. In the CT case, the cooperators consume additional energy in receiving the broadcast packet initiated by the source and in conducting CT. The transition kernel of $\mathbf{e}(t)$ is given in (25).

( 4) System Dynamics. The transition matrix of the system states $\mathbf{s}\triangleq \{\mathcal{L},\mathbf{e},\mathbf{q}\}$ can be obtained by the following theorem.

Theorem 1.

The transition kernel of the EC system, $p\{{\mathbf{s}}^{\prime }\mid \mathbf{s},a\}$ , is equal to the product of (5), (13), and (14):

$$\begin{array}{c}P\left\{{\mathbf{q}}^{\prime }\mid {\mathcal{L}}^{\prime },\mathcal{L},\mathbf{q},a\right\}=\left\{\begin{array}{cc}{P}_{c_{1a}}+P_{c_{1b}}\hfill & \mathit{\text{if\hspace{0.17em}\hspace{0.17em}}}{\mathbf{q}}^{\prime }=\mathbf{q},\hfill \\ P_{c_{2a}}\left(j\right)+P_{c_{2b}}\left(j\right)\hfill & \mathit{\text{if\hspace{0.17em}\hspace{0.17em}}}{\mathbf{q}}^{\prime }=\mathbf{q}+{\mathbf{I}}_j,\hspace{0.17em}\hspace{0.17em}j\in V,\hfill \\ P_{c_3}\left(i\right)\hfill & \mathit{\text{if\hspace{0.17em}\hspace{0.17em}}}{\mathbf{q}}^{\prime }=\mathbf{q}-{\mathbf{I}}_i,\hspace{0.17em}\hspace{0.17em}i\in V,\hfill \\ P_{c_{4a}}\left(i,j\right)+P_{c_{4b}}\left(i,j\right)\hfill & \mathit{\text{if\hspace{0.17em}\hspace{0.17em}}}{\mathbf{q}}^{\prime }=\mathbf{q}-{\mathbf{I}}_i+{\mathbf{I}}_j,\hspace{0.17em}\hspace{0.17em}i,j\in V,\hfill \\ P_{c_5}\left(i,j\right)\hfill & \mathit{\text{if\hspace{0.17em}\hspace{0.17em}}}{\mathbf{q}}^{\prime }=\mathbf{q}-{\mathbf{I}}_i+2{\mathbf{I}}_j,\hspace{0.17em}\hspace{0.17em}i,j\in V,\hspace{0.17em}\hspace{0.17em}i\ne j,\hfill \\ P_{c_6}\left(i,j,k\right)\hfill & \mathit{\text{if\hspace{0.17em}\hspace{0.17em}}}{\mathbf{q}}^{\prime }=\mathbf{q}-{\mathbf{I}}_i+{\mathbf{I}}_j+{\mathbf{I}}_k,\hspace{0.17em}\hspace{0.17em}i,j,k\in V,\hspace{0.17em}\hspace{0.17em}i\ne j\ne k,\hfill \\ 0\hfill & \mathit{\text{otherwise}},\hfill \end{array}\right.\end{array}$$

(13)

$$\begin{array}{c}P\left\{{\mathbf{e}}^{\prime }\mid {\mathcal{L}}^{\prime },\mathcal{L},\mathbf{e},a\right\}=\left\{\begin{array}{cc}\mathbb{1}\left\{{\mathbf{e}}^{\prime }=\mathbf{e}-{\mathbf{I}}_i{e}_{tx}-{\mathbf{I}}_j{e}_{rx}\right\}\hfill & if\Delta \mathcal{L}=\left(i,j,⌀\right)\in {\left(\mathcal{L}\cup a\right)}^{\mathit{\text{so}}},\hspace{0.17em}\hspace{0.17em}j\in V^{\ast },\hfill \\ \mathbb{1}\left\{{\mathbf{e}}^{\prime }=\mathbf{e}-{\mathbf{I}}_i{e}_{\mathit{\text{init}}}^{\mathit{\text{CT}}}-{\displaystyle \sum }_{k\in \mathcal{H}\left(i\right)}{\mathbf{I}}_k{e}_{co}^{\mathit{\text{CT}}}\right\}\hfill & if\Delta \mathcal{L}=\left(i,0,\mathcal{H}\left(i\right)\right)\in {\left(\mathcal{L}\cup a\right)}^{\mathit{\text{vo}}},\hspace{0.17em}\hspace{0.17em}i\in V,\hfill \\ \mathbb{1}\left\{{\mathbf{e}}^{\prime }=\mathbf{e}\right\}\hfill & if\Delta \mathcal{L}=⌀,\hfill \\ 0\hfill & otherwise.\hfill \end{array}\right.\end{array}$$

(14)

Proof.

According to the chain rule, $p\{{\mathbf{s}}^{\prime }\mid \mathbf{s},a\}$ can be expressed as $p({\mathbf{q}}^{\prime }\mid A)·p({\mathbf{e}}^{\prime }\mid B)·p({\mathcal{L}}^{\prime }\mid C)$ , where $A=\{({\mathcal{L}}^{\prime },{\mathbf{e}}^{\prime },\mathcal{L},\mathbf{q},\mathbf{e},a)\}$ , $B=\{({\mathcal{L}}^{\prime },\mathcal{L},\mathbf{q},\mathbf{e},a)\}$ , and $C=\{(\mathcal{L},\mathbf{q},\mathbf{e},a)\}$ . Further, since ${\mathbf{q}}^{\prime }$ and $({\mathbf{e}}^{\prime },\mathbf{e})$ are conditionally independent, given $({\mathcal{L}}^{\prime },\mathcal{L},\mathbf{q},a)$ , the first element in the product, $p({\mathbf{q}}^{\prime }\mid A)$ , is reduced to (13). Applying similar arguments to $p({\mathbf{e}}^{\prime }\mid B)$ and $p({\mathcal{L}}^{\prime }\mid C)$ , Theorem 1 is proved.

4.1.4. Expected Total Rewards

During a time slot, the system obtains a reward $g(\mathbf{s},a)=1$ if a packet was delivered to the Sink, and $g(\mathbf{s},a)=0$ otherwise, where s is the current state. Then, with the state space S and the termination states $S_t$ , for $\mathbf{s}\in S\setminus S_t$ , we have

$$\begin{array}{c}g\left(\mathbf{s},a\right)\triangleq E\left[r\right]={\displaystyle \sum }_{\begin{array}{c}\underset{d\left(l\right)=0}{l\in \Psi \left(\mathcal{L}\cup a\right)}\end{array}}P\left\{z^l\mid \mathcal{L},a\right\}\mathbb{1}\left\{\Psi \left(\mathcal{L}\cup a\right)\ne ⌀\right\}.\end{array}$$

(15)

And, for $\mathbf{s}\in S_t$ , $g(\mathbf{s},a)=0$ . The expected total rewards of the process starting from an initial state s, under the policy π, are denoted by

$$\begin{array}{c}{\mathcal{J}}^{\mathbf{\pi }}\left(\mathbf{s}\right)={\displaystyle \sum }_{t=1}^{\infty }g\left(\mathbf{s}\left(t\right),a\left(t\right)\right).\end{array}$$

(16)

4.1.5. MDP Formulation

A transmission policy is a series of decision rules $\mathbf{\pi }=[a(1),a(2),\dots ]$ , where $a(t):\mathcal{S}(t)\to \mathcal{A}\left(\mathcal{S}(t)\right)$ . The maximum lifetime ${\mathcal{J}}^{\ast }(\mathbf{s})$ is given by

$$\begin{array}{c}{\mathcal{J}}^{\ast }\left(\mathbf{s}\right)=\underset{\mathbf{\pi }}{\max }{\mathcal{J}}^{\mathbf{\pi }}\left(\mathbf{s}\right).\end{array}$$

(17)

The optimal lifetime ${\mathcal{J}}^{\ast }(\mathbf{s})$ is the unique solution of the Bellman's optimality equation [16]:

$$\begin{array}{c}\mathcal{J}\left(\mathbf{s}\right)=\underset{a}{\max }\left\{g\left(\mathbf{s},a\right)+{\displaystyle \sum }_{{\mathbf{s}}^{\prime }}P\left\{{\mathbf{s}}^{\prime }\mid \mathbf{s},a\right\}\mathcal{J}\left({\mathbf{s}}^{\prime }\right)\right\}.\end{array}$$

(18)

A policy ${\mathbf{\pi }}^{\ast }$ is optimal if it achieves the maximum expected lifetime for all starting states; that is,

$$\begin{array}{c}{\mathcal{J}}^{{\mathit{\pi }}^{\ast }}\left(\mathbf{s}\right)={\mathcal{J}}^{\ast }\left(\mathbf{s}\right),\hspace{1em}\forall \mathbf{s}\in \mathcal{S}\setminus {\mathcal{S}}_t.\end{array}$$

(19)

4.2. Energy Harvesting Networks

The link set dynamics in the EH network are the same as in the EC network. In this section, we give a more general expression for the queue system, relaxing the low-traffic assumption in the previous sections. EH networks may provide better support for high traffic load than EC networks [43].

In a uniform form, the evolutions of q and e are governed by the following balance equations:

$$\begin{array}{c}{\psi }^{\left(i\right)}\left(t+1\right)={\psi }^{\left(i\right)}\left(t\right)+{\mathbf{R}}^{\left(i\right)}{\mathbf{M}}^{\left(i\right)}\left(t\right){\Gamma }^{\left(i\right)}\left(t\right)+{\sigma }^{\left(i\right)}\left(t\right),\end{array}$$

(20)

where $i\in \{q,e\}$ , ${\psi }^{(q)}=\mathbf{q}$ , and ${\psi }^{(e)}=\mathbf{e}$ ; and ${\mathbf{R}}^{(i)}$ and ${\mathbf{M}}^{(i)}$ will be defined below.

4.2.1. Queue Length Dynamics

In (20), ${\mathbf{M}}^{(q)}$ is a $\left|E\right|\times \left|E\right|$ diagonal matrix. The diagonal element ${\mathbf{M}}_l^{(q)}=1$ if the transmission of link l is completed and successful (a finished transmission, while always consuming energy, is not necessarily successfully decoded by the receiver, because of collision due to hidden terminals; in case of failure, the packet remains in the queue of the transmitter), and ${\mathbf{M}}_l^{(q)}=0$ otherwise. Note that it is assumed ${\sum }_{l\in E}‍{\mathbf{M}}_l^{(q)}\le 1$ ; that is, at most one link can finish transmission. ${\Gamma }^{(q)}$ is the transmission schedule vector, which satisfies ${\Gamma }_l^{(q)}=1$ , $\forall l\in \mathcal{L}\cup a$ , and ${\Gamma }_l^{(q)}=0$ otherwise. ${\mathbf{R}}^{(q)}$ is the $N\times \left|E\right|$ routing matrix, whose element in the ith row and lth column is

$$\begin{array}{c}{r}_{il}^{\left(q\right)}=\left\{\begin{array}{cc}1\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}d\left(l\right)=i\text{\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}Node\hspace{0.17em}\hspace{0.17em}}i\text{\hspace{0.17em}\hspace{0.17em}is\hspace{0.17em}\hspace{0.17em}not\hspace{0.17em}\hspace{0.17em}Sink},\hfill \\ -1\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}s\left(l\right)=i,\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}\right.\end{array}$$

(21)

${\sigma }^{(q)}$ is an $N\times 1$ vector representing the number of self-generated packets and ${\sigma }_i^{(q)}=f_i$ .

Let ${\mathbf{I}}_i$ represent a vector with its ith element being $1$ and other elements being $0$ , $1\le i\le N$ , and let ${\mathbf{I}}_0$ be the zero vector. Denote $V^{(q)}(t)$ as the nodes whose queue is not full at time t and define the difference

$$\begin{array}{c}{\mathcal{D}}^{\left(q\right)}=\left\{\begin{array}{cc}{\mathbf{q}}^{\prime }-\mathbf{q}\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{\mathcal{D}}^{\left(L\right)}=⌀,\hfill \\ {\mathbf{q}}^{\prime }-\mathbf{q}+{\mathbf{I}}_{s\left({\mathcal{D}}^{\left(L\right)}\right)}-{\mathbf{I}}_{d\left({\mathcal{D}}^{\left(L\right)}\right)}\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{\mathcal{D}}^{\left(L\right)}\in {\Phi }^{\text{so}}\left\{\mathcal{L}\cup a\right\},\hfill \\ {\mathbf{q}}^{\prime }-\mathbf{q}+{\mathbf{I}}_{s\left({\mathcal{D}}^{\left(L\right)}\right)}\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{\mathcal{D}}^{\left(L\right)}\in {\Phi }^{\text{vo}}\left\{\mathcal{L}\cup a\right\}.\hfill \end{array}\right.\end{array}$$

(22)

Then considering exogenous packet arrivals, the queue status is updated as follows:

$$\begin{array}{c}p\left\{{\mathbf{q}}^{\prime }\mid {\mathcal{L}}^{\prime },\mathcal{L},\mathbf{q},a\right\}=\left\{\begin{array}{cc}{\displaystyle \prod }_{i\in V}\left({\alpha }_i\mathbb{1}\left\{{\mathcal{D}}_i^{\left(q\right)}=1\right\}+\left(1-{\alpha }_i\right)\mathbb{1}\left\{{\mathcal{D}}_i^{\left(q\right)}=0,\hspace{0.17em}\hspace{0.17em}i\in V^{\left(q\right)}\right\}\right)\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{\mathcal{D}}^{\left(L\right)}\in \left\{{\Phi }^{\text{so}}\left\{\mathcal{L}\cup a\right\},{\Phi }^{\text{vo}}\left\{\mathcal{L}\cup a\right\},⌀\right\},\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}\right.\end{array}$$

(23)

where $\mathbb{1}\left\{·\right\}$ is the indicator function.

4.2.2. Energy Evolutions

Using analysis similar to that for the queue in Section 4.2.1, denote $V^{(e)}(t)$ as the nodes whose battery is not full at time t and define

$$\begin{array}{c}{\mathcal{D}}^{\left(e\right)}=\left\{\begin{array}{cc}{\mathbf{e}}^{\prime }-\mathbf{e}\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{\mathcal{D}}^{\left(L\right)}=⌀,\hfill \\ {\mathbf{e}}^{\prime }-\mathbf{e}+{\mathbf{I}}_{s\left({\mathcal{D}}^{\left(L\right)}\right)}·e_{tx}+{\mathbf{I}}_{d\left({\mathcal{D}}^{\left(L\right)}\right)}·e_{rx}\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{\mathcal{D}}^{\left(L\right)}\in E^{\mathrm{s}\mathrm{o}},\hfill \\ {\mathbf{e}}^{\prime }-\mathbf{e}+{\mathbf{I}}_{s\left({\mathcal{D}}^{\left(L\right)}\right)}{e}_{\text{init}}^{\mathrm{C}\mathrm{T}}+{\displaystyle \sum }_{k\in \mathcal{C}\left({\mathcal{D}}^{\left(L\right)}\right)}{\mathbf{I}}_k{e}_{\mathrm{co}}^{\mathrm{C}\mathrm{T}}\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{\mathcal{D}}^{\left(L\right)}\in E^{\mathrm{v}\mathrm{o}}.\hfill \end{array}\right.\end{array}$$

(24)

Then, considering the influence of energy harvesting, we get

$$\begin{array}{c}p\left\{{\mathbf{e}}^{\prime }\mid {\mathcal{L}}^{\prime },\mathcal{L},\mathbf{e},a\right\}=\left\{\begin{array}{cc}{\displaystyle \prod }_{i\in V}\left(h_i\mathbb{1}\left\{{\mathcal{D}}_i^{\left(e\right)}=1\right\}+\left(1-h_{\mathrm{i}}\right)\mathbb{1}\left\{{\mathcal{D}}_i^{\left(e\right)}=0,\hspace{0.17em}\hspace{0.17em}i\in V^{\left(e\right)}\right\}\right)\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{\mathcal{D}}^{\left(L\right)}\in \left\{E^{\text{so}},E^{\text{vo}},⌀\right\},\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}\right.\end{array}$$

(25)

4.2.3. System Dynamics

Similar to the EC network, the transition matrix of the system states $\mathbf{s}\triangleq \{\mathcal{L},\mathbf{e},\mathbf{q}\}$ for EH network can be obtained by the following theorem.

Theorem 2.

The transition kernel of the EH system, $p\{{\mathbf{s}}^{\prime }\mid \mathbf{s},a\}$ , is equal to the product of (5), (23), and (25).

4.3. Performance Criterion and Reward Function

Let $v^{\pi }(\mathbf{s})$ denote the expected total reward, given an initial state s and a policy π (a series of decisions),

$$\begin{array}{c}{v}^{\pi }\left(\mathbf{s}\right)=\underset{T\to \infty }{\lim }{\mathbb{E}}_{\mathbf{s}}^{\pi }\left[{\displaystyle \sum }_{t=1}^{T}r\left(\mathbf{s}\left(t\right),a\left(t\right)\right)\right],\end{array}$$

(26)

where T is the decision-making horizon length.

The definition of the reward function, $r(\mathbf{s},a)$ , depends on application requirements and can be subjective to network designers. We define two reward functions as follows:

$$\begin{array}{c}{r}^{\left(1\right)}\left(\mathbf{s},a\right)={\displaystyle \sum }_{\begin{array}{c}\underset{d\left(l\right)=0}{l\in \Phi \left(\mathcal{L}\cup a\right)}\end{array}}P\left\{z^l\mid \mathcal{L},a\right\}\mathbb{1}\left\{\Phi \left(\mathcal{L}\cup a\right)\ne ⌀\right\},\end{array}$$

(27)

$$\begin{array}{c}{r}^{\left(2\right)}\left(\mathbf{s},a\right)={\omega }_1{r}^{\left(1\right)}\left(\mathbf{s},a\right)-{\omega }_2\mathbb{1}\left\{\mathbf{s}\in S_D\right\}.\end{array}$$

(28)

In (27), a unit of reward is obtained if a packet is successfully received by the Sink. Equation (28) is a weighted sum of the rewards from delivered packets and the penalty from a QoS degradation (e.g., $S_D$ are states where any node's energy drops to zero).

Since the state space and the action space are finite, with a finite reward function, (26) can be expressed as

$$\begin{array}{c}{v}_{\lambda }^{\pi }\left(s\right)={\mathbb{E}}_s^{\pi }\left\{{\displaystyle \sum }_{t=1}^{\infty }{\lambda }^{t-1}r\left(\mathbf{s}\left(t\right),a\left(t\right)\right)\right\},\end{array}$$

(29)

where it is assumed that T is geometrically distributed with parameter λ, $0\le \lambda <1$ . Thus, the expected value of $\hspace{0.17em}\hspace{0.17em}T$ is $1/(1-\lambda )$ . The parameter λ is the discount factor, which measures the present value of one unit of reward received one period in the future [16].

4.4. Optimality Equations

The optimality equation for the expected total discounted reward criteria given an initial state, a.k.a. Bellman's equation, is given as follows:

$$\begin{array}{c}{v}_{\lambda }^{\ast }\left(\mathbf{s}\right)=\underset{a\in \mathcal{A}\left(s\right)}{\max }\left\{r\left(\mathbf{s},a\right)+{\displaystyle \sum }_{{\mathbf{s}}^{\prime }}\lambda p\left\{{\mathbf{s}}^{\prime }\mid \mathbf{s},a\right\}{v}_{\lambda }^{\ast }\left({\mathbf{s}}^{\prime }\right)\right\}.\end{array}$$

(30)

5. Computational Methods

In this section, we present the offline computational algorithms for the EC and EH networks. The methods seek to find the global optimum and do not require any information exchange among the nodes.

5.1. Energy Constrained Networks

Considering that the total energies are nonincreasing, we group the states according to the sum energy (in decreasing order) and solve the stages backwards. The transmission sets are computed using the Bron-Kerbosch algorithm [42], which is widely used and considered as one of the fastest. In each stage, the computation of the energy and queue state spaces is similar to the Bin-Ball problem (the Bin-Ball problem refers to enumerating the ways of allocating $n_1$ balls into $n_2$ bins. In our problem, the “bins” are the nodes and the “balls” are the total energies or queued packets). For instance, the number of ways to allocate $n_1$ units of energies into $n_2$ nodes with minimum energy requirement of $\theta =e_{\text{th}}+1$ (in $\mathcal{S}\setminus {\mathcal{S}}_t$ ) is given by

$$\begin{array}{c}{\displaystyle \sum }_{l=0}^{n_2}{\left(-1\right)}^l\left(\begin{array}{c}{n}_2\\ l\end{array}\right)\left(\begin{array}{c}{n}_2+n_1-n_2\theta -l\left(e^{\max }+1-\theta \right)-1\\ n_2-1\end{array}\right).\end{array}$$

(31)

For an integer m (stage), the set of states $S_m$ is defined as

$$\begin{array}{c}{S}_m=\left\{\left(\mathcal{L},\mathbf{e},\mathbf{q}\right)\mid {\displaystyle \sum }_{i\in V}{\mathbf{e}}_i=m\right\},\hspace{1em}m=N\theta ,\dots ,N{e}^{\max }.\end{array}$$

(32)

For each $S_m$ and $\mathbf{s}\in S_m$ , we have

$$\begin{array}{c}{\mathcal{J}}^{\ast }\left(\mathbf{s}\right)=\underset{a\in \mathcal{A}\left(\mathbf{s}\right)}{\max }\left\{g\left(\mathbf{s},a\right)+{\displaystyle \sum }_{{\mathbf{s}}^{\prime }}P\left\{{\mathbf{s}}^{\prime }\mid \mathbf{s},a\right\}{\mathcal{J}}^{\ast }\left({\mathbf{s}}^{\prime }\right)\right\}=\underset{a\in \mathcal{A}\left(\mathbf{s}\right)}{\max }\left\{g\left(\mathbf{s},a\right)+{\displaystyle \sum }_{{\mathbf{s}}^{\prime }\in S_{N\theta }\cup \cdots \cup S_{m-1}}P\left\{{\mathbf{s}}^{\prime }\mid \mathbf{s},a\right\}{\mathcal{J}}^{\ast }\left({\mathbf{s}}^{\prime }\right)+{\displaystyle \sum }_{{\mathbf{s}}^{\prime }\in S_m}P\left\{{\mathbf{s}}^{\prime }\mid \mathbf{s},a\right\}{\mathcal{J}}^{\ast }\left({\mathbf{s}}^{\prime }\right)\right\},\hspace{1em}\mathbf{s}\in S_m.\end{array}$$

(33)

Therefore, we formulate a mixed integer linear programming (MILP) problem (Subproblem-1):

$$\begin{array}{c}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}\hspace{1em}{\displaystyle \sum }_{\mathbf{s}\in S_m}\mathcal{J}\left(\mathbf{s}\right)\\ \mathrm{s}.\mathrm{t}.\hspace{1em}\mathcal{J}\left(\mathbf{s}\right)-{\displaystyle \sum }_{{\mathbf{s}}^{\prime }\in S_m}P\left({\mathbf{s}}^{\prime }\mid \mathbf{s},a\right)\mathcal{J}\left({\mathbf{s}}^{\prime }\right)\ge b\left(\mathbf{s},a\right),\\ \hspace{1em}\mathcal{J}\left(\mathbf{s}\right)\in Z^{+},\hspace{0.17em}\hspace{0.17em}\mathbf{s}\in S_m,\end{array}$$

(34)

where $b(\mathbf{s},a)$ are obtained from previously computed stages:

$$\begin{array}{c}b\left(\mathbf{s},a\right)=g\left(\mathbf{s},a\right)+{\displaystyle \sum }_{{\mathbf{s}}^{\prime }\in S_{N\theta }\cup \cdots \cup S_{m-1}}P\left\{{\mathbf{s}}^{\prime }\mid \mathbf{s},a\right\}{\mathcal{J}}^{\ast }\left({\mathbf{s}}^{\prime }\right).\end{array}$$

(35)

The computational method is summarized in Algorithm 1.

Algorithm 1: SSP-MILP algorithm.

input:

$N\theta$ = Minimum Sum energy in states $\mathcal{S}$ ∖ ${\mathcal{S}}_t$ ;

$N{e}^{\max }$ = Maximum Sum energy in states $\mathcal{S}$ ∖ ${\mathcal{S}}_t$ ;

$\mathcal{J}(\mathbf{s})$ = $0$ , $\forall \mathbf{s}\in {\mathcal{S}}_t$ ;

output: $\mathcal{J}(\mathbf{s})$ , $\forall \mathbf{s}\in \mathcal{S}$ ∖ ${\mathcal{S}}_t$ ;

begin

 Find transmission link set using Bron-Kerbosch alg.;

  $m\leftarrow N\theta$ ;

 while   $m\le N{e}^{\max }$   do

  Solve Bin-ball problem to obtain $S_m$ .

  Solve MILP for Subproblem-1 to obtain $\mathcal{J}(\mathbf{s})$ , $\mathbf{s}\in S_m$ .

   $m\leftarrow m+1$ .

 end

end

5.2. Energy Harvesting Networks

The following theorem provides the basis for a linear programming (LP) approach to solve the MDP problem. The proof can be found in [16].

Theorem 3.

Suppose there exists v, for which $\mathbf{v}\ge \mathbb{T}(\mathbf{v})$ ; then $\mathbf{v}\ge {\mathbf{v}}_{\lambda }^{\ast }$ . $\mathbb{T}$ is the nonlinear operator defined as $\mathbb{T}(\mathbf{v})\equiv {\mathrm{s}\mathrm{u}\mathrm{p}}_{\mathrm{a}\mathrm{l}\mathrm{l}\hspace{0.17em}a}\{{\mathbf{r}}_a+\lambda {\mathbf{P}}_a\mathbf{v}\}$ .

From the observation in Theorem 3, the primal LP is constructed as follows.

Primal LP. Consider the following:

$$\begin{array}{c}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}\hspace{1em}{\displaystyle \sum }_{\mathbf{s}\in S}\eta \left(\mathbf{s}\right)v\left(\mathbf{s}\right)\\ \mathrm{s}.\mathrm{t}.\hspace{1em}v\left(\mathbf{s}\right)-{\displaystyle \sum }_{{\mathbf{s}}^{\prime }\in S}\lambda P\left({\mathbf{s}}^{\prime }\mid \mathbf{s},a\right)v\left({\mathbf{s}}^{\prime }\right)\ge r\left(\mathbf{s},a\right)\end{array}$$

(36)

for all $a\in \mathcal{A}\left(\mathbf{s}\right)$ and all $\mathbf{s}\in S$ . $\eta (\mathbf{s})$ are chosen to be positive scalars that satisfy ${\sum }_{\mathbf{s}\in S}‍\eta (\mathbf{s})=1$ .

However, it is more informative to solve this model using its dual LP [16], which has less rows in the constraints matrix.

Dual LP. Consider the following:

$$\begin{array}{c}\mathrm{M}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}\hspace{1em}{\displaystyle \sum }_{\mathbf{s}\in S\hspace{0.17em}}{\displaystyle \sum }_{a\in \mathcal{A}\left(\mathbf{s}\right)}r\left(\mathbf{s},a\right)x\left(\mathbf{s},a\right)\\ \mathrm{s}.\mathrm{t}.\hspace{1em}{\displaystyle \sum }_{a\in \mathcal{A}\left(\mathbf{s}\right)}x\left(\mathbf{s},a\right)-{\displaystyle \sum }_{{\mathbf{s}}^{\prime }\in S\hspace{0.17em}}{\displaystyle \sum }_{a\in \mathcal{A}\left({\mathbf{s}}^{\prime }\right)}\lambda P\left(\mathbf{s}\mid {\mathbf{s}}^{\prime },a\right)x\left({\mathbf{s}}^{\prime },a\right)=\eta \left(\mathbf{s}\right)\end{array}$$

(37)

and $x(\mathbf{s},a)\ge 0$ for $a\in \mathcal{A}(\mathbf{s})$ and $\mathbf{s}\in S$ .

Note that the primal LP has ${\sum }_{\mathbf{s}\in \mathcal{S}}‍\left|\mathcal{A}\left(\mathbf{s}\right)\right|$ rows and $\left|\mathcal{S}\right|$ columns, while the dual LP has $\left|\mathcal{S}\right|$ rows and ${\sum }_{\mathbf{s}\in \mathcal{S}}‍\left|\mathcal{A}\left(\mathbf{s}\right)\right|$ columns. In [44], it is shown that the value function for each state in the primal LP is equal to the dual price (the dual price (a.k.a. shadow price) of a constraint is instantaneous change in the objective function by relaxing the constraint by one unit) corresponding to the constraint associated with the state in the dual LP. Thus, the value function associated with a given initial state of the primal LP is obtained by firstly solving the dual problem and then finding the dual price of the corresponding constraint.

5.3. Computation Complexity of Large-Scale MDP

It has been shown that MDP is solvable in polynomial time by successive approximation methods such as value iteration, without the existence of a known strongly polynomial time algorithm [45]. Notoriously, when dealing with large-scale MDP, the computation of iterations may appear prohibitive. In contrast, computation based on large-scale linear programming is fast solvable with the state-of-the-art LP solver, in polynomial time depending only on the size of A (the constraint matrix in $Ax\le b$ ) [46]. Also, the practical fast ϵ-approximation algorithm is most useful for solving LP with exponentially many constraints to trade off solution accuracy for time [47].

In this paper, we do not target finding the optimal policy. However, we note that standard policy iteration algorithms and their various modifications can be applied offline to obtain the optimal stationary policies, at the cost of iterations over the whole state space. In addition, policies obtained from this network model would require a network centralizer to conduct the execution, which is less appealing than a distributed approximate optimal policy and which is out of the scope of this paper.

6. Numerical Evaluations

We present some numerical results on the optimal lifetime for the EC network and the expected total discounted reward for the EH network, under the 2-hop funnel topology (Nodes A, B, C, and D and the Sink), as in Figure 1. Besides the computational limitation inherent to MDP, the motivation for choosing a small network is as follows. It is observed from our previous work [10] that under duty-cycling, a large network is typically reduced into isolated small networks in the time domain, to reduce interference and collisions. This observation renders the possibility to analyze a small tree topology towards the whole network.

6.1. Energy Constrained Networks

In this section, we evaluate the optimal lifetime of the EH network and impacts from certain parameters. Note that the funnel topology (with the existence of the bottleneck Node C) captures the essence of the energy hole problem.

Figure 3 depicts the optimal lifetime for non-CT and CT networks. We observe the lifetimes are linearly increasing with the battery capacity. We also observe that the performance of CT network is significantly higher than that of the non-CT network. For example, with battery capacity of $10$ units, the lifetime improvement factor of CT network is $1.89$ with $1$ cooperator. We assume that either one or two cooperators are required to reach the Sink. This is motivated by the fact that theory [32] predicts that one cooperator can double the range for certain path loss exponents, but our testbed experiment shows that at least two cooperators are needed to double the range [8]. Considering our VMISO link model in (1), we can show that two helpers ( $\left|\mathcal{C}\left(s\left(l\right)\right)\right|=2$ ) are required for double range instead of one helper, if the path loss exponent increases by a factor $\log (G_3)/\log (G_2)$ . In practice, the number of cooperators required to conduct CT range extension is regulated by the channel status and different CT physical layer techniques; an efficient protocol should select just enough cooperators to avoid energy overuse. The lifetime with larger battery capacities can be obtained from linear regression. From the slope of the extrapolated lines (not shown), the lifetime growth rate with respect to battery capacity is $1.0$ for non-CT, $1.6$ for CT with 2 cooperators, and $2.1$ for CT with $1$ cooperator. Therefore, we conclude that CT gives less of an advantage in terms of lifetime, if more cooperators are required for the same factor of range extension.

Figure 3

Optimal lifetime for non-CT network and CT network, as battery capacity varies. $N=4$ (number of nodes). $N_H=1$ or 2 (required number of cooperators). $e_{\text{th}}=1$ , $e_{tx}=e_{rx}=1$ , $e_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}^{\mathrm{C}\mathrm{T}}=1$ , $e_{\mathrm{c}\mathrm{o}}^{\mathrm{C}\mathrm{T}}=2$ , $q_{\max }=1$ .

Figure 4 shows the computed lifetime for different packet arrival rates (PARs). The PAR is normalized with the packet length, that is, the number of new exogenous packets per node during a packet duration. The PARs where the curves start to become flat represent the PAR thresholds, beyond which the numerical results are accurate. Figure 4 demonstrates that our model is accurate for a very large range of arrival rates (<0.4). For example, if a node generates a packet every 100 seconds and the packet transmission time is 100 ms, then the PAR is 0.001, well below $0.4$ . In other words, the queue model assumption is valid as long as the application generates less than $4$ packets per node per second.

Figure 4

Lifetime for non-CT network and CT network, as the normalized packet arrival rate (PAR) varies. $N=4$ , $N_H=2$ . $e_{\text{th}}=1$ , $e_{tx}=e_{rx}=1$ , $e_{\text{init}}^{\text{CT}}=1$ , $e_{\mathrm{co}}^{\text{CT}}=2$ , $q^{\mathrm{m}\mathrm{a}\mathrm{x}}=1$ .

6.2. Energy Harvesting Networks

In this section, we evaluate the expected total discounted reward and the impacts from certain parameters. The reward function in (27) is used for plotting Figures 5, 6, and 7, while (28) is used for Figure 8.

Figure 5

The effect of discount factor λ. $\alpha =0.1$ , $\gamma =0.01$ , $e_{\max }=4$ , $q_{\max }=1$ .

Figure 6

The effect of energy harvesting rate γ. $\alpha =0.1$ , $\lambda =0.9999$ or $0.99999$ , $e_{\max }=4$ , $q_{\max }=1$ .

Figure 7

The effect of packet arrival rate α. $\lambda =0.9999$ , $\gamma =0.02$ , $e_{\max }=4$ .

Figure 8

The effect of weight factor ${\omega }_1$ . $\alpha =0.1$ , $\lambda =0.99999$ , $\gamma =0.02$ , $e_{\max }=4$ .

Figure 5 shows the expected total discounted reward of non-CT and CT networks, versus different discount factors λ. The decision-making horizon is indicated by $1/(1-\lambda )$ time units. The system parameters are $\alpha =0.1$ , $e_{\mathrm{m}\mathrm{a}\mathrm{x}}=4$ , $q_{\mathrm{m}\mathrm{a}\mathrm{x}}=1$ , $e_{tx}=e_{rx}=e_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}^{\mathrm{C}\mathrm{T}}=1$ , and $e_{\mathrm{c}\mathrm{o}}^{\mathrm{C}\mathrm{T}}=2$ , where $e_{tx}$ and $e_{rx}$ represent transmitting/receiving energy, $e_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}^{\mathrm{C}\mathrm{T}}$ represents the energy consumption of initiating CT, and $e_{\mathrm{c}\mathrm{o}}^{\mathrm{C}\mathrm{T}}$ is the energy consumption of a cooperator. We observe that the rewards increase as λ increases due to an increase in the number of decision epochs, which allows prolonged operation of the network in the optimal states prescribed by the optimal policy. While CT network always outperforms non-CT network ( $55\%$ to $61\%$ improvement), the curve of the CT network also grows with a steeper rate.

Figure 6 depicts the expected total discounted reward of non-CT and CT networks versus energy harvesting rate γ, with $\lambda =0.9999$ or $\lambda =0.99999$ . Other system parameters are the same as used in Figure 5. We observe that, for both non-CT and CT network, the rewards increase linearly when γ is below $0.01$ . Then the curves grow with a slower rate, before the performances of non-CT and CT networks start to converge when the value of γ reaches a threshold value ${\gamma }_{\mathrm{t}\mathrm{h}}$ . For $\lambda =0.9999$ , ${\gamma }_{\text{th}}=0.04$ , and for $\lambda =0.99999$ , ${\gamma }_{\mathrm{t}\mathrm{h}}=0.07$ . This threshold is related to both energy harvesting rate and battery capacity. High harvesting rate, with which packet transmission opportunities are less constrained by energy level, diminishes the benefits of CT. However, in practice, when harvesting rate is not large enough, the CT network still provides a significant gain over the non-CT network, as shown in the figure.

Figure 7 shows the expected total discounted reward of non-CT and CT networks versus exogenous packet arrival rate α, with $q_{\max }=1$ or $q_{\max }=2$ . The system parameters are $\lambda =0.9999$ and $\gamma =0.02$ . It can be seen that when the offered traffic load (α) increases, while not saturating the networks, the rewards show a linear increase. When α is large enough, the rewards are bottlenecked by the queue capacity. We also observe that the reward of non-CT shows slight instability between $\alpha =0.01$ and $0.02$ , between which the rewards slightly decrease, for both $q_{\max }=1$ and $q_{\max }=2$ . We did not explain this observation. However, the CT network does not show this behavior.

Figure 8 gives an example that the calculated benchmark of performance is highly affected by the choice of reward function. The network designer can choose a reward function according to a particular QoS requirement. The x-axis is the preference weighting factor ${\omega }_1$ in the reward function in (28). Note that ${\omega }_2=1-{\omega }_1$ represents the significance of an occurrence of QoS degradation. When a certain application deems a state with a node having zero energy as a severe situation, both non-CT and CT networks will be conservative in initiating a packet transmission, limiting the packets delivered in a given time horizon. On the other hand, the value functions increase as the significance relaxes. Again, CT network still outperforms non-CT network in both magnitude and rate in all cases.

In Figure 9, we evaluate the effect of the CT overhead on the expected total discounted reward of the CT network, with different preference weighting factors ${\omega }_1=\mathrm{0.1,0.5,0.9}$ . The expected total discounted rewards of the non-CT network, which are constant, are also plotted for reference. The CT overhead is modeled as an extended packet transmission time, which is a result of the extra preamble symbols in front of the packet payload to achieve diversity and synchronization for the purpose of range extension [12, 48]. Specifically, the horizontal axis in Figure 9 is the CT overhead, quantified in terms of the fraction of the total packet length, from $0$ to $10\%$ . For ${\omega }_1=0.9$ , at the overhead rate of $0.1$ ( $10\%$ ), the CT network outperforms the non-CT network by $38.2\%$ . For the CT network, comparing with the overhead rate of $0$ , only $2.3\%$ degradation is observed at the overhead rate of $0.1$ . We note that, in a more practical setting with control packets accounted for, we expect to see a bigger degradation value. However, fortunately it is shown that, even with control packets accounted for, CT still has a significant advantage over non-CT in both small networks [9] and large-scale multihop networks [13].

Figure 9

The effect of CT overhead rate. Weight factors ${\omega }_1=\mathrm{0.1,0.5,0.9}$ , $\alpha =0.1$ , $\lambda =0.99999$ , $\gamma =0.02$ , $e_{\max }=4$ .

7. Conclusions

Energy harvesting (EH) and energy constrained (EC) multihop wireless sensor networks are modeled under a unified Markovian decision process framework. The model extends the literature by encompassing, from the network interaction point of view, energy harvesting, routing, MAC, and cooperative transmission. The performance of EC and EH networks is characterized by finite-horizon and infinite-horizon problems, respectively. The EC model captures the traffic imbalance due to the “energy hole,” which is a known bottleneck problem limiting network lifetime. In the EH network, the expected total discounted reward model allows one to formulate different performance metrics of interest. The model of EC networks is solved by mixed integer linear programing, capitalizing on the stochastic shortest path nature of the problem, while the model of EH networks is solved by a dual linear programming approach. Numerical results evaluate the sensitivity of several network parameters on the optimal performance of non-CT and CT networks and show the nontrivial impacts of cooperative transmission.