Heat-Diffusion: Pareto Optimal Dynamic Routing for Time-Varying Wireless Networks

TL;DR

The paper introduces Heat-Diffusion, a dynamic routing policy for wireless networks that optimally balances throughput, routing cost, and delay using only current network conditions, and links it to heat equation mathematics.

Contribution

It develops a novel, throughput and delay optimal routing policy that operates with current network state information and establishes a connection to heat equation analysis.

Findings

HD achieves Pareto optimal tradeoff between routing cost and delay.

HD is throughput optimal and minimizes average delay.

Fluid limit of HD follows the graph heat equation.

Abstract

A dynamic routing policy, referred to as Heat-Diffusion (HD), is developed for multihop uniclass wireless networks subject to random traffic, time-varying topology and inter-channel interference.The policy uses only current condition of queue occupancies and channel states, with requiring no knowledge of traffic and topology.Besides throughput optimality, HD minimizes an average quadratic routing cost defined by endowing each channel with a time-varying cost factor. Further, HD minimizes average network delay in the class of routing policies that base decisions only on current condition of traffic congestion and channel states. Further, in this class of routing policies, HD provides a Pareto optimal tradeoff between average routing cost and average network delay, meaning that no policy can improve either one without detriment to the other. Finally, HD fluid limit follows graph combinatorial heat equation, which can open a new way to study wireless networks using heat calculus, a very active area of pure mathematics.

Figures

Tables

Table 1. Table 1: Algorithmic structure of HD versus V-parameter BP in a uniclass network.

Weighting	$\stackrel{widehat}{f_{ij}}\left(n\right)$	BP	$\min \{{\mu }_{ij}\left(n\right),q_i\left(n\right)\}$
		HD	$\min \left\{{\varphi }_{ij}\left(n\right)q_{ij}\left(n\right){}_{}{}^{+},{\mu }_{ij}\left(n\right)\right\}$
	$w_{ij}\left(n\right)$	BP	${\mu }_{ij}\left(n\right)(q_{ij}\left(n\right)-V{\rho }_{ij}\left(n\right){\mu }_{ij}\left(n\right)){}^{+}$
		HD	$2{\varphi }_{ij}\left(n\right)q_{ij}\left(n\right)\stackrel{widehat}{f_{ij}}\left(n\right)-\stackrel{widehat}{f_{ij}}{\left(n\right)}^2$
Scheduling			$𝝅\left(n\right)=\arg {\max }_{𝝅\in \Pi }{\sum }_{ij\in ℰ}{\pi }_{ij}{w}_{ij}\left(n\right)$
Forwarding

Equations

$$\overline{x}:=\underset{\tau \to \infty }{\mathrm{limsup}}1/\tau {\sum }_{n=0}^{\tau -1}\mathbb{E}\{x(n)\}.$$

$$\mathrm{Minimize}\overline{Q}:={\sum }_{i\in \mathcal{V}}\overline{q_i}.$$

$$\mathrm{Minimize}\;\;\;\,\overline{R}:=\!\sum\nolimits_{ij\in\cal E}\;\overline{\!\rho_{ij}\bigl{(}f_{ij}\bigr{)}{{}^{2}}}\;.$$

$$\mathrm{Minimize}(1-\beta )\overline{Q}+\beta \overline{R}$$

$$w_{ij}(n):={\mu }_{ij}(n)q_{ij}(n){}^{+}.$$

$$\mathbf{\pi }(n)=\arg \underset{\mathbf{\pi }\in \Pi }{\max }{\sum }_{ij\in \mathcal{E}}{\pi }_{ij}{w}_{ij}(n)$$

$$w_{ij}(n):=\mu_{ij}(n)\bigl{(}\,q_{ij}(n)-V\rho_{ij}(n)\,\mu_{ij}(n)\,\bigr{)}{{}^{+}}$$

$$\begin{gathered}\smash{\smash{\widehat{f_{ij}}}}(n):=\min\bigl{\{}\phi_{ij}(n)q_{ij}(n){{}^{+}}\!,\>\mu_{ij}(n)\bigr{\}}\\ \phi_{ij}(n):=(1\!-\!\beta)/\vartheta_{ij}+\beta/\rho_{ij}(n)\end{gathered}$$

$$w_{ij}(n):=2{\varphi }_{ij}(n)q_{ij}(n)\widehat{f_{ij}}(n)-\widehat{f_{ij}}(n{)}^2.$$

$$\mathbf{\pi }(n)=\arg \underset{\mathbf{\pi }\in \Pi }{\max }{\sum }_{ij\in \mathcal{E}}{\pi }_{ij}{w}_{ij}(n)$$

$$f_{ij}(n)=\{\begin{array}{ll}\widehat{f_{ij}}(n)& \text{if }{\pi }_{ij}(n)=1\\ 0& \text{otherwise}\end{array}$$

$$\boldsymbol{q}_{\circ}(n):=\bigl{[}\,q_{1}(n),\ldots,q_{d-1}(n),q_{d+1}(n),\ldots,q_{|\cal V|}(n)\,\bigr{]}.$$

$$\boldsymbol{a}_{\circ}(n):=\bigl{[}\,a_{1}(n),\ldots,a_{d-1}(n),a_{d+1}(n),\ldots,a_{|\cal V|}(n)\,\bigr{]}.$$

$$\boldsymbol{f}(n):=\bigl{[}\,f_{1}(n),\ldots,f_{|\cal E|}(n)\,\bigr{]}$$

$$({\mathbf{B}}_{\circ }\mathbf{f}{)}_i(n)={\sum }_{b\in \mathrm{out}(i)}{f}_{ib}(n)-{\sum }_{a\in \mathrm{in}(i)}{f}_{ai}(n)$$

$${\mathbf{q}}_{\circ }(n+1)={\mathbf{q}}_{\circ }(n)+{\mathbf{a}}_{\circ }(n)-{\mathbf{B}}_{\circ }\mathbf{f}(n).$$

$$\mathbf{\Phi }(n):=\mathrm{diag}(\mathbf{\varphi }(n))$$

$$D(\mathbf{f},{\mathbf{q}}_{\circ },n):=2\mathbf{f}(n{)}^{\top }\mathbf{\Phi }(n){\mathbf{B}}_{\circ }{}^{\top }{\mathbf{q}}_{\circ }(n)-\mathbf{f}(n{)}^{\top }\mathbf{f}(n).$$

$$\overline{D}(\overline{\mathbf{f}},\overline{{\mathbf{q}}_{\circ }}):=2{\overline{\mathbf{f}}}^{\top }\overline{\mathbf{\Phi }}{\mathbf{B}}_{\circ }{}^{\top }\overline{{\mathbf{q}}_{\circ }}-{\overline{\mathbf{f}}}^{\top }\overline{\mathbf{f}}.$$

$$s:=\mathbb{P}\bigl{\{}\boldsymbol{S}(n)\!=\!\boldsymbol{S}\bigr{\}}=\limsup_{\tau\to\infty}{1}/{\tau}\sum\nolimits_{n=0}^{\tau-1}{\mathbb{I}}_{\boldsymbol{S}(n)=\boldsymbol{S}}$$

$$\displaystyle\mathbb{E}\bigl{\{}\boldsymbol{\mu}(n)\}=$$

$$\displaystyle\mathbb{E}\bigl{\{}\boldsymbol{\rho}(n)\}=$$

$$\displaystyle\mathbb{E}\bigl{\{}\boldsymbol{\mu}(n)\}=$$

$$\displaystyle\mathbb{E}\bigl{\{}\boldsymbol{\rho}(n)\}=$$

$$\overline{a_i}={\sum }_{b\in \mathrm{out}(i)}\overline{f_{ib}}-{\sum }_{a\in \mathrm{in}(i)}\overline{f_{ai}},\forall i\in \mathcal{V}\setminus \{d\}$$

$$\boldsymbol{M}_{\!\circ}(n):=\bigl{(}\boldsymbol{B}_{\circ}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\bigr{)}{{}^{-1}}\boldsymbol{B}_{\circ}\boldsymbol{\Phi}(n)\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}.$$

$${\mathbf{B}}_{\circ }{}^{\top }{\mathbf{M}}_{\circ }(n)\mathbf{x}=\mathbf{\Phi }(n){\mathbf{B}}_{\circ }{}^{\top }\mathbf{x}.$$

$$\boldsymbol{x}^{\!\top}\big{(}\boldsymbol{M}_{\!\circ}(n)^{\!\top}\!+\boldsymbol{M}_{\!\circ}(n)\big{)}\,\boldsymbol{y}\leqslant\eta\;\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\!\circ}(n)\,\boldsymbol{y}\,.$$

$$W(n):={\mathbf{q}}_{\circ }(n{)}^{\top }{\mathbf{M}}_{\circ }(n){\mathbf{q}}_{\circ }(n).$$

$$\displaystyle\Delta W(n)=\bigl{(}$$

$$\displaystyle+\bigl{(}\boldsymbol{a}_{\circ}(n)-\boldsymbol{B}_{\circ}\boldsymbol{f}(n)\bigr{)}^{\!\top}\boldsymbol{M}_{\!\circ}(n)\bigl{(}\boldsymbol{a}_{\circ}(n)-\boldsymbol{B}_{\circ}\boldsymbol{f}(n)\bigr{)}.$$

$$\Delta W⩽\eta ({\mathbf{a}}_{\circ }-{\mathbf{B}}_{\circ }\mathbf{f}{)}^{\top }{\mathbf{M}}_{\circ }{\mathbf{q}}_{\circ }+({\mathbf{a}}_{\circ }-{\mathbf{B}}_{\circ }\mathbf{f}{)}^{\top }{\mathbf{M}}_{\circ }({\mathbf{a}}_{\circ }-{\mathbf{B}}_{\circ }\mathbf{f})$$

Full text

1 Introduction

Throughput optimality, which means utilizing the full capacity of a wireless network, is critical to respond to increasing demand for wireless applications. The seminal work in [43] showed that the link queue-differential, channel rate-based Back-Pressure (BP) algorithm is throughput optimal under very general conditions on arrival rates and channel state probabilities. Follow-up works showed that the class of throughput optimal routing policies is indeed large [39, 15, 38, 32]. The challenge is then to develop one that, in addition, is optimal relative to some other important routing objectives.

We propose Heat-Diffusion (HD), a throughput optimal routing policy that operates under the same general conditions and with the same algorithmic structure, complexity and overhead as BP, while also holding the following important qualities: (i) HD minimizes the average quadratic routing cost $\overline{R}$ in the sense of Dirichlet. Endowing each wireless link with a time-varying cost factor, we define average Dirichlet routing cost as the product of the link cost factors and the square of the average link flow rates. Such a generic routing cost may reflect different topology-based penalties, e.g., channel quality, routing distance and power usage, even a cost associated with greedy hyperbolic embedding [41]. (ii) HD minimizes average total queue congestion $\overline{Q}$, which is proportional to average network delay by Little’s Theorem, within the class of routing algorithms that use only current queue occupancies and current channel states, possibly together with the knowledge of arrival/channel probabilities. (iii) In the same class, HD operates on the Pareto boundary of performance region built on the average network delay $\overline{Q}$ and the average quadratic routing cost $\overline{R}$ and can be made to move along this boundary by changing a control parameter $\beta$ that compromises between the two objectives $\overline{Q}$ and $\overline{R}$ (see Fig. 1).

*Related works—*The study of BP schemes has been a very active research area with wide-ranging applications and many recent theoretical results. In packet switches, congestion-based scheduling [15, 16, 27] was extended to admit more general functions of queue lengths with particular interest on $\alpha$-weighted schedulers using $\alpha$-exponent of queue lengths [39]. As another extension in packet switches, [38] introduced Projective Cone Schedulers (PCS) to allow scheduling with non diagonal weight assignments. The work in [32] generalized PCS using a tailored “patch-work” of localized piecewise quadratic Lyapunov functions.

In wireless networks, shadow queues enabled BP to handle multicast sessions with reduced number of actual queues that need to be maintained [12]. Replacing queue-length by packet-age, [22] introduced a delay-based BP policy. To improve BP delay performance, [21] proposed place-holders with Last-In-First-Out (LIFO) forwarding. Adaptive redundancy was used in [1] to reduce light traffic delay in intermittently connected mobile networks. Using graph embedding, [41] combined BP with a greedy routing algorithm in hyperbolic coordinates to obtain a throughput-delay tradeoff.

There have been several reductions of BP to practice in the form of distributed wireless protocols of pragmatically implemented and experimentally evaluated [31, 30, 20]. Some attempts have also been made to adopt the BP framework for handling finite queue buffers [46].

Similar to BP, also HD rests on a centralized scheduling with a computational complexity that can be prohibitive in practice. Fortunately, much progress has recently been made to ease this difficulty by deriving decentralized schedulers with the performance of arbitrarily close to the centralized version as a function of complexity [28, 11, 24].

*Contributions—*We derive HD from combinatorial analogue of classical heat equation on smooth manifolds, which leads to the following key contributions:

(Fluid) Translating “queue occupancy measured in packets” to “heat quantity measured in calories,” the fluid limit of interference HD flow mimics a suitably-weighted non-interference heat flow, in agreement with the second Principle of thermodynamics. In doing so, we introduce a new paradigm that might be called “wireless network thermodynamics,” which builds a rigorous connection between wireless networking and well-studied domains of physics and mathematics.

(Cost) HD reduces the Dirichlet routing cost to its minimum feasible value among all stabilizing routing algorithms. To the best of our knowledge, this is the first time a feasible routing algorithm asserts the strict minimization of a cost function subject to network stability, i.e., bounded queue occupancies and network delay. This is while the drift-plus-penalty approach of [34], as the best-known alternative, can get only close to the minimum of this routing cost at the expense of infinitely large network delay.

(Delay) HD minimizes average queue lengths, and so average network delay, within the class of routing algorithms that act based only on current condition of queue occupancies and channel states, including those with the perfect knowledge of arrival/channel probabilities. This important class contains stationary randomized algorithms [34], original BP policy [43], and most BP derivations [39, 15, 38, 32, 41, 16, 12, 22, 21, 1, 31, 30, 20, 46, 11, 24, 23].

(Pareto) In the class of algorithms defined in (Delay), let the performance region built on average delay and the Dirichlet routing cost be convex. Then HD operates on the Pareto boundary of this region while the optimal tradeoff can solely be controlled by a routing parameter independently of network topology and traffic. This means that no other policy in this class can make a better compromise between these two routing objectives and that any deviation from HD operation leads to the degradation of at least one of them.

(Complexity) Last but not least, HD enjoys the same algorithmic structure, complexity and overhead as BP, giving it the same wide-reaching impact. This also provides an easy way to leverage all advanced improvements of BP to further enhance HD quality. At the same time, it simplifies the way to practice via a smooth software transition from BP to HD.

*Continuation—*The infant idea of HD algorithm first appeared in [3], very different indeed from what is called HD in this paper. The results on minimum network delay are extended to multiclass wireless networks in [6]. By developing the idea of mapping a wireless network onto a nonlinear resistive network, the results on minimum routing cost are extended to multiclass wireless networks in [4]. By extending the principles of classical thermodynamics to routing and resource allocation on wireless networks, the concept of “wireless network thermodynamics” is fully established in [7].

*Organization—*After preliminaries in the next section, HD policy is introduced in Sec. 3 followed by some illustrative examples. Section 4 presents HD key property – a foundation to all HD features. Using Lyapunov theory, throughput optimality is proven in Sec. 5. We show in Sec. 6 that HD minimizes average network delay in a class of routing policies. Physics-oriented model of heat process on directed graphs is proposed in Sec. 7. Using fluid limit theory, Sec. 8 shows that in limit, HD packet flow resembles combinatorial heat flow on its underlying directed graph. Using heat calculus, Sec. 9 shows that HD strictly minimizes the Dirichlet routing cost. HD Pareto optimal performance is discussed in Sec. 10. The paper is concluded in Sec. 11.

*Notation—*We denote vectors by bold lowercase and matrices by bold capital letters. By $\boldsymbol{0}$ we denote the vector of all zeros, by $\boldsymbol{1}$ the vector of all ones, and by $\boldsymbol{I}$ the identity matrix. On arrays: $\min$ and $\max$ are taken entrywise; $\preccurlyeq$ and $\succcurlyeq$ express entrywise comparisons; and $\odot$ denotes the Schur product. For $\boldsymbol{v}$ as a vector, $\boldsymbol{v}^{\top}\!$ denotes its transpose, $\mathrm{diag}(\boldsymbol{v})$ its diagonal matrix expansion, $\|\boldsymbol{v}\|$ its Euclidean norm, and $\boldsymbol{v}{{}^{+}}\!:=\max\{\boldsymbol{0},\,\boldsymbol{v}\}$. For $\cal S$ as a set, $|{\cal S}|$ denotes its cardinality. We use $\mathbb{I}$ as the scalar indicator function, and ${\mathbb{I}}{\mathbb{I}}_{\boldsymbol{v}\succ\boldsymbol{0}}$ as the vector indicator function that its entry $i$ takes the value 1 if $v_{i}\!>\!0$, and 0 otherwise. By $\dot{x}(t)$ we denote the time derivative of $x(t)$. For a variable $x$ related to a directed edge $\ell$ from node $i$ to node $j$, we use notations $x_{\ell}$ and $x_{ij}$ interchangeably. We use $\mathrm{in}(i)$ and $\mathrm{out}(i)$ to denote the sets of nodes neighbor to node $i$ with respectively incoming links to and outgoing links from node $i$.

Note: To keep continuity and enhance readability of the manuscript, proofs are all placed in Appendix.

2 Preliminaries

Consider a uniclass wireless network that operates in slotted time with normalized slots $n\!\in\!\{0,1,2,\cdots\}$. The network is described by a simple, directed connectivity graph with set of nodes $\cal V$ and directed edges $\cal E$. New packets, all with the same destination at node $d\in\cal V$, randomly arrive into different nodes, requiring a multihop routing to reach the destination. Wireless channels may change due to node mobility or surrounding conditions. Assuming the sets $\cal V$ and $\cal E$ change much slower than channel states, we fix them during the time of our interest; then a temporarily unavailable link (due to, e.g., obstacle effect and channel fading) is characterized by zero link capacity. Extended mobility that can lead to permanent change in network topology is not considered here. We assume that channel states remain fixed during a timeslot, while they may change across slots.

In wireless networks, transmission over a channel can happen only if certain constraints are imposed on transmissions over the other channels. An interference model specifies these restrictions on simultaneous transmissions. We consider a family of interference models under which a node cannot transmit to more than one neighbor at the same time. Thus, in a most general case, a node may receive packets from several neighbors while sending packets over one of its outgoing links. Interference constraints used by all well-known network and link layer protocols, including the general K-hop interference models, fall in this family.

Definition 1**.**

Given an interference model, a maximal schedule is such a set of wireless channels that no two channels interfere with each other and no more channel can be added to the set without violating the model constraints.

We describe a maximal schedule with a scheduling vector $\boldsymbol{\pi}\!\in\!\{0,1\}^{|\cal E}$ where $\pi_{ij}=1$ if channel $ij$ is included, and $\pi_{ij}=0$ otherwise.

Definition 2**.**

Given a connectivity graph $({\cal V},{\cal E})$, scheduling set $\Pi$ is the collection of all maximal scheduling vectors.

Definition 3**.**

With $\mathbb{E}$ denoting expectation, the expected time average of a discrete-time stochastic process $x(n)$ is defined as

[TABLE]

Definition 4**.**

A queuing network is stable if queue at each node $i$ and at each slot $n$, denoted as $q_{i}(n)$, has a bounded time average expectation, viz., $\overline{q_{i}}<\infty$.

Definition 5**.**

Given a wireless network, an arrival vector $\boldsymbol{a}(n)$ is stabilizable if there exists a routing policy that can make the network stable under $\boldsymbol{a}(n)$.

For a link $ij$, its capacity $\mu_{ij}(n)$, which is frequently called transmission rate in literature, counts the maximum number of packets the link can transmit at slot $n$. The link actual-transmission $f_{ij}(n)$, on the other hand, counts the number of packets genuinely sent over the link at slot $n$. Each link is also endowed with a cost factor $\rho_{ij}(n)\geqslant 1$ that represents the cost of transmitting one packet over the link at slot $n$; for example, $\rho_{ij}=\mathrm{ETX}_{ij}$, with $\mathrm{ETX}$ as defined in [13], or a cost associated with greedy embedding [41].

2.1 Problem Statement

For a constrained uniclass network described above, we propose HD routing policy that solves the three stochastic optimization problems as follows. It is important to note that these optimization problems must be solved at network layer alone, which makes it totally different from cross-layer optimization [29, 35, 17, 42] that aims to control congestion by controlling arrival rates into network layer. With no control on arrivals, the basic assumption here is that arrival rates lie within network capacity region, making the routing system stabilizable. Obviously, nothing prevents one to either install a flow controller on top of HD or develop an HD-based Network Utility Maximization (NUM) protocol.

(Delay) Average network delay minimization:

[TABLE]

Solving this problem for a general case requires the Markov structure of network topology process, plus arrival and channel state probabilities. Then in theory, the solution is obtained through dynamic programming for each possible topology along with solving a Markov decision problem. By even having all this required information, the number of queue backlogs and channel states increase exponentially with the size of network, which makes dynamic programming and Markov decision theory prohibitive in practice. In fact, even for the case of a single channel, it is hard to implement the resulting stochastic algorithms [8]. While having a practical solution for a general case seems dubious, we show in Th. 3 that HD policy solves this problem within an important class of routing algorithms, without requiring any of the above-mentioned information or dealing with any dynamic programing or Markov decision process.

(Cost) Average quadratic routing cost minimization:

[TABLE]

The loss function $\overline{R}$, by concept, spreads out traffic with a weighted bias towards lower penalty links that reminds the optimal diffusion processes in physics, such as heat flow and electrical current [33]. It is shown in [18, 34] that a stationary randomized algorithm can solve this problem. While such an algorithm exists in theory, it is intractable in practice as it requires a full knowledge of traffic and channel state probabilities. Further, assuming all of the probabilities could be accurately estimated, the network controller still needs to solve a dynamic programming for each topology state, where the number of states grows exponentially with the number of channels. Nonetheless, we show in Th. 8 that HD policy solves this problem without requiring any knowledge of traffic and channel state probabilities or dealing with any dynamic programming.

(Pareto) Pareto optimal performance:

[TABLE]

where $\beta\in[0,1]$ is a control parameter to determine relative importance between average delay and average routing cost, which naturally plays the role of Lagrange multiplier too. To our knowledge, this is the first time in literature that such a multi-objective optimization problem is addressed in the level of solely network layer. While even the related single-objective optimization problems are not easy to manage, we show in Th. 9 that within the same class of routing policies mentioned in (Cost), HD policy solves problem (4) subject to convex Pareto boundary on the feasible $(\overline{Q},\overline{R})$ region, with requiring no knowledge of traffic and topology.

2.2 Back-Pressure (BP) Policy

At each slot $n$, the original BP [43] observes queue backlogs $q_{i}(n)$ at network layer and estimates channel capacities $\mu_{ij}(n)$ to make a routing decision as follows.

1. 1)

BP weighting: For every link $ij$ find link queue-differential $q_{ij}(n):=q_{i}(n)-q_{j}(n)$ and weight the link with

[TABLE] 2. 2)

BP scheduling: Find a scheduling vector such that

[TABLE]

where ties are broken arbitrarily. 3. 3)

BP forwarding: Over each activated link with $w_{ij}(n)>0$ transmit packets at full capacity $\mu_{ij}(n)$. If there is no enough packets at node $i$, transmit null packets.

2.3 V-Parameter BP Policy

Thus far, the drift-plus-penalty approach [18, 34], which we refer to as V-parameter BP hereafter, has been the only feasible approach to decreasing (not minimizing) a generic routing penalty at network layer. We take the V-parameter BP as a yardstick as to how HD performs. To incorporate average routing cost $\overline{R}$ into the original BP, the V-parameter BP adds a usage cost to each link queue-differential via replacing the link weight of BP by

[TABLE]

where $V\in[0,\infty)$ trades queue occupancy for routing penalty, while $V\!=0$ recovers the original BP.

The V-parameter BP yields a Dirichlet routing cost within $O(1/V)$ of its minimum feasible value to the detriment of growing average delay of $O(V)$ relative to that of original BP [34]. Thus, the policy is not able to achieve minimum routing cost subject to finite network delay, i.e., delay grows to infinity as routing cost is pushed towards its minimum. Another issue is that the resulting tradeoff depends on both $V$ and the network, with two negative consequences: (i) The same value of $V$ leads to different levels of tradeoff in different networks, and (ii) The level of tradeoff in the same network varies by topology and arrival rates, making it difficult to find a proper $V$ in practice.

3 Heat-Diffusion (HD) Policy

To provide a convenient way of unifying our proposed scheme with the large body of previous works on BP, we design HD with the same algorithmic structure, complexity and overhead, in both computation and implementation, as BP.

3.1 HD Algorithm

At each slot $n$, HD policy observes link queue-differentials $q_{ij}(n):=q_{i}(n)-q_{j}(n)$ at network layer and estimates channel capacities $\mu_{ij}(n)$ and channel cost factors $\rho_{ij}(n)$ to make a routing decision as follows.

1. 1)

HD weighting: For every link $ij$ first calculate the number of packets it would transmit if it were activated as

[TABLE]

where $\vartheta_{ij}=1$ if node $j$ is the final destination, i.e., $j=d$, and $\vartheta_{ij}=2$ otherwise. The Lagrange control parameter $\beta$ is as defined in (4) to make a tradeoff between queue occupancy and routing penalty, and the hat notation denotes a predicted value which would not necessarily be realized. Then determine the link weight as

[TABLE] 2. 2)

HD scheduling: Find a scheduling vector, in the same way as BP, using the max-weight scheduling, such that

[TABLE]

where ties are broken arbitrarily. 3. 3)

HD forwarding: Over each activated link transmit $\smash{\widehat{f_{ij}}}(n)$ number of packets, viz.,

[TABLE]

where $f_{ij}(n)$ represents the number of packets genuinely sent over link $ij$ at slot $n$.

It is critical to discriminate among actual link transmissions $f_{ij}(n)$, link transmission predictions $\smash{\widehat{f_{ij}}}(n)$ and link capacities $\mu_{ij}(n)$. Also notice that $\widehat{f_{ij}}(n)$ in (6) could be non-integer. In practice, the final number of packets to be transmitted over links can be rounded to the nearest integer to $\smash{\widehat{f_{ij}}}(n)$ with no important influence on the performance. To be more precise, however, every node may algebraically add the packet residuals sent on each of its ongoing links so as to make a compensation as soon as the sum hits either $1$ or $-1$.

Table 1 compares HD and V-parameter BP algorithms, which emphasizes the same algorithmic structure, computational complexity and overhead signaling.

Remark 1*.*

(i) Since $\rho_{ij}(n)\geqslant 1$ by assumption, we get $0<\phi_{ij}(n)\leqslant 1$ for all $\beta\!\in\![0,1]$. (ii) If $q_{ij}(n)\leqslant 0$, we get $\smash{\widehat{f_{ij}}}(n)=0$ due to (6) and $w_{ij}(n)=0$ from (7); in this case, even if the link were scheduled by (8), still no packet would be transmitted over it. (iii) If $q_{ij}(n)>0$, we get $q_{ij}(n){{}^{+}}\!=q_{ij}(n)$ and since $\smash{\widehat{f_{ij}}}(n)\leqslant\phi_{ij}(n)q_{ij}(n)$ due to (6), the link weight (7) still remains positive. (iv) In light of $q_{ij}(n){{}^{+}}\!\leqslant q_{i}(n)$ and $\phi_{ij}(n)\leqslant 1$, the value of $\smash{\widehat{f_{ij}}}(n)$ never exceeds the number of packets in the transmitting node $i$.

Remark 2*.*

For $\beta=0$, HD policy reduces to the initial adiabatic-based HD policy proposed in [3], where the packet forwarding follows a thermally adiabatic, and so insulated, heat process on each link.

Remark 3*.*

In a special case that all links are of the same capacity, i.e., $\mu_{ij}(n)=\mu(n)$, and all link queue-differentials remain less than it, i.e., $q_{ij}(n)<\mu(n)$, HD policy with $\beta=0$ and $\alpha$-weighted policy of [39] with $\alpha=2$ turn to be equivalent. Packet switches are well suited to this special case. It was suggested in [39] that a smaller $\alpha$ may lead to a lower network delay, with a non-proven conjecture that heavy traffic delay is minimized when $\alpha\to 0$. A discussion of this was given in [23] along with some counterexamples. Even if the conjecture were true, note that for a multihop routing problem, the requirement of $q_{ij}(n)\!<\mu_{ij}(n)$ would imply the network not to be in a heavy traffic condition.

3.2 Highlights of HD Design

H1: While BP is derived by link capacity $\mu_{ij}(n)$, HD emphasizes on actual number of transmittable packets $\smash{\widehat{f_{ij}}}(n)$, though it also implicitly takes the link capacity into account through (6). Thus, HD allocates resources based only on genuinely transmittable packets, without counting on null packets as being practiced in BP schemes.

H2: The link weight (7), which itself directly controls the scheduling optimization problem, is taken quadratic in the link queue-differential $q_{ij}(n)$, where for $\phi_{ij}(n)q_{ij}(n)\leqslant\mu_{ij}(n)$ is simplified to $w_{ij}(n)=\phi_{ij}(n){{}^{2}}q_{ij}(n){{}^{2}}$. This contrasts with BP weighting $w_{ij}(n)=\mu_{ij}(n)q_{ij}(n)$ which is linear in $q_{ij}(n)$. The quadratic weight is central to HD key property (Th. 1) which is fundamental to other HD qualities.

H3: Varying the penalty factor $\beta$ makes a universal tradeoff in performance that depends neither on network nor on arrivals with the following significant results:

• •

HD is throughput optimal for all $\beta\in[0,1]$ (Th. 2).

• •

At $\beta=0$, the average total queue $\overline{Q}$, and so average network delay, decrease to their minimum feasible values within the class of routing policies that rely only on present queue backlogs and current channel states (Th. 3).

• •

Raising $\beta$ adds to average delay in return for a lower routing cost, where the exclusive merit of HD is to provide the best tradeoff between these two criteria (Th. 9).

• •

At $\beta=1$, the average routing cost $\overline{R}$ reaches its minimum (Th. 8) through an optimal tradeoff with average network delay. Note that in V-parameter BP, network delay grows to infinity as routing cost is pushed towards its minimum.

H4: Unlike BP that forwards the highest possible number of packets over activated links, HD controls packet forwarding by limiting it to $\phi_{ij}(n)q_{ij}(n)$ with $\phi_{ij}$ changing between 0 and 1 as a function of $\beta$, $\vartheta_{ij}$ and $\rho_{ij}$. This reduces queue oscillations by decreasing unnecessary packet forwarding across links, which itself reduces total power consumption and routing penalty. Thus, it is not surprising to see that $\phi_{ij}$ is decreasing, and so as to have a higher impact, by increasing $\beta$ that means more emphasis on routing penalty. Forwarding a portion of link queue-differentials rather than filling up link capacities also complies with resembling heat flow on the underlying directed graph (Th. 5) that in effect minimizes time average routing cost in light of Dirichlet’s principle (Th. 8).

Figure 1 provides a graphical comparison between operation of HD for ${\beta\in[0,1]}$ and V-parameter BP for $V\in[0,\infty)$. The performance region is restricted to the set of all $\overline{Q}$ achievable by the class of all routing policies that act based only on present queue backlogs and current channel states, and is assumed to have a convex Pareto boundary.

3.3 Illustrative Examples

In order to focus merely on the policy itself, we take everything deterministic in our examples here, resting assure that the results purely show the policy performance not contaminated by stochastic effects. We however know that all HD properties are analytically proven for stochastic arrivals and random topologies under very general conditions.

Two-queue downlink: Consider a base station that transmits data to two downlink users, where at most one link can be activated at each timeslot. Let link 1 be of constant capacity $\mu_{1}=3$ (packets/slot) and link 2 of time-varying capacity $\mu_{2}\geqslant 2$. Assume one packet to arrive for each user at every timeslot. It is then easy to verify that for $\mu_{2}<1.5$, the given arrival goes beyond the network capacity region.

For $q_{1}(0)=q_{2}(0)=0$, Fig. 2 compares the performance of HD with $\beta=0$ and original BP. The left side panel depicts timeslot evolution of $q_{1}(n)+q_{2}(n)$ for $\mu_{2}=18$. The right side panel shows the steady-state average of total queue length as a function of $\mu_{2}$. For $2\leqslant\mu_{2}\leqslant 5$, both HD and BP perform the same. For $\mu_{2}\geqslant 5$, however, average total queue length increases linearly in $\mu_{2}$ under BP, while HD holds the optimal performance for all admissible link capacities. This exemplifies H1 in the previous subsection, i.e., the efficiency of link scheduling based on actual transmittable packets rather than link capacities.

Lossy link network: Consider the 4-node network of Fig. 3 with lossy links and subject to 1-hop interference model, i.e., two links with a common node cannot be activated at the same time. The links are labeled with both ETX and capacity, where ETX is a quality metric defined as the expected number of data transmissions required to send a packet without error over a link [13]. Assume that at every timeslot a single packet arrives at node 1 destined for node $d$. Following [31], let us take $\rho_{ij}=\nobreak\mathrm{ETX}_{ij}$.

For zero initial conditions, Fig. 3 compares the performance of HD with V-parameter BP. While HD easily stabilizes total queued packets at 1 for any $\beta>\nobreak 0$, trying with different values of $V$ indicates the weakness of V-parameter BP in aptly supporting the arrival. This simplistically shows one of the impacts of entering link cost factor $\rho_{ij}$ as a multiplicand in the HD weighting formula (7) rather than an addend in the V-parameter BP weighting formula (5).

Power minimization: Consider the sensor network of Fig. 4 subject to 1-hop interference model. Suppose that each link $ij$ has a noise intensity $N_{ij}\in[1,5]$ which is randomly assigned at first and keeps fixed during the simulation. For each link, we adopt Shannon capacity $\mu_{ij}=\Omega_{ij}\log_{2}(1+P_{ij}/N_{ij})$ with $P_{ij}$ as power transmission and $\Omega_{ij}$ as bandwidth. At every timeslot, two packets arrive at nodes 1, 2, 3 and 4, destined for node $d$. The aim is to minimize total $\smash{\rho_{ij}(f_{ij}){{}^{2}}}$ with $\rho_{ij}:=P_{ij}/\mu_{ij}$, which implicitly minimizes total power consumption in the network. For simplicity, let us fix $P_{ij}=15$ and $\Omega_{ij}=5$ for all links so that the capacity on each link is decided only by its noise intensity.

Figure 4 displays timeslot evolution of total queue length for HD with $\beta=0$ and for the original BP ($V\!=0$). Average queue congestion is minimized at about 50 packets under HD, compared with over 100 packets under original BP. Further, little steady-state oscillations in total queue congestion under HD contrary to its large variations under BP verifies H4 in the previous subsection.

In minimizing average routing cost, Fig. 5 displays timeslot evolution of total power consumption for HD with $\beta=1$ and for V-parameter BP with $V=10$. Note that while the total power consumption and the average routing cost are not identical, they are highly correlated with each other. Smaller steady-state oscillations in total power under HD endorses both H1 and H4 in the previous subsection, showing the defect of link capacity-driven scheduling and maximum packet forwarding by BP.

Figure 6 displays the tradeoff between queue congestion and power usage in HD as a function of $\beta$ and in V-parameter BP as a function of $V$. The results verify H3 in the previous subsection and concur with the graphical illustration of HD Pareto optimal performance depicted by Fig. 1. They also match the timeslot evolution results displayed in Fig. 4 for total queue length at $\beta=0$ and $V=0$, and in Fig. 5 for total power consumption at $\beta=1$ and $V=10$. Note the rapid growth of queue lengths in V-parameter BP when average power usage is pushed downwards, indicating the fact that the V-parameter BP cannot reach the minimum routing cost subject to network stability, i.e., bounded queue lengths.

4 Key Property of HD Policy

Consider a general uniclass queuing network with a single destination node $d$. As before, let $q_{i}(n)$ be the number of existing packets at node $i$ at slot $n$. State variables of the system can then be represented by the following vector:

[TABLE]

Note that $q_{d}(n)\equiv 0$ is discarded from state variables.

Notation 1*.*

We use subscript $\circ$ to denote reduced vectors or matrices obtained by discarding the entries corresponding to the destination node $d$.

Let a stochastic process $a_{i}(n)$ represent the number of exogenous packets arriving into node $i$ at slot $n$. Discard $a_{d}(n)\equiv 0$ and compose the vector of node arrivals as

[TABLE]

Also compose the vector of link actual transmissions as

[TABLE]

where, as before, $f_{ij}(n)$ represents the number of packets actually sent over link $ij$ at slot $n$.

Given a directed graph $(\cal V,\cal E)$, let $\boldsymbol{B}$ denote the node-edge incidence matrix, in which $B_{i\ell}$ is 1 if node $i$ is the tail of directed edge $\ell$, $-1$ if $i$ is the head, and 0 otherwise.111 In combinatorial geometry, one can view graph as a 1-complex, where $\boldsymbol{B}$ is its 1-incidence matrix that describes the correlation between all oriented 1-cells (edges) and 0-cells (nodes) in the complex. Then $\boldsymbol{B}_{\circ}$ denotes a reduction of $\boldsymbol{B}$ that discards the row related to the destination node $d$, which is referred to as reduced incidence matrix. One can verify that $\boldsymbol{B}_{\circ}\boldsymbol{f}(n)$ is a node vector, in which the entry corresponding to node $i$ reads the net outflow as

[TABLE]

Using the above notation, the $\boldsymbol{f}$-controlled, stochastic state dynamics of a uniclass queuing network is captured by

[TABLE]

Note that the link capacities $\mu_{ij}(n)$ vary by channel states, while the link actual transmissions $\boldsymbol{f}_{ij}(n)$ are assigned by a routing policy subject to $0\leqslant f_{ij}(n)\leqslant\min\{q_{i}(n),\,\mu_{ij}(n)\}$. This difference explains why despite traditional notation in literature, there is no need for $(\cdot){{}^{+}}$ operation in the queue equation (10).

In the wake of (10), the next theorem formalizes the HD main characteristic, which is central to the proof of Th. 2 on HD throughput optimality, Th. 3 on HD average network delay minimization, Th. 5 on connection between HD fluid limit and combinatorial heat equation, and Th. 8 on HD average quadratic routing cost minimization. Before proceeding to the theorem, let us define the link weight matrix as

[TABLE]

where $\boldsymbol{\phi}(n)$ represents the vector composed of $\phi_{ij}(n)$ as defined in (6).

Theorem 1** (HD Key Property).**

Consider a uniclass wireless network constraint by capacity, directionality and interference. At every timeslot $n$ and for all $\beta\in[0,1]$, HD policy maximizes the $\boldsymbol{f}$-controlled functional

[TABLE]

Consider the long-term average of functional $D(\boldsymbol{f},\boldsymbol{q}_{\circ},n)$ defined as

[TABLE]

Next assumption is being used in the analytical proofs of HD properties, stating that the greedy maximization of $D(\boldsymbol{f},\boldsymbol{q}_{\circ},n)$ at each timeslot leads to its maximum long-term average. The assumption implies that one can apply the Bellman’s principle of optimality, and so dynamic programming, to maximize $\overline{D}$. It also implicitly means no overlapping among slot-based substructures of $\overline{D}$ maximization problem.

Assumption 1**.**

Consider a uniclass wireless network constraint by capacity, directionality and interference. Given a combination of network topology and traffic rates, timeslot maximization of $D(\boldsymbol{f},\boldsymbol{q}_{\circ},n)$ is an optimal substructure for global maximization of $\overline{D}(\,\overline{\boldsymbol{f}},\,\overline{\boldsymbol{q}_{\circ}\!}\,)$.

In practice, almost every wireless mesh network meets this assumption. As an example that fails the requirement though, consider the case where exogenous packets arrive only to one node, say $a$, which is connected directly to the final destination. Assume that all links are bidirectional with unit cost factors and infinite capacities, and so link interference is the only network constraint. Obviously, depleting the whole queue into the destination maximizes $D$ to $q_{a}(n)^{2}$ at each timeslot. To maximize $\overline{D}$, however, a portion of traffic must be forwarded through other paths that connect node $a$ to the destination.

5 HD Throughput Optimality

Let the stochastic process \boldsymbol{S}(n)=\bigl{(}S_{1}(n),\cdots,S_{|{\cal E}|}(n)\bigr{)} represent channel states at slot $n$, describing all uncontrollable factors that affect wireless link capacities and cost factors. We assume that $\boldsymbol{S}(n)$ evolves according to an ergodic stationary process and takes values in a finite set $\cal S$. Thus, by Birkhoff’s ergodic theorem, each state $\boldsymbol{S}\in\cal{S}$ has a probability of

[TABLE]

where $\sum_{\boldsymbol{S}\in\cal S}s=1$. Then the expected link capacities and cost factors are obtained as

[TABLE]

where $\boldsymbol{\mu}(n)$ and $\boldsymbol{\rho}(n)$ represent the vectors composed of link capacities $\mu_{ij}(n)$ and link cost factors $\rho_{ij}(n)$, respectively.

Note that the existence of probability distribution (14) or expected values (15) and (16) by no means imply that they are known to a routing policy. Specifically, HD performs without knowing any of these information. Nonetheless, the ergodicity of $\boldsymbol{S}(n)$ along with the law of large numbers imply

[TABLE]

meaning that the expectations converge to the long-term averages. Thus, a routing policy could estimate $\mathbb{E}\{\boldsymbol{\mu}(n)\}$ and $\mathbb{E}\{\boldsymbol{\rho}(n)\}$ by observing timeslot variables $\boldsymbol{\mu}(n)$ and $\boldsymbol{\rho}(n)$ for a long enough period of time, at least in theory. This justifies the existence of stationary randomized policies that base their routing decisions only on arrival statistics and channel state probabilities, but fully independent of queue occupancies.

5.1 Characteristic of Network Capacity Region

Consider a uniclass wireless network that is described by a connectivity graph $(\cal{V},\cal{E})$, a destination node $d$, and an ergodic stationary channel state process $\boldsymbol{S}(n)$.

Definition 6**.**

Given a routing policy, its stability region is the set of all arrival vectors that it can stably support, i.e., make the network stable under those arrivals.

Definition 7**.**

Given a network layer, its capacity region is the union of stability regions achieved by all routing policies, including those which are possibly unfeasible.

It can be shown that for any network, its capacity region is convex and compact and so is closed and bounded [18].

Definition 8**.**

A routing policy is throughput optimal if it can stabilize the entire network capacity region, i.e., secure queue stability under all stabilizable arrival vectors.

An arrival vector $\boldsymbol{a}_{\circ}(n)$ is in the network capacity region, i.e., stabilizable, if and only if there exists a set of link actual transmissions $\boldsymbol{f}(n)$ that satisfy

[TABLE]

constrained by link capacities and interference. Under an ergodic channel state process, this basically reads the long-term average flow conservation at the nodes. In a matrix form, (17) can equivalently be shown by $\overline{\boldsymbol{a}_{\circ}\!}\,=\boldsymbol{B}_{\circ}\overline{\boldsymbol{f}}$.

Remark 4*.*

Link actual transmissions $\boldsymbol{f}(n)$ are not fixed, but depend on routing policy. Further, there could potentially exist infinite number of routing policies that meet (17) for any stabilizable $\boldsymbol{a}_{\circ}(n)$. Among them are the ones that use the simple probability concept of distributing packets randomly so that the desired time averages (17) can be achieved. These stationary randomized policies are prohibitive in practice as they typically require perfect knowledge of arrival statistics and channel state probabilities along with an expensive computation. Nonetheless, the fact that these queue-independent policies exist plays a crucial role in the analytical proof of HD properties in this and next section.

5.2 HD Throughput Optimality for all $\beta$

To prove network stability under HD policy, as well as some other HD properties in next sections, we are compelled to choose unorthodox Lyapunov candidates based on the following nonsymmetric system matrix:

[TABLE]

Handling Lyapunov arguments turns to be a lot more challenging, since the easy way of working with symmetric positive definite matrices ceases to exist here. Nonetheless, the specific structure of $\boldsymbol{M}_{\!\circ}(n)$ makes the following lemmas possible.

Lemma 1**.**

Given a connected uniclass wireless network, $\boldsymbol{M}_{\!\circ}(n)$ is pseudo positive definite in the sense that all of its eigenvalues are positive and $\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\!\circ}(n)\,\boldsymbol{x}\geqslant 0$ for any vector $\boldsymbol{x}\in{\mathbb{R}}^{|{\cal V}|-1}$, with equality if and only if $\boldsymbol{x}=\boldsymbol{0}$.

Lemma 2**.**

Given a connected uniclass wireless network, for any vector $\boldsymbol{x}\in{\mathbb{R}}^{|{\cal V}|-1}$, the following identity holds:

[TABLE]

Lemma 3**.**

Given a connected uniclass wireless network, there exists such a scalar $1\leqslant\eta\leqslant 3$ that for any vectors $\boldsymbol{x},\boldsymbol{y}\in{\mathbb{R}}^{|{\cal V}|-1}$, the following inequality holds:

[TABLE]

To analyze the HD throughput optimality, consider the Lyapunov candidate

[TABLE]

Though $W(n)$ is indeed an energy function in light of Lem. 1, due to the nonsymmetric weighting matrix $\boldsymbol{M}_{\!\circ}(n)$, it has no trivial interpretation of a specific energy in the system. Nonetheless, it clearly penalizes high queue differentials across links, compelling a more even distribution of packets over the network. It also incites transmission over the links of lower cost factors, leading to a less expensive routing decision. Note that either at $\beta\!=\!0$ or for the case that all links are of the same cost factor, $\boldsymbol{\Phi}(n)$ is simplified to a scaled identity matrix that leads to $\boldsymbol{M}_{\!\circ}(n)=\boldsymbol{\Phi}(n)$, which in turn reduces $W(n)$ to the sum of squares of queue lengths – a familiar Lyapunov function in most of previous results in literature.

Let $\Delta W(n):=W(n+1)-W(n)$ be the Lyapunov drift. Substituting for $\boldsymbol{q}_{\circ}(n+1)$ from (10) leads to

[TABLE]

Let us drop timeslot variable $(n)$ for ease of notation. Applying Lem. 3 to the first line of the above drift equation yields

[TABLE]

with $1\leqslant\eta\leqslant 3$. Let us replace $\boldsymbol{f}^{\top}\!\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{M}_{\!\circ}\boldsymbol{q}_{\circ}$ by $\boldsymbol{f}^{\top}\!\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{q}_{\circ}$ in light of Lem. 2, add and subtract the term $\frac{1}{2}\,\eta\,\boldsymbol{f}^{\top}\!\boldsymbol{f}$, and use the $D(\boldsymbol{f},\boldsymbol{q}_{\circ},n)$ expression in (12) to obtain

[TABLE]

Taking conditional expectation from the latter given the current queue backlogs $\boldsymbol{q}_{\circ}(n)$ and knowing that the term $\eta\,\boldsymbol{f}^{\!\top}\!\boldsymbol{f}$ has a zero lower bound lead to

[TABLE]

where the conditional expectation is with respect to the randomness of arrivals, channel states and routing decision – in case of a randomized routing algorithm.

Observe that \boldsymbol{M}_{\!\circ}(n)=\bigl{(}\boldsymbol{B}_{\circ}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\bigr{)}{{}^{-1}}\boldsymbol{B}_{\circ}\boldsymbol{\Phi}(n)\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}} is a function only of control parameter $\beta$ and link cost factors $\rho_{ij}(n)$. Since arrivals are independent of both $\beta$ and $\rho_{ij}$, we get

[TABLE]

At the same time, both $\beta$ and $\rho_{ij}$ are independent of $\boldsymbol{q}_{\circ}$, so is $\boldsymbol{M}_{\!\circ}$, which means $\mathbb{E}\{\boldsymbol{M}_{\!\circ}|\boldsymbol{q}_{\circ}\}=\mathbb{E}\{\boldsymbol{M}_{\!\circ}\}$. On the other hand, since the network layer routing controller has no impact on arrivals, $\boldsymbol{a}_{\circ}(n)$ turns to be an independent system variable that is not influenced by anything, which implies $\mathbb{E}\{\boldsymbol{a}_{\circ}{\!\!\!{}^{\top}}|\boldsymbol{q}_{\circ}\}=\mathbb{E}\{\boldsymbol{a}_{\circ}{\!\!\!{}^{\top}}\}$. Putting these results together yields

[TABLE]

Given the current queue backlogs $\boldsymbol{q}_{\circ}(n)$, let $\boldsymbol{f}^{\star}(n)$ be the link actual transmissions provided by HD policy. As compared to any alternative transmission decision $\boldsymbol{f}(n)$, Th. 1 secures $D(\boldsymbol{f}^{\star}\!,\boldsymbol{q}_{\circ},n)\geqslant D(\boldsymbol{f},\boldsymbol{q}_{\circ},n)$ for all $\beta$ and at each slot $n$. Considering this with the equality (19) of Lem. 2 implies

[TABLE]

Taking conditional expectation given current queues yields

[TABLE]

As one alternative transmission decision $\boldsymbol{f}(n)$ to be compared with the $\boldsymbol{f}^{\star}(n)$ provided by HD policy, consider the case where $\boldsymbol{f}(n)$ is produced by a routing algorithm which makes independent, stationary and randomized transmission decisions at each slot $n$ based only on arrivals and link capacities and so independent of both queue backlogs and link cost factors [18]. Let us fix $\boldsymbol{f}(n)$ for such an algorithm and refer to it as $\boldsymbol{f^{\prime}}(n)$. Using equality $\mathbb{E}\{\boldsymbol{M}_{\!\circ}|\boldsymbol{q}_{\circ}\}=\mathbb{E}\{\boldsymbol{M}_{\!\circ}\}$ and considering that $\boldsymbol{f^{\prime}}(n)$ is independent from $\boldsymbol{q}_{\circ}(n)$ and $\boldsymbol{M}_{\!\circ}(n)$, we obtain

[TABLE]

Exploiting this and (22) in (21) leads to the following Lyapunov drift inequality which is evaluated under HD policy given current queue backlogs at slot $n$:

[TABLE]

Investigating $\Gamma(n)$, note that (i) all arrivals are of finite mean and variance, (ii) each link actual-transmission is at most equal to the link capacity which is finite, and so both $\boldsymbol{f}^{\star}(n)$ and $\boldsymbol{f^{\prime}}(n)$ have finite upper bounds, and (iii) $\boldsymbol{M}_{\!\circ}(n)$ is a pseudo positive definite matrix in the sense of Lem. 1 with finite entries (recall $\phi_{ij}(n)\leqslant 1$). Thus, the expected value of $\Gamma(n)$ is finite at each slot $n$, and so there exists a finite positive scalar $\Gamma_{\!\max}$ such that $\mathbb{E}\{\Gamma(n)\,|\boldsymbol{q}_{\circ}(n)\}\leqslant\Gamma_{\!\max}$. Utilizing this in the Lyapunov drift inequality yields

[TABLE]

In the wake of (23), the next theorem is proven by showing that $\mathbb{E}\{\Delta W\,|\boldsymbol{q}_{\circ}\}$ is always negative for all $\beta\in[0,1]$. (Refer to the Appendix for the end of the proof.)

Theorem 2** (HD Throughput Optimality).**

Over any uniclass wireless network, HD policy with any $\beta\in[0,1]$ is throughput optimal, meaning that it guarantees network stability under all stabilizable arrival vectors.

6 HD Minimum Delay at $\beta=0$

Pareto optimal performance of HD policy stands on two pillars: minimization of the average queue congestion $\overline{Q}$ with $\beta=0$, and minimization of the average routing cost $\overline{R}$ with $\beta=1$. This section settles the first pillar based on a timeslot analysis. The result of this section is analytically proven under the general K-hop interference model, where two wireless links can be activated at the same time if they are at least K$+1$ hops away from one another. For example, in the 1-hop interference model, links with the exclusive nodes may be scheduled at the same time. Let us start with two lemmas (proof in the appendix) that help us analyze the final delay minimization in Th. 3.

Lemma 4**.**

At $\beta=0$ and under the K-hop interference model, timeslot maximization of the functional $D(\boldsymbol{f},\boldsymbol{q}_{\circ},n)$ in (12) is equivalent to timeslot maximization of

[TABLE]

It is critical to understand that Lem. 4 does not claim about the same maximum values for functionals $D$ and $G$, which is obviously not true, but about the same maximizing control action $\boldsymbol{f}(n)$ at each slot $n$. Another point is that while at each timeslot, HD maximizes $D$ for all $\beta\in[0,1]$, it maximizes $G$ for only $\beta=0$.

Lemma 5**.**

Consider a uniclass wireless network under an arrival rate $\overline{\boldsymbol{a}_{\circ}\!}\,$ that is stabilized by a routing policy, resulting in average queue occupancies $\overline{\boldsymbol{q}_{\circ}\!}\,$ and average link actual transmissions $\overline{\boldsymbol{f}}\,$. Then the following identity holds:

[TABLE]

where for two random variables $\boldsymbol{X}$ and $\boldsymbol{Y}$, $\mathbb{C}\mathrm{ov}\{\boldsymbol{X},\boldsymbol{Y}\}:=\mathbb{E}\{\boldsymbol{X}^{\!\top\!}\boldsymbol{Y}\}-\mathbb{E}\{\boldsymbol{X}\}^{\!\top}\mathbb{E}\{\boldsymbol{Y}\}$ and $\mathbb{V}\mathrm{ar}\{\boldsymbol{X}\}:=\mathbb{C}\mathrm{ov}\{\boldsymbol{X},\boldsymbol{X}\}$.

To gain an insight into this lemma, consider a constant arrival vector which makes the right-hand side of (25) vanished. In light of $\mathbb{C}\mathrm{ov}\{\boldsymbol{B}_{\circ}\boldsymbol{f},\boldsymbol{q}_{\circ}\}=\mathbb{C}\mathrm{ov}\{\boldsymbol{f},\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{q}_{\circ}\}$, equality (25) then implies that a stabilizing routing decision with a higher average total variance of link forwardings necessarily results in a higher average total covariance between link forwardings and link queue-differentials. For example, compared with BP that saturates activated links to their capacity limits, HD with a more conservative packet forwarding results in less variations in link actual transmissions. The lemma then claims that HD leads to a smaller correlation between link forwardings and link queue-differentials, which is confirmed by comparing HD and BP algorithms (see H4 in Sec. 3.2).

Definition 9**.**

We specify $\cal D$-class routing policies as a collection of all dynamic routing policies that make timeslot routing decisions based only on current queue occupancies and channel states and so independent of arrival statistics and channel state probabilities.

By allowing as many routes as possible, $\cal D$-class routing policies tend to distribute traffic all over the network. This class includes all opportunistic max-weight schedulers that do not incorporate the Markov structure of topology process into their decisions, including BP [43] and most of its derivations [39, 15, 38, 32, 41, 16, 12, 22, 21, 1, 31, 30, 20, 46, 11, 24, 23]. The class also encompasses all offline stationary randomized algorithms (possibly unfeasible) that make routing decisions as pure functions only of observed channel states, and so independent of queue occupancies, by typically using the knowledge of arrival statistics and channel state probabilities.

Theorem 3** (HD Minimum Delay).**

Consider a uniclass wireless network that meets Assum. 1 under a stabilizable arrival rate. Within $\cal D$-class routing policies and under the K-hop interference model, HD with $\beta=0$ minimizes the average total queue congestion $\overline{Q}$ as defined in (2), which is proportional to average network delay by Little’s Theorem.

7 Classical vs Combinatorial Heat Process

To formulate heat diffusion on graph, we use the theory of combinatorial geometry, where the notion of chains-cochains on a combinatorial domain provides a genuine counterpart for differential forms in classical geometry. Details are found in [2] and references therein.

7.1 Heat Equations on Manifolds

On a smooth manifold $\cal M$ charted in local coordinates $\boldsymbol{z}$, consider $Q(\boldsymbol{z},t)$ as spatial distribution of temperature, $\boldsymbol{F}(\boldsymbol{z},t)$ as heat flux, and $A(\boldsymbol{z},t)$ as scalar field of heat sources (with minus for sinks). The law of heat conservation entails

[TABLE]

Fick’s law relates the diffusive flux to the concentration, postulating that the heat flux goes from warm regions of high concentration to cold regions of low concentration, with a magnitude that is proportional to the concentration gradient:

[TABLE]

where $\sigma(\boldsymbol{z})$ is thermal diffusivity that quantifies how fast heat moves through the material. Combining (26) and (27) together, we obtain

[TABLE]

To have a unique solution, besides time initial condition, one must prescribe $Q$ conditions on a boundary $\partial{\cal M}$.

7.2 Heat Equations on Undirected Graphs

In the context of combinatorial geometry, let us view a graph as a simplicial 1-complex and transfer elements of classical heat equations to this cell complex. In doing so, the smooth manifold $\cal M$ is replaced by a 0-chain vector representing the discrete domain, the pointwise functions $Q(\boldsymbol{z},t)$ and $A(\boldsymbol{z},t)$ are respectively replaced by 0-cochain vectors $\boldsymbol{q}(t)$ and $\boldsymbol{a}(t)$ (node variables), the line integral $\boldsymbol{F}(\boldsymbol{z},t)$ is replaced by 1-cochain vector $\boldsymbol{f}(t)$ (edge variable), and the thermal diffusivity $\sigma$ is replaced by a vector of edge weights $\boldsymbol{\sigma}$.

As a 1-complex, the graph structure is fully described by its node-edge incidence matrix $\boldsymbol{B}$. 222 The incidence matrix defined in Sec. 4 for a directed graph has the same structure except that the edge directions are substituted for the arbitrarily assigned algebraic topological edge orientations here. Then on an undirected graph with a node $d$ as the heat sink, combinatorial analogue of the classical heat equations (26)–(28) are obtained as

[TABLE]

Note that the boundary $\partial{\cal M}$ on the manifold collapses to the single node $d$ on the graph at the fixed zero temperature, which absorbs all the heat generated by the heat sources $\boldsymbol{a}(t)$.

Enforcing boundary condition $q_{d}(t)\!=\!0$, one can eliminate the sink $d$ from (29)–(31), which yields the reduced set of continuous-time graph heat equations as

[TABLE]

where as before, subscript $\circ$ denotes a reduced vector or matrix that discards the entries corresponding to the destination node $d$. The linear operator $\boldsymbol{L}_{\circ}$ is called the Dirichlet Laplacian with respect to the node $d$, which is a symmetric and diagonally dominant matrix. Further, it can be shown that for any connected graph, $\boldsymbol{L}_{\circ}$ is positive definite.

7.3 Heat Equations on Directed Graphs

On a directed graph, the combinatorial heat conservation (29) remains unchanged, but the Fick’s law (30) must be modified to allow flow in only one direction. Let arbitrarily assigned edge orientations concur with edge directions. Like the undirected case, one can drop the sink node $d$ from equations by fixing $q_{d}(t)=0$ as boundary condition. Then we get the reduced set of continuous-time heat equations on an uncapacitated directed graph as

[TABLE]

We refer to $\boldsymbol{\vec{L}}_{\circ}$ as nonlinear Dirichlet Laplacian that acts on a directed graph and, unlike $\boldsymbol{L}_{\circ}$, is an operand-dependent operator that retains neither linearity nor symmetry.

Remark 5*.*

For the first time, heat diffusion on directed graphs is formulated via a nonlinear Laplacian. This is in agreement with the recent work in [36] showing that heat diffusion on Finsler manifolds, the natural counterparts of directed graphs in continuous domain, leads to a nonlinear Laplacian. In the graph literature, different linear Laplacians have been proposed for directed graphs (see [9, Sec. 3] for a review). While successful to address some purely graphical issues, they are not able to convey the physics of the diffusion process, nor the intrinsic nonlinearity due to the one-way flow restrictions.

Given finite heat sources, heat equations on a connected undirected graph always lead to finite temperatures at the nodes. However, for (37) to have a finite solution, each nonzero heat source needs to connect to the sink through at least one directed path. If this basic condition does not hold, the network flow problem has indeed no solution in the sense that there is no way to transfer all commodities, which is heat in our case, to the destination.

Definition 10**.**

A nonzero heat source is feasible if it connects to sink through at least one directed path, with the path being directed from source to sink for a positive heat source and from sink to source for a negative heat source. A vector of heat sources is feasible if each of its nonzero components is feasible.

8 Wireless Network Thermodynamics

Though defined on a directed graph, the heat equations (34)–(37) still represent a deterministic, continuous-time process with no link interference. The latter, particularly, makes the wireless problem quite intractable. Nonetheless, this section advocates a genuine diffusion on stochastic, slotted-time, interference networks by showing that under HD routing policy, the long-term average dynamics of an interference wireless network comply with non-interference combinatorial heat equations on a suitably-weighted directed graph.

8.1 HD Fluid Limit

Fluid limit of a stochastic process is the limiting dynamics obtained by scaling in time and amplitude. Under very mild conditions, it is shown that these scaled trajectories converge to a set of deterministic equations called fluid model. Using such a deterministic model, one can analyze rate-level, rather than packet-level, behavior of the original stochastic process. Details are found in [14, 10] and references therein.

Fluid limit: Let $\boldsymbol{X}(\omega,t)$ be a realization of a continuous-time stochastic process $\boldsymbol{X}$ along a sample path $\omega$. Define the scaled process $\boldsymbol{X}^{r}(\omega,t):=\boldsymbol{X}(\omega,rt)/r$ for any $r>0$. A deterministic function $\boldsymbol{\tilde{X}}(t)$ is a fluid limit if there exist a sequence $r$ and a sample path $\omega$ such that $\;\lim_{r\to\infty}\boldsymbol{X}^{r}(\omega,t)\to\boldsymbol{\tilde{X}}(t)$ uniformly on compact sets. For a stable flow network, the existence of fluid limits is guaranteed if exogenous arrivals are of finite variance. It is further shown that each fluid limit is Lipschitz-continuous, and so differentiable, almost everywhere with respect to Lebesgue measure on $[0,\infty)$.

Cumulative process: Note that the fluid theorem is defined for continuous-time stochastic processes, while a wireless network is a slotted-time process. To resolve this issue, we derive a first-order continuous-time approximation of wireless network dynamics using its cumulative model. Let $\boldsymbol{a}_{\circ}^{\mathrm{tot}}(n)$ and $\boldsymbol{f}^{\mathrm{tot}}(n)$ be respectively the vector of cumulative node arrivals and link transmissions up to slot $n$. In light of dynamic equation (10) and by assuming the initial conditions $\boldsymbol{a}_{\circ}^{\mathrm{tot}}(0)\!=\!\boldsymbol{0}$ and $\boldsymbol{f}^{\mathrm{tot}}(0)\!=\!\boldsymbol{0}$, we obtain

[TABLE]

Let $\smash{\widehat{f_{ij}}}(n)$ be the predicted number of packets that link $ij$ would transmit if it were activated at slot $n$ and compose the vector $\smash{\smash{\widehat{\boldsymbol{f}(}}}n)$. Also let $T_{\boldsymbol{\pi}}(n)$ be the cumulative number of timeslots, until slot $n$, in which the scheduling vector $\boldsymbol{\pi}\!\in\!\Pi$ has been selected. Assuming the initial condition $T_{\boldsymbol{\pi}}(0)=0$, one can verify that

[TABLE]

The first parenthesis equals 1 if the scheduling vector $\boldsymbol{\pi}$ has been selected at slot $k$, and 0 otherwise. The term $(\boldsymbol{\pi}\odot\smash{\smash{\widehat{\boldsymbol{f}(}}}k))$ represents the number of packets that could be transmitted over each link if the scheduling vector $\boldsymbol{\pi}$ were selected. Note that a routing policy needs to determine each entry of $\smash{\widehat{\boldsymbol{f}(}}k)$ and select a scheduling vector $\boldsymbol{\pi}\!\in\!\Pi$ at each timeslot.

General fluid equations: Given a sample path $\omega$, we extend a slotted-time process to be continuous-time via linear interpolation in each timeslot interval $(n,n+1)$. With no loss of generality, let exogenous arrivals occur at the beginning of each timeslot so that $\boldsymbol{a}_{\circ}^{\mathrm{tot}}(t)$ represents cumulative arrivals by time $t$. Assuming normalized timeslots with the period of time unit, (38) directly provides the first set of stochastic general fluid equations as

[TABLE]

with $\overline{\boldsymbol{a}_{\circ}\!}\,$ being the time average expectation of the random arrivals $\boldsymbol{a}_{\circ}(n)$. The second set of general fluid equations are obtained from the time derivative of (39) as

[TABLE]

Note that (42) entails the existence of a $\delta>0$ such that

[TABLE]

for any $t^{\prime}\!\in[t,t+\delta]$. This states the fact that if a link has a positive flow of packets at time $t$, the number of packets transmitted by the link in an interval $[t,t^{\prime}]\subset[t,t+\delta]$ is equal to the amount of time the link has been activated during $[t,t^{\prime}]$ multiplied by its transmission rate prediction at time $t$.

Particular fluid equations: While (40)–(44) hold for any stable network operating under an arbitrary non-idling control policy, each policy determines $\smash{\smash{\widehat{\boldsymbol{f}(}}t)}$ and $T_{\boldsymbol{\pi}}(t)$ in its own particular way. Referring to (6), HD policy enforces

[TABLE]

where $\overline{\boldsymbol{\Phi}}$ represents the time average expectation of $\boldsymbol{\Phi}(n)$ as defined in (11). Note that the existence of $\overline{\boldsymbol{\mu}}$ and $\overline{\boldsymbol{\rho}}$ is secured by (15) and (16). Referring now to (7) and (8), HD policy determines the scheduling vector $\boldsymbol{\pi}(t)$ by solving the following max-weight optimization problem:

[TABLE]

where $\boldsymbol{w}(t)$ is the vector of weights assigned by HD policy to each link at time $t$.

For a comparison, observe that the original BP solves the same max-weight optimization problem (46) to find a scheduling vector $\boldsymbol{\pi}(t)$, but it enforces $\smash{\widehat{\boldsymbol{f}(}}t)$ and $\boldsymbol{w}(t)$ to be

[TABLE]

Theorem 4** (HD Fluid Model).**

On a uniclass wireless network stabilized by Pareto optimal HD policy, every fluid limit \boldsymbol{\tilde{X}}(t)=\bigl{(}\boldsymbol{\tilde{q}}_{\circ}(t),\boldsymbol{\tilde{f}}{{}^{\mathrm{tot}}(t)},{\tilde{T}}_{\boldsymbol{\pi}}(t)\bigr{)} satisfies HD fluid model, which is defined as the collection of deterministic continuous-time equations (40)–(47).

Remark 6*.*

It is important to discriminate between fluid limit and fluid model of a discrete-time stochastic process. The former is the scaled process of the first-order continuous-time approximation for an arbitrary realization of the stochastic process, while the latter is a (set of) fully deterministic, continuous-time equation(s). Consider now a wireless network under HD routing policy, where the discrete-time stochastic processes $\boldsymbol{q}_{\circ}(n)$, $\boldsymbol{f}{{}^{\mathrm{tot}}(n)}$ and $T_{\boldsymbol{\pi}}(n)$ have respectively the continuous-time fluid limits $\boldsymbol{\tilde{q}}_{\circ}(t)$, $\boldsymbol{\tilde{f}}{{}^{\mathrm{tot}}(t)}$ and ${\tilde{T}}_{\boldsymbol{\pi}}(t)$. Then Th. 4 states that for large enough scaling factors, the fluid limit of every realization converges to a set of deterministic, continuous-time functions $\boldsymbol{q}_{\circ}(t)$, $\boldsymbol{f}{{}^{\mathrm{tot}}(t)}$ and $T_{\boldsymbol{\pi}}(t)$ which solve the HD fluid model equations (40)–(47).

8.2 Thermodynamic-Like Packet Routing

Consider a uniclass wireless network with packets being routed under HD policy (microscopic flow). At each timeslot, HD policy activates a particular set of links to transmit a specific number of packets over them. Obviously, each link transmits packets at some slots and is switched off at some other slots. Let us now look at the limit flow on each link, defined as the total number of packets transmitted over the link during a large period of time divided by the time duration. We claim that observing average packet flow in limit (macroscopic flow), it takes the form of heat flow on the underlying directed graph with suitably-weighted edges.

Consider a thermal graph with the same node-edge incidence matrix $\boldsymbol{B}_{\circ}$ and the edge thermal diffusivity $\sigma_{ij}=\,\overline{\!\phi_{ij}\!}\,$. Associate with each arrival $a_{i}(n)$ on the wireless network a static heat source of intensity $\overline{a_{i}}$ on the graph and fix zero temperature at the destination node. The flow of heat on this directed graph is governed by (34)–(37), which provides the wireless network with a static reference thermal model as

[TABLE]

Note that $\,\overline{\!\phi_{ij}\!}\,$ depends not only on the link cost factor $\,\overline{\!\rho_{ij}\!}\,$, but also on the penalty factor $\beta$, where varying $\beta$ leads to different edge weights and so different graph topologies.

Recall that $T_{\boldsymbol{\pi}}(t)$ represents the cumulative time until $t$ in which the scheduling vector $\boldsymbol{\pi}\in\Pi$ has been selected. Obviously, each scheduling policy leads to its own specific $T_{\boldsymbol{\pi}}(t)$. For example, under HD policy, $T_{\boldsymbol{\pi}}(t)$ is determined by the HD scheduling (46)–(47), while the original BP determines it according to (48)–(49).

Definition 11**.**

Under a sequence of wireless link scheduling, the effective capacity on each link is the time average expectation of capacity made genuinely available on that link:

[TABLE]

where $\boldsymbol{\mu}_{\mathrm{eff}}$ denotes the vector of effective link capacities.

Observe that the classical heat equations (26)–(28), and their combinatorial counterparts (32)–(33), take no limit in either flow direction or flow capacity. Then note that while (34)–(37) extend heat equations to directed graphs, they still consider no capacity limits on branches. In fact, the underneath assumption is that the flow of heat on each directed edge follows the Fick’s law of diffusion, not intervened by the edge capacity.

Assumption 2**.**

Given an arrival rate $\overline{\boldsymbol{a}_{\circ}\!}\,$, there exists at least one sequence of wireless link scheduling under which the effective link capacities meet the requirement of reference heat flow (50), which is stated by $\boldsymbol{f}^{\mathrm{opt}}\!\preccurlyeq\boldsymbol{\mu}_{\mathrm{eff}}$.

While $\boldsymbol{\mu}_{\mathrm{eff}}$ is a network characteristic and independent of arrivals, satisfaction of Assum. 2 does depend on arrivals. Further, for a given arrival rate, there could be a large number of link scheduling sequences that meet the requirement.

Theorem 5** (Wireless Network Thermodynamics).**

Consider a uniclass wireless network that meets Assum. 1 and 2 under a stabilizable arrival rate. Then the HD fluid model (40)–(47) asymptotically converges to the thermal model (50)–(51). In particular, HD fluid model with $\beta=0$ complies with heat equations on an unweighted directed graph, and with $\beta=1$ to those on a weighted directed graph with $\sigma_{ij}=1/\,\overline{\!\rho_{ij}\!}\,$.

Remark 7*.*

Assum. 2 examines if it is possible in principle to stabilize the wireless network such that its fluid limit follows uncapacitated heat equations. We fully revoke this assumption in [4] by developing diffusion equations on capacitated directed graphs and showing that the fluid equations (40)–(47) still respect them with no need of satisfying Assum. 2. In fact, we solve a more complicated diffusion problem in [4], where not only directed edges have limited capacity, but flows with different destinations need to be carried over the network, which raises the challenge of optimal designation of edge capacities to each of them.

9 HD Minimum Routing Cost at $\beta=1$

To establish the second pillar of HD Pareto optimality, this section shows, via Dirichlet’s principle, that average quadratic routing cost is minimized under HD policy with $\beta\!=\!1$. In fact, we show a more general result that HD with any $\beta\in[0,1]$ solves the following $\beta$-dependent optimization problem:

[TABLE]

where $\beta\!=\!1$ leads to $\phi_{ij}=1/\rho_{ij}$, which recovers (3) on minimizing the average quadratic routing cost $\overline{R}$.

Remark 8*.*

At $\beta\!=\!0$, we get $\phi_{ij}=0.5$ for all links, implying that $\sum\nolimits_{ij}\overline{\!(f_{ij})^{2}}$ is minimized by HD with $\beta\!=\!1$. The average total queue congestion $\overline{Q}$ is also minimized by HD with $\beta=0$ (see Th. 3). This entails that the two objective functions are minimized by the same timeslot control action $\boldsymbol{f}(n)$, which makes the ground for our results on HD weak Pareto optimality for uniform link costs in Sec. 10.

9.1 Classical Dirichlet Principle

Consider the classical heat diffusion equations (26)–(28) subject to constant heat sources $A(\boldsymbol{z})$. In a steady-state thermal conduction, the amount of heat entering any region of manifold $\cal M$ is equal to the amount of heat leaving out the region. Thus, while partial derivatives of temperature with respect to space may have either zero or nonzero values, all time derivatives of temperature at any point on $\cal M$ will remain uniformly zero. This leads to the classical Poisson equation

[TABLE]

which formulates stationary heat transfer by substituting zero for the time derivative of temperature in (28). Dirichlet’s principle then states that the Poisson equation has a unique solution that minimizes the Dirichlet energy

[TABLE]

among all twice differentiable functions $Q(\boldsymbol{z})$ that respect the boundary conditions on $\partial{\cal M}$.

9.2 Combinatorial Dirichlet Principle

To derive the combinatorial analogue of Poisson equation on undirected graphs, one identifies the classical $\mathrm{div}$ with the boundary operator $\boldsymbol{B}$ and the classical gradient $\nabla$ with the minus333 In vector calculus, the gradient of a scalar field is positive in the direction of increase of the field. On a graph, on the other hand, we take the gradient of a node variable positive in the direction of decrease of the variable. By the same reason, the classical Laplace operator is a negative semi-definite operator, while the graph Laplacian is a positive semi-definite matrix.

of coboundary operator $\boldsymbol{B}^{\top}\!$. Fixing $q_{d}(t)=0$ yields

[TABLE]

which correctly realizes (33) in steady-state. Note that the equation has no time variable $(t)$, since it represents the steady-state condition. It is not difficult to see that, like the classical case, this equation has a unique solution that minimizes the combinatorial Dirichlet energy

[TABLE]

The proof of Dirichlet’s principle is much simpler in the combinatorial case. In fact, as $\boldsymbol{L}_{\circ}$ is positive definite, $E_{D}(\boldsymbol{q}_{\circ})$ is convex and so has a minimum at the critical point where its first order variation vanishes, which readily leads to the combinatorial Poisson equation (53).

9.3 Nonlinear Dirichlet Principle

Essentially, the Poisson equation on a directed graph should capture the steady-state behavior of combinatorial nonlinear diffusion process (37) subject to constant heat sources $\boldsymbol{a}_{\circ}$. This leads to the following nonlinear Poisson equation:

[TABLE]

Difficulty arises from the fact that contrary to linear Laplacian $\boldsymbol{L}_{\circ}$ on undirected graphs that is a symmetric positive definite matrix, $\boldsymbol{\vec{L}}_{\circ}$ is an operand-dependent operator that retains neither linearity nor symmetricity. Thus, the easy way of proving Dirichlet’s principle on undirected graphs ceases to exist here, as one can not claim that $\boldsymbol{\vec{L}}_{\circ}\boldsymbol{q}_{\circ}$ in (54) is the directional derivative of $\frac{1}{2}\,\boldsymbol{q}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{\vec{L}}_{\circ}\boldsymbol{q}_{\circ}$. Nonetheless, we extend the concept of combinatorial Dirichlet principle to directed graphs by the next theorem.

Theorem 6** (Nonlinear Dirichlet Principle).**

Given a feasible $\boldsymbol{a}_{\circ}$ on a directed graph, the nonlinear Poisson equation (54) has a unique solution that minimizes the nonlinear Dirichlet energy

[TABLE]

Remark 9*.*

Though Dirichlet’s principle on undirected graphs has been known for long time, its extension to directed graphs is completely new to literature. As a model of heat flow on directed graphs, one can conceptualize a resistive network with a diode added to each edge [4]. Electrical current – the counterpart of combinatorial heat flux – moves along negative gradient of voltage, but only under the condition of respecting the diode direction. Another example is a piping network of liquid/gas with a check valve on each line. Again, the liquid/gas flows along negative gradient of pressure, while each check valve allows the flow in only one direction.

9.4 Quadratic Routing Cost Minimization

The framework of Th. 6 is not yet aligned with what we need for the optimization problem (52). The next theorem resolves this incongruity by showing that minimizing the functional (55) is indeed the dual of minimizing network energy dissipation, known as Thomson’s principle, on the directed graph with zero duality gap.

Theorem 7** (Nonlinear Thomson Principle).**

Minimizing the nonlinear Dirichlet energy (55) subject to the nonlinear Poisson equation (54) is equivalent to minimizing total energy dissipation on the graph subject to flow conservation at the nodes, stated by

[TABLE]

where $\boldsymbol{f}\!\succcurlyeq\boldsymbol{0}$ is imposed by network directionality. Further, temperatures at the nodes play the natural role of the Lagrange multipliers in the dual of the optimization problem (56).

It is worth comparing the minimization problem (56) with the celebrated law of least energy dissipation on resistive networks. In essence, Th. 7 extends the law to directed graphs, or to nonlinear resistive-diode networks for that matter [4]. The upshot is then due to the connection between heat diffusion on capacitated directed graphs and HD fluid limit, which brings together circuit theory and wireless networking under one umbrella.

Theorem 8** (HD Minimum Routing Cost).**

Consider a uniclass wireless network that meets Assum. 1 and 2 under a stabilizable arrival rate. Then HD policy solves the $\beta$-dependent optimization problem (52). In particular, HD policy with $\beta\!=\!1$ minimizes the quadratic routing cost $\overline{R}$ as defined in (3).

In light of Th. 5, every expected time average value on a stochastic wireless network governed by HD policy follows the corresponding stationary value produced by nonlinear heat equations on the suitably weighted underlying directed graph. In particular, the $\beta$-dependent objective function in (52) complies with the total energy dissipation $\vec{E}_{R}(\boldsymbol{f})$ on the graph weighted by $\sigma_{ij}=\,\overline{\!\phi_{ij}\!}\;$. By the same token, the average quadratic routing cost $\overline{R}$ complies with $\vec{E}_{R}(\boldsymbol{f})$ on the graph weighted by $\sigma_{ij}=1/\,\overline{\!\rho_{ij}\!}\;$.

Remark 10*.*

As Rem. 7 explained, Assum. 2 ensures that the link capacity constraints of wireless network do not intervene the Fick’s law on its underlying directed graph. Again, this assumption is fully revoked in [4, 7] by developing Dirichlet’s principle on capacitated directed graphs and showing that HD fluid model still complies with it.

10 Pareto Optimality

Minimizing average network delay and minimizing average routing cost are often conflicting objectives, meaning that as one decreases the other has to increase. This naturally leads to a multi-objective optimization framework. Then the favorite operating points lie on the Pareto boundary that corresponds to equilibria from which any deviation results in performance degradation in at least one objective. In other words, a Pareto optimal solution is a state of allocation of resources from which it is impossible to reallocate so as to make any one objective better off without making at least another objective worse off.

10.1 Strong Pareto Optimality for Nonuniform Link Costs

We have shown that HD with $\beta=0$ minimizes the average network delay $\overline{Q}$ among all $\cal D$-class routing policies – solving the optimization problem (2). We have also shown that HD with ${\beta=1}$ strictly minimizes the quadratic routing cost $\overline{R}$ among all stabilizing routing algorithms – solving the optimization problem (3). Consider now the region of operation built on joint variables $(\overline{Q},\overline{R})$ in which $\overline{Q}$ is achievable by $\cal D$-class routing policies (possibly unfeasible). The next theorem shows that HD policy operates on the Pareto boundary of this $(\overline{Q},\overline{R})$ region by altering $\beta\in[0,1]$ – solving the multi-objective optimization problem (4).

Theorem 9** (HD Pareto Optimality).**

Consider a uniclass wireless network that meets Assum. 1 and 2 under a stabilizable arrival rate and the K-hop interference model. Suppose that the operating region built on all possible joint variables $(\overline{Q}\,,\overline{R})$ with $\overline{Q}$ produced by a $\cal D$-class routing policy is convex. Then HD policy operates on the Pareto boundary of $(\overline{Q}\,,\overline{R})$ region by altering $\beta\in[0,1]$.

It is worth to note that in the case of non-convex Pareto boundary, HD with $\beta\in[0,1]$ still covers the points on convex parts of the boundary, though some Pareto optimal points lie on non-convex parts [40, 26, 25]. The convexity is jeopardized in the presence of a positive correlation between $\overline{Q}$ and $\overline{R}$, e.g., if the routing cost is defined as it could grow by the increase of queue occupancy.

Remark 11*.*

To the best of our knowledge, this is the first time a network layer routing policy provides Pareto optimal performance with respect to average network delay and routing cost, without requiring any knowledge of traffic and topology.

10.2 Weak Pareto Optimality for Uniform Link Costs

Recall from (6) that $\phi_{ij}(n)=(1\!-\!\beta)/\vartheta_{ij}+\beta/\rho_{ij}(n)$ with $\vartheta_{ij}=1$ if node $j$ is the final destination, i.e., $j=d$, and $\vartheta_{ij}=2$ otherwise. When the cost factors in all links converge to $\vartheta_{ij}$, we get $\phi_{ij}=1/\vartheta_{ij}$ for any $\beta$, and so the performance of HD policy turns to be independent of the penalty factor $\beta$. Considering this observation along with Th. 3 implies that the average network delay $\overline{Q}$ must be minimized for all $\beta\in[0,1]$. Considering it along with Th. 8, on the other hand, implies that the routing cost $\overline{R}$ must also be minimized for all $\beta\in[0,1]$. Holding these two requirements at the same time entails that $\overline{Q}$ and $\overline{R}$ must be minimized together, which equivalently means that the Pareto boundary of $(\overline{Q}\,,\overline{R})$ region must shrink into one single point. Such an operating point is called weakly Pareto optimal in the sense that no tradeoff is allowed as it is impossible to strictly improve at least one operating objective. The upshot is formalized by the next corollary.

Corollary 1**.**

Consider a uniclass wireless network under the same condition of Th. 9. Suppose the cost factors for all wireless links converge to $\vartheta_{ij}$, defined in (6). Then the Pareto boundary of $(\overline{Q}\,,\overline{R})$ region shrinks to a point at which HD policy operates for all $\beta\in[0,1]$.

Under uniform cost factor condition for all links, Fig. 7 provides a graphical illustration of the feasible region built on $(\overline{Q}\,,\overline{R})$. It emphasizes HD operation at the weakly Pareto optimal point for all $\beta\in[0,1]$ in comparison with the performance of V-parameter BP for $V\in[0,\infty)$.

11 Conclusion

We have introduced a network layer routing policy, called Heat-Diffusion (HD), for uniclass wireless networks that (i) is throughput optimal, (ii) minimizes average quadratic routing cost, (iii) minimizes average network delay within an important class of routing policies, (iv) provides a Pareto optimal tradeoff between average network delay and quadratic routing cost, and (v) enjoys the same algorithmic structure, complexity and overhead as Back-Pressure (BP) routing policy. Further, HD policy is strongly connected to the world of heat calculus in mathematics, which we believe opens the door to a rich array of theoretical techniques to analyze and optimize wireless networking. For example, such a connection provides a new way of analyzing the impact of wireless network topology on stability and capacity region [44] or on delay/routing energy performance [45]. A decentralized HD protocol has been pragmatically implemented and experimentally evaluated in [19] for data collection in wireless sensor networks, including a comparative analysis of its performance with respect to the Backpressure Collection Protocol [31]. In [37], a HD-based delay-aware framework is designed for joint dynamic routing and link-scheduling in multihop wireless networks.

Though motivated by wireless networks, the HD framework can be extended in various ways to other application areas. Among them is packet scheduling in high speed switches with a lot of attention in recent years. Resource allocation problems in manufacturing and transportation also fall within the scope of the model we considered here.

Acknowledgment

The authors thank Professor Michael Neely at USC, for his helpful feedback and discussions on this work.

Appendix A

Proof of Theorems and Lemmas

Note that in the proofs we often drop timeslot variable $(n)$ for ease of notation and concision.

Proof of Theorem 1 (HD Key Property)

One can verify that

[TABLE]

Let us temporarily relax all network constraints. Then each link-related component of $D(\boldsymbol{f},\boldsymbol{q}_{\circ},n)$ turns to be strictly concave. For each link $ij$, by taking the first derivative with respect to $f_{ij}$, we find the maximizing link transmission $f_{ij}^{\mathrm{opt}}\!=\phi_{ij}q_{ij}$. Considering the link constraints that $f_{ij}$ must be non-negative and at most equal to the link capacity yields

[TABLE]

which follows $\smash{\widehat{f_{ij}}}$ in (6). Considering the link interference constraint, on the other hand, enforces to activate the links that contribute most to the $D$ maximization. Then assuming that an interference model does not let a node transmit to more than one neighbor at the same time, the latter directly leads to the max-weight scheduling (8) alongside the HD weighting (7), which concludes the proof. $\blacksquare$

Proof of Lemma 1

Define $\boldsymbol{\Delta}:=\boldsymbol{B}_{\circ}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}$ and $\boldsymbol{\Delta}_{\phi}:=\boldsymbol{B}_{\circ}\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}$, which are both positive definite matrices. Since $\boldsymbol{\Delta}_{\phi}^{\!1/2}\boldsymbol{\Delta}^{\!-1}\boldsymbol{\Delta}_{\phi}^{\!1/2}$ is congruent to $\boldsymbol{\Delta}^{\!-1}$ which has only positive eigenvalues, by Sylvester’s law of inertia, $\boldsymbol{\Delta}_{\phi}^{\!1/2}\boldsymbol{\Delta}^{\!-1}\boldsymbol{\Delta}_{\phi}^{\!1/2}$ has only positive eigenvalues too. The latter is similar to $\boldsymbol{M}_{\!\circ}$, and so they have the same eigenvalues, proving that $\boldsymbol{M}_{\!\circ}$ has only positive eigenvalues.

We now show that $\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\!\circ}\,\boldsymbol{x}\geqslant 0$. Letting $\boldsymbol{v}:=\boldsymbol{\Delta}^{\!-1}\boldsymbol{x}$ and substituting for $\boldsymbol{M}_{\!\circ}$, it suffices to show that

[TABLE]

Doing another change of variable, let $\boldsymbol{f}:=\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{v}$ that represents an edge vector in which $f_{ij}=v_{i}-v_{j}$, $\forall\,ij\in\cal E$. Recall that $\boldsymbol{B}_{\circ}$ is a signed node-edge incidence matrix with arbitrarily assigned algebraic topological edge orientations. Let us assign edge orientations such that $f_{ij}\!\geqslant 0$, $\forall\,ij\in\cal E$. Then to fulfill (57), it suffices to show that $\boldsymbol{f}^{\top}\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{B}_{\circ}\,\boldsymbol{f}\geqslant 0$ subject to $\boldsymbol{f}\succcurlyeq\boldsymbol{0}$, which reads $f_{ij}\!\geqslant 0$, $\forall\,ij\in\cal E$. To this end, we equivalently show that minimum cost in the following optimization problem is non-negative:

[TABLE]

Let us construct the Lagrangian dual problem

[TABLE]

with $\boldsymbol{\lambda}$ being the vector of Lagrange multipliers. Since the primal variable $\boldsymbol{f}$ is continuously differentiable, so the Lagrangian $\cal L$, and thus the minimum occurs where $\nabla_{\!\boldsymbol{f}}\,{\cal L}=\boldsymbol{0}$, which leads to

[TABLE]

Substituting $\boldsymbol{f}^{\mathrm{opt}}$ in (58) and noting that both $\boldsymbol{f}^{\mathrm{opt}}$ and $\boldsymbol{\lambda}$ are entrywise non-negative, we obtain

[TABLE]

By the weak duality theorem, the minimum of the primal problem is greater than or equal to the maximum of the dual problem. Thus, (59) entails $\min_{\boldsymbol{f}\succcurlyeq\boldsymbol{0}}\,\boldsymbol{f}^{\top}\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{B}_{\circ}\boldsymbol{f}\!\geqslant 0$, which equally means $\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\!\circ}\,\boldsymbol{x}\geqslant 0$.

It remains to show that $\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\!\circ}\,\boldsymbol{x}=0$ only if $\boldsymbol{x}=0$, which is equivalent to show that matrix $\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\!+\boldsymbol{M}_{\!\circ}$ is positive definite. Since $\boldsymbol{x}^{\!\top}(\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\!+\boldsymbol{M}_{\!\circ})\,\boldsymbol{x}\!\geqslant 0$ already guarantees that $\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\!+\boldsymbol{M}_{\!\circ}$ is positive semi-definite, it suffices to show that $\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\!+\boldsymbol{M}_{\!\circ}$ has no zero eigenvalue. Let us assume it does, which implies the existence of an eigenvector $\boldsymbol{\nu}\neq\boldsymbol{0}$ such that

[TABLE]

Because $\boldsymbol{M}_{\!\circ}$ is the product of two positive definite matrices, $\boldsymbol{\nu}\neq\boldsymbol{0}$ entails $\boldsymbol{M}_{\!\circ}\boldsymbol{\nu}\neq\boldsymbol{0}$, which leads to $(\boldsymbol{M}_{\!\circ}\boldsymbol{\nu})^{\!\top}\boldsymbol{M}_{\!\circ}\boldsymbol{\nu}+(\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{\nu})^{\!\top}\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{\nu}>0$. Utilizing (60) in the latter results in

[TABLE]

which is not true as $\smash{(\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\!-\boldsymbol{M}_{\!\circ})^{2}}$ is a symmetric positive semi-definite matrix. Therefore, $\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\!+\boldsymbol{M}_{\!\circ}$ has no zero eigenvalue and so is symmetric positive definite. $\blacksquare$

Proof of Lemma 2

By definition of $\boldsymbol{M}_{\!\circ}$, we already have $\boldsymbol{B}_{\circ}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{M}_{\!\circ}\,\boldsymbol{x}=\boldsymbol{B}_{\circ}\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{x}$, which could easily be seen by substituting $\boldsymbol{M}_{\!\circ}$ from (18). Thus, to prove the claim, it suffices to show that for any vectors $\boldsymbol{x}$ and $\boldsymbol{y}$, equality $\boldsymbol{B}_{\circ}\boldsymbol{y}=\boldsymbol{B}_{\circ}\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{x}$ entails $\boldsymbol{y}=\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{x}$. To this end, we utilize the properties of heat equations on undirected graphs (see Sec. 7.2).

Consider a thermal graph with reduced node-edge incidence matrix $\boldsymbol{B}_{\circ}$ and edge thermal diffusivity matrix $\boldsymbol{\Phi}$ and let the destination node be fixed at zero temperature. As the first scenario, let us envision $\boldsymbol{y}$ as the vector of heat fluxes through the branches, implying that $\boldsymbol{B}_{\circ}\boldsymbol{y}$ represents the vector of heat sources injected at the nodes (see (32) and (33) under constant heat sources.) As the second scenario, envision $\boldsymbol{x}$ as the vector of temperatures at the nodes, implying that $\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{x}$ represents the vector of heat fluxes through the branches and $\boldsymbol{B}_{\circ}\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{x}$ represents the vector of heat sources injected at the nodes.

Assuming $\boldsymbol{B}_{\circ}\boldsymbol{y}=\boldsymbol{B}_{\circ}\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{x}$ means that the thermal graph is charged by the same configuration of heat sources in both scenarios above. It follows that the vector of temperatures at the nodes are also the same as the Dirichlet Laplacian is positive definite in (33). Hence, in both scenarios the vector of heat fluxes through the branches must be equal, because $\boldsymbol{B}_{\circ}$ has full row rank in (32). This entails $\boldsymbol{y}=\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{x}$, concluding the proof. $\blacksquare$

Proof of Lemma 3

Let us replace $\boldsymbol{M}_{\!\circ}+\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}$ by $2\boldsymbol{M}_{\!\circ}+(\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\!-\boldsymbol{M}_{\!\circ})$. Doing some matrix manipulation, we need to show that there exists such $1\leqslant\eta\leqslant 3$ that for arbitrary vectors $\boldsymbol{x}$ and $\boldsymbol{y}$,

[TABLE]

To this end, it suffices to show $\left|\boldsymbol{x}^{\!\top}(\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\!-\boldsymbol{M}_{\!\circ})\,\boldsymbol{y}\right|\leqslant\left|\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\!\circ}\,\boldsymbol{y}\right|$, which then makes the inequality (61) true for $\eta=1$ in case of $\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\!\circ}\,\boldsymbol{y}\leqslant 0$, and for $\eta=3$ in case of $\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\!\circ}\,\boldsymbol{y}>0$. This is equivalent to show that the following inequality holds:

[TABLE]

By little algebra, the latter can be rephrased as

[TABLE]

To prove the above inequality, it suffices to show that the minimum objective value in the following optimization problem is non-negative:

[TABLE]

where the constraints are enforced in light of Lem. 1. The Lagrangian dual problem, with $\lambda_{x}$ and $\lambda_{y}$ as the Lagrange multipliers, is found as

[TABLE]

Imposing the first order conditions $\nabla_{\!\boldsymbol{x}}\,{\cal L}=\boldsymbol{0}$ and $\nabla_{\!\boldsymbol{y}}\,{\cal L}=\boldsymbol{0}$ leads to

[TABLE]

Let us plug these two equations into the Lagrangian ${\cal L}$ and utilize the identities $\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{y}=\boldsymbol{y}^{\!\top}\boldsymbol{M}_{\!\circ}\boldsymbol{x}:=a$ and $\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\!\circ}\boldsymbol{y}=\boldsymbol{y}^{\!\top}\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{x}:=b$ with $a$ and $b$ being scalars. One can easily confirm the following identities:

[TABLE]

Then by little algebra, the Lagrangian can be transformed to

[TABLE]

Since $\boldsymbol{M}_{\!\circ}\!-\!\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}$ is skew-symmetric, both $\boldsymbol{x}^{\!\top}(\boldsymbol{M}_{\!\circ}-\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}})\boldsymbol{x}$ and $\boldsymbol{y}^{\!\top}(\boldsymbol{M}_{\!\circ}-\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}})\boldsymbol{y}$ vanish. In light of $\boldsymbol{x}^{\!\top}\boldsymbol{M}_{\!\circ}\boldsymbol{x}>0$ and $\boldsymbol{y}^{\!\top}\boldsymbol{M}_{\!\circ}\boldsymbol{y}>0$, the Lagrangian dual problem reads

[TABLE]

This entails $\boldsymbol{x}^{\!\top}(2\,\boldsymbol{M}_{\!\circ}\!-\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}})\,\boldsymbol{y}\boldsymbol{y}^{\!\top}\boldsymbol{M}_{\!\circ}\,\boldsymbol{x}\geqslant 0$ by the weak duality theorem that the maximum of the dual problem provides a lower bound for the minimum of the primal problem. $\blacksquare$

Proof of Theorem 2 (HD Throughput Optimality)(Cont.)

To simplify the proof, we assume arrivals to be i.i.d. over timeslots. For non-i.i.d. arrivals with stationary ergodic processes of finite mean and variance, a similar analysis can be done using $N$-slot Lyapunov drift [18], where the queue evolution (10) is modified to

[TABLE]

One can view $N$ as the time required for the system to reach “near steady state,” noting that in the i.i.d. case, the steady state is reached on each and every timeslot, and so $N=1$.

Back to the proof for i.i.d. arrivals, suppose that $\overline{\boldsymbol{a}_{\circ}\!}\,$ is interior to the network capacity region $\cal C$. Thus, there exists an $\epsilon>0$ such that $\overline{\boldsymbol{a}_{\circ}\!}\,+\epsilon\boldsymbol{1}\in\cal C$. Since the stationary randomized algorithm that generates $\boldsymbol{f^{\prime}}(n)$ is throughput optimal [18], it can stabilize the arrival $\overline{\boldsymbol{a}_{\circ}\!}\,+\epsilon\boldsymbol{1}$ at each timeslot. The i.i.d. assumption on arrivals then leads to

[TABLE]

implying that both $\boldsymbol{a}_{\circ}$ and $\boldsymbol{f}^{\boldsymbol{\prime}}$ reach their steady states on each and every timeslot. Plugging this into the Lyapunov drift inequality (23) yields

[TABLE]

Let us assume that there exists a $\mu\!>0$, which is explored later, such that $\boldsymbol{1}^{\!\top}\mathbb{E}\{\boldsymbol{M}_{\!\circ}\}\,\boldsymbol{q}_{\circ}\geqslant\mu\,\boldsymbol{1}^{\!\top}\!\boldsymbol{q}_{\circ}$. Using this in the latter drift inequality leads to

[TABLE]

Thus, $\mathbb{E}\{\Delta W\,|\boldsymbol{q}_{\circ}\}<0$ for any \sum_{i}q_{i}>\Gamma_{\!\max}\bigl{/}(\eta\,\mu\,\epsilon). Then in light of Theorem 2 in [27], the queuing system is stable and so $\overline{\boldsymbol{a}_{\circ}\!}\,$ is in the HD stability region. This implies that any arrival rate $\overline{\boldsymbol{a}_{\circ}\!}\,$ being interior to the network capacity region is stabilized by HD with any $\beta\in[0,1]$, meaning that HD is throughput optimal for all $\beta\in[0,1]$.

We now show that there exists such a $\mu\!>0$ that satisfies $\boldsymbol{1}^{\!\top}\mathbb{E}\{\boldsymbol{M}_{\!\circ}\}\,\boldsymbol{q}_{\circ}\geqslant\mu\,\boldsymbol{1}^{\!\top}\!\boldsymbol{q}_{\circ}$. Let us temporarily ignore the expectation and solve the problem for $\boldsymbol{M}_{\!\circ}$. The claim is trivial for $\boldsymbol{q}_{\circ}=\boldsymbol{0}$, and so we assume $\boldsymbol{q}_{\circ}\neq\boldsymbol{0}$. Further, $\boldsymbol{q}_{\circ}$ represents the vector of queue occupancies on the wireless network that are always non-negative. Let $\|\boldsymbol{q}_{\circ}\|_{1}$ represent the $\ell^{1}$ norm of $\boldsymbol{q}_{\circ}$, defined as the sum of all queue occupancies. With no loss of generality, one may normalize $\|\boldsymbol{q}_{\circ}\|_{1}$ to one. The problem can then be rephrased as finding a $\mu>0$ such that

[TABLE]

The latter is a standard linear programming problem. Using simplex method, the minimum lies on a vertex of the simplex, where the vertices of the simplex are the natural basis elements $\boldsymbol{e}_{j}:j=1,...,|{\cal V}|$. Thus, the $\mu$ is to be sought such that

[TABLE]

This immediately leads to $\mu=\min_{j}\,\boldsymbol{1}^{\!\top}\boldsymbol{M}_{\circ}\,\boldsymbol{e}_{j}$.

It remains to show that $\boldsymbol{1}^{\!\top}\boldsymbol{M}_{\circ}\,\boldsymbol{e}_{j}>0$ for every natural basis $\boldsymbol{e}_{j}$. By Lem. 3, there exists such a $1\leqslant\eta\leqslant 3$ as

[TABLE]

which implies $(\eta-1)\,\boldsymbol{1}^{\!\top}\boldsymbol{M}_{\!\circ}\,\boldsymbol{e}_{j}\geqslant\boldsymbol{1}^{\!\top}\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\,\boldsymbol{e}_{j}$. The right hand side is always positive by the next electrical circuit argument, which implies $(\eta-1)\,\boldsymbol{1}^{\!\top}\boldsymbol{M}_{\circ}\,\boldsymbol{e}_{j}>0$. The latter entails $\eta>1$ and $\boldsymbol{1}^{\!\top}\boldsymbol{M}_{\circ}\,\boldsymbol{e}_{j}>0$ as we desired.

To argue that $\boldsymbol{1}^{\!\top}\boldsymbol{M}_{\circ}{\!\!\!{}^{\top}}\,\boldsymbol{e}_{j}=\boldsymbol{e}_{j}^{\top}\boldsymbol{M}_{\!\circ}\,\boldsymbol{1}>0$, by substituting $\boldsymbol{M}_{\!\circ}$ from (18), we need to show

[TABLE]

Let us associate node-edge incidence matrix $\boldsymbol{B}_{\circ}$ with a resistive network and $\boldsymbol{e}_{j}$ with the vector of independent current sources attached to the nodes. Then the vector $\boldsymbol{v}:=(\boldsymbol{B}_{\circ}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}){{}^{-1}}\boldsymbol{e}_{j}$ reads the voltages induced at the nodes. Since $\boldsymbol{e}_{j}$ implies that electrical current is injected into the network by a single current source at the node $j$, the resulting voltage at each node is non-negative ($\boldsymbol{v}\succcurlyeq\boldsymbol{0}$). Further, the voltages at the node $j$ and at least at one of the nodes neighbor to ground (destination node) are always positive. On the other hand, the elements of each row of the Dirichlet Laplacian $\boldsymbol{B}_{\circ}\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}$ sum to zero, except for those rows representing the nodes neighbor to ground, which always sum to a positive value. (Recall that $\boldsymbol{B}_{\circ}$ is obtained from $\boldsymbol{B}$ by discarding the row related to ground.) This implies that in the vector $\boldsymbol{u}:=(\boldsymbol{B}_{\circ}\boldsymbol{\Phi}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}})\,\boldsymbol{1}$, the components related to the nodes neighbor to ground are positive, and others are zero. Considering the conditions of $\boldsymbol{u}$ and $\boldsymbol{v}$ together, we get $\boldsymbol{e}_{j}^{\top}\boldsymbol{M}_{\!\circ}\,\boldsymbol{1}=\boldsymbol{v}^{\top}\boldsymbol{u}>0$. Replacing $\boldsymbol{M}_{\!\circ}$ by $\mathbb{E}\{\boldsymbol{M}_{\!\circ}\}$, the same argument leads to $\mu=\min_{j}\,\boldsymbol{1}^{\!\top}\mathbb{E}\{\boldsymbol{M}_{\!\circ}\}\,\boldsymbol{e}_{j}>0$, which concludes the proof. $\blacksquare$

Proof of Lemma 4

Considering the maximum of $G(\boldsymbol{f})$ subject to $\boldsymbol{f}\succcurlyeq\boldsymbol{0}$, the Lagrangian dual problem is obtained as

[TABLE]

with $\boldsymbol{\lambda}$ being the vector of Lagrange multipliers. Using the first order condition $\nabla_{\!\boldsymbol{f}}\,{\cal L}=\boldsymbol{0}$, we get

[TABLE]

Plugging $\boldsymbol{\lambda}$ from (64) into the Lagrangian ${\cal L}$ leads to

[TABLE]

By the weak duality theorem, the maximum of the primal problem is smaller than or equal to the minimum of the dual problem. Further, the duality gap is zero as ${\cal L}(\boldsymbol{\lambda},\boldsymbol{f})$ is a convex functional, which leads to

[TABLE]

Substituting the latter into the $G$ functional entails

[TABLE]

One can verify, on the other hand, that $\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\!\boldsymbol{B}_{\circ}\boldsymbol{f}$ is an edge vector, in which the entry corresponding to edge $ij$ reads

[TABLE]

Under the K-hop interference model, two wireless links that share a common node cannot be scheduled in the same timeslot. Thus, for each scheduled link $ij\in{\cal E}$, we get the net outflow for node $i$ as $(\boldsymbol{B}_{\circ}\boldsymbol{f})_{i}=f_{ij}$. When node $j$ is not the final destination, we get $(\boldsymbol{B}_{\circ}\boldsymbol{f})_{j}=-f_{ij}$, and when it is, we get $(\boldsymbol{B}_{\circ}\boldsymbol{f})_{j}=0$. (Recall that $\boldsymbol{B}_{\circ}$ is a reduction of $\boldsymbol{B}$ by discarding the row related to the final destination.) Using these identities in (66), we get $(\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{B}_{\circ}\boldsymbol{f})_{ij}=\vartheta_{ij}f_{ij}$ with $\vartheta_{ij}$ defined in (6). Substituting the latter into (65), we obtain

[TABLE]

Considering the link constraints that $f_{ij}$ must be non-negative and at most equal to the link capacity yields

[TABLE]

which follows $\smash{\widehat{f_{ij}}}$ in (6) with $\beta=0$. Next is to activate the links that contribute most to the $G$ maximization that directly leads to the max-weight scheduling (8) alongside the HD weighting (7) with $\beta=0$, concluding the proof. $\blacksquare$

Proof of Lemma 5

Consider $W(n):=\boldsymbol{q}_{\circ\!}(n)^{\!\top}\boldsymbol{q}_{\circ\!}(n)$ as the classical quadratic Lyapunov candidate and take expectation from the Lyapunov drift $\Delta W(n)=W(n+1)-W(n)$ to obtain

[TABLE]

where the equality holds at each timeslot and expectation is with respect to the randomness of arrivals, channel states and (possibly) routing decision. Let $\boldsymbol{g}:=\boldsymbol{a}_{\circ}-\boldsymbol{B}_{\circ}\boldsymbol{f}+2\,\boldsymbol{q}_{\circ}$, sum over timeslots 0 until $\tau-1$, divide by $\tau$ and take a $\limsup$ of $\tau\to\infty$ from both sides of (67) to obtain the following expected time average equation:

[TABLE]

where we utilized $\limsup_{\tau\to\infty}(W(\tau)-W(0))/\tau=0$, as the routing policy stabilizes $\overline{\boldsymbol{a}_{\circ}\!}\,$ and so keeps $W(n)$ finite with probability 1 at each timeslot.

It remains to show that the left-hand side of (68) vanishes. Observe that $\boldsymbol{g}(n)$ is entrywise non-negative and finite. Thus, there exist constant vectors $\boldsymbol{g}_{\min}$ and $\boldsymbol{g}_{\max}$ such that

[TABLE]

Hence, the left-hand side of (68) is bounded from below to $(\,\overline{\boldsymbol{a}_{\circ}\!}\,-\boldsymbol{B}_{\circ}\overline{\boldsymbol{f}}\,)^{\!\top}\boldsymbol{g}_{\min}$ and from above to $(\,\overline{\boldsymbol{a}_{\circ}\!}\,-\boldsymbol{B}_{\circ}\overline{\boldsymbol{f}}\,)^{\!\top\!}\boldsymbol{g}_{\max}$. Further, as $\overline{\boldsymbol{a}_{\circ}\!}\,$ is stabilized by the routing policy, the feasibility condition in (17) entails $\overline{\boldsymbol{a}_{\circ}\!}=\boldsymbol{B}_{\circ}\overline{\boldsymbol{f}}$, implying that the left-hand side of (68) vanishes. $\blacksquare$

Proof of Theorem 3 (HD Minimum Delay)

To simplify the proof, we assume arrivals are i.i.d. over timeslots, with the understanding that it can easily be modified to yield similar result for non-i.i.d. arrivals, using the $N$-slot analysis derived from (62).

Consider an arrival rate $\overline{\boldsymbol{a}_{\circ}\!}\,$ interior to the stability region of a $\cal D$-class routing policy, which we refer to it as “generic”. Let the timeslot quantities $\boldsymbol{f}(n)$ and $\boldsymbol{q}_{\circ}(n)$ be produced by such a generic routing policy. If this generic routing policy also maximizes the $G$ functional (24) at each slot $n$, by Assum. 1, it will result in the same $\overline{Q}$ as that of HD policy at $\beta=0$. Thus, we assume the $G$ obtained by the generic policy is not maximal. Then for a sufficiently small $\epsilon>0$, there exists a routing algorithm (possibly unfeasible) that can stabilize the arrival $\overline{\boldsymbol{a}_{\circ}\!}\,+\epsilon\boldsymbol{1}$ while making $G(\boldsymbol{f},\boldsymbol{q}_{\circ},n)$ not less than that of the generic routing policy at each slot $n$. Let us refer to this algorithm as “fictitious,” as we do not intend to know how it really works. To rest assure that such an algorithm exists, one may endow it with the ability of perfectly predicting all future events with no uncertainty.

Let $\boldsymbol{f^{\prime}}(n)$ represent the vector of link actual transmissions produced by the fictitious algorithm at slot $n$ given $\boldsymbol{q}_{\circ}(n)$. Taking expectation from $G(\boldsymbol{f^{\prime}},\boldsymbol{q}_{\circ},n)\geqslant G(\boldsymbol{f},\boldsymbol{q}_{\circ},n)$ and considering $\mathbb{E}\{\boldsymbol{B}_{\circ}{\boldsymbol{f}^{\prime}}\}=\mathbb{E}\{\boldsymbol{B}_{\circ}\boldsymbol{f}\}+\epsilon\boldsymbol{1}$ due to the feasibility condition (17) and the i.i.d. arrivals, we obtain

[TABLE]

which holds for each timeslot. Summing over timeslots 0 until $\tau-1$, dividing by $\tau$ and taking a $\limsup$ of $\tau\to\infty$ from both sides lead to the following expected time average inequality:

[TABLE]

Let us exploit Lem. 5 in the second and third lines and apply the identities $\mathbb{C}\mathrm{ov}\{\boldsymbol{a}_{\circ}+\epsilon\boldsymbol{1},\boldsymbol{q}_{\circ}\!-\!\boldsymbol{B}_{\circ}\boldsymbol{f^{\prime}}\}=\mathbb{C}\mathrm{ov}\{\boldsymbol{a}_{\circ},\boldsymbol{q}_{\circ}\!-\!\boldsymbol{B}_{\circ}\boldsymbol{f^{\prime}}\}$ and $\mathbb{V}\mathrm{ar}\{\boldsymbol{a}_{\circ}+\epsilon\boldsymbol{1}\}=\mathbb{V}\mathrm{ar}\{\boldsymbol{a}_{\circ}\}$ to obtain

[TABLE]

Since $\boldsymbol{f}$ produced by the generic routing policy is independent of arrival statistics, we get $\mathbb{C}\mathrm{ov}\{\boldsymbol{a}_{\circ},\boldsymbol{B}_{\circ}\boldsymbol{f}\}=0$. Replacing $\boldsymbol{1}\!^{\top}\,\overline{\!\boldsymbol{q}_{\circ}\!}\,$ by the $\overline{Q}$ expression as defined in (2), we then obtain

[TABLE]

Consider this time HD policy at $\beta=0$ with the timeslot quantities of $\boldsymbol{q}_{\circ}^{\star}(n)$ and $\boldsymbol{f}^{\star}(n)$. Let again $\boldsymbol{f^{\prime}}(n)$ be produced by the fictitious algorithm at each slot $n$ to stabilize the arrival $\overline{\boldsymbol{a}_{\circ}\!}\,+\epsilon\boldsymbol{1}$, but this time, given $\boldsymbol{q}_{\circ}^{\star}(n)$. In light of Lem. 4, $G(\boldsymbol{f^{\prime}}\!,\boldsymbol{q}_{\circ}^{\star},n)\leqslant G(\boldsymbol{f}^{\star}\!,\boldsymbol{q}_{\circ}^{\star},n)$ at each slot $n$. Performing the similar steps of taking expectation, exploiting $\mathbb{E}\{\boldsymbol{B}_{\circ}\boldsymbol{f^{\prime}}\}=\mathbb{E}\{\boldsymbol{B}_{\circ}\boldsymbol{f}^{\star}\}+\epsilon$, translating the results into the expected time average form, using the fact that $\mathbb{C}\mathrm{ov}\{\boldsymbol{a}_{\circ},\boldsymbol{B}_{\circ}\boldsymbol{f}^{\star}\}=0$ as $\boldsymbol{f}^{\star}$ is independent of arrival statistics, and applying Lem. 5 by knowing that HD policy is throughput optimal and so stabilizes $\overline{\boldsymbol{a}_{\circ}\!}\,$, we obtain

[TABLE]

Comparing (69) and (70) along with $\epsilon>0$ lead to $\overline{Q^{\star}\!}\leqslant\overline{Q}$. This means the average network delay under HD policy with $\beta=0$ remains less than or equal to that under any other $\cal D$-class routing policy, which was called “generic” here. $\blacksquare$

Proof of Theorem 4 (HD Fluid Model)

The proof follows the exact same line of argument proposed in [14, Theorem 2.3.2] and [10, Proposition 4.12]. $\blacksquare$

Proof of Theorem 5 (Wireless Network Thermodynamics)

Let $\boldsymbol{q}_{\circ}^{\star}(t)$ and $\boldsymbol{f}^{\star}(t)$ denote the HD fluid model variables. Consider the continuous-time Lyapunov function

[TABLE]

where \overline{\boldsymbol{M}}_{\!\circ}=\bigl{(}\boldsymbol{B}_{\circ}\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\bigr{)}{{}^{-1}}\boldsymbol{B}_{\circ}\,\overline{\boldsymbol{\Phi}}\,\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}} represents the time average expectation of matrix $\boldsymbol{M}_{\!\circ}(n)$ as defined in (18). Taking time derivative from $Y(t)$, we obtain

[TABLE]

Exploiting Lem. 3 in the latter leads to

[TABLE]

for an $1\leqslant\eta\leqslant 3$. As a positive coefficient, $\eta$ has no impact on the Lyapunov argument and can simply be omitted, but for the sake of consistency we prefer to keep it in here.

To find an appropriate expression for $\boldsymbol{\dot{q}}_{\circ}^{\star}(t)$, let us begin by plugging (41) in (40) and taking time derivative to obtain

[TABLE]

Note in (45) that the entry of $\smash{\widehat{\boldsymbol{f}(}}t)$ corresponding to link $ij$ specifies the number of packets the link will send per unit time if it is activated at time $t$. Then $\boldsymbol{f}(t)$ identifies the vector of rate of actual transmissions realized at time $t$. Assume now that the entry of $\boldsymbol{f}(t)$ corresponding to link $ij$ at time $t$ is equal to $x\!\geqslant 0$, i.e., at time $t$ the link transmits $x$ number of packets per unit time. Then it should be obvious that the same entry of $\boldsymbol{\dot{f}}{{}^{\mathrm{tot}}}(t)$ at time $t$ must also be equal to $x$. In light of $\lim_{\delta\to 0}\boldsymbol{f}^{\mathrm{tot}}(t+\delta)=\boldsymbol{f}^{\mathrm{tot}}(t)+\delta\boldsymbol{f}(t)$, this can be explained more formally by the classical definition of limit as

[TABLE]

Further, (50)–(51) imply $\overline{\boldsymbol{a}_{\circ}\!}=\boldsymbol{\vec{L}}{{}_{\circ}^{\mathrm{opt}}}\boldsymbol{q}_{\circ}^{\mathrm{opt}}=\boldsymbol{B}_{\circ}\boldsymbol{f}^{\mathrm{opt}}\!$. Exploiting these latter identities in (72) yields

[TABLE]

Returning to the Lyapunov argument, let us substitute (73) in (71) and utilize equality (19) in Lem. 2 to obtain

[TABLE]

Multiplying both sides by two, adding and subtracting the term $\boldsymbol{f}^{\star}(t)^{\!\top}\!\boldsymbol{f}^{\star}(t)+\boldsymbol{f}^{\mathrm{opt}\top}\!\boldsymbol{f}^{\mathrm{opt}}$ on the left-hand side, and recasting the terms lead to

[TABLE]

Characterizing (74a) and (74b) on the wireless network given $\boldsymbol{q}_{\circ}^{\star}(t)$, they respectively read $-D(\boldsymbol{f}^{\star}\!,\boldsymbol{q}_{\circ}^{\star},t)$ and $D(\boldsymbol{f}^{\mathrm{opt}}\!,\boldsymbol{q}_{\circ}^{\star},t)$. Under Assum. 1 and in light of the HD fluid equations (45) and (47), the immediate result of Th. 4 is that at each time $t$, HD fluid limit maximizes the $D$ functional compared to any alternative forwarding that satisfies wireless network constraints. The $\boldsymbol{f}^{\mathrm{opt}}$ obviously meets the directionality constraints due to the structure of the reference thermal model. It also meets the capacity constraints due to Assum. 2. Hence, $\eqref{ProofTh5_8a}+\eqref{ProofTh5_8b}\leqslant 0$ which leads to

[TABLE]

We now characterize (75a) and (75b) on the reference thermal model and let

[TABLE]

It can be shown that given $\boldsymbol{q}^{\mathrm{opt}}_{\circ}$, the maximum of $H$ occurs at $\boldsymbol{f}=\boldsymbol{f}^{\mathrm{opt}}$ produced by heat flow. To see this, rephrase $H$ as

[TABLE]

where directionality constraints entail $f_{ij}\!\geqslant 0$. To maximize $H$, one then needs to assign $f_{ij}=0$ if $q_{ij}^{\mathrm{opt}}\!\leqslant 0$, and $f_{ij}\!=\,\overline{\!\phi_{ij}\!}\;q_{ij}^{\mathrm{opt}}$ otherwise. Putting this back in a matrix form, we arrive at the same expression as $\boldsymbol{f}^{\mathrm{opt}}$ in (50). Further, the maximizing $\boldsymbol{f}$ is unique by the reason that a given $\boldsymbol{q}^{\mathrm{opt}}_{\circ}$ leads to a unique $\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{q}^{\mathrm{opt}}_{\circ}$, and so to unique $q_{ij}^{\mathrm{opt}}$ components, as the matrix $\boldsymbol{B}_{\circ}$ has full row rank. From $H(\boldsymbol{f}^{\mathrm{opt}})\geqslant H(\boldsymbol{f}^{\star}(t))$, we then obtain $\eqref{ProofTh5_9a}+\eqref{ProofTh5_9b}\leqslant 0$, which by $1\leqslant\eta\leqslant 3$ yields $\dot{Y}(t)\leqslant 0$.

Let $\Omega$ be the largest invariant set in the set of all $\boldsymbol{q}_{\circ}^{\star}(t)$ trajectories for which $\dot{Y}(t)=0$. Since $Y(t)$ is a non-negative and radially unbounded function with $\dot{Y}(t)\leqslant 0$, LaSalle’s invariance principle states that every trajectory $\boldsymbol{q}_{\circ}^{\star}(t)$ asymptotically converges to $\Omega$. It remains to show that $\Omega$ contains only the trivial trajectory of $\boldsymbol{q}_{\circ}^{\star}=\boldsymbol{q}_{\circ}^{\mathrm{opt}}$. If $\dot{Y}=0$, then (75) entails $H(\boldsymbol{f}^{\star}(t))=H(\boldsymbol{f}^{\mathrm{opt}})$. We previously showed as well that $\boldsymbol{f}^{\mathrm{opt}}$ maximizes $H$ and is unique, which implies

[TABLE]

The intentionally dropped time variable $(t)$ in (76) emphasizes that $\boldsymbol{f}^{\star}(t)$ turns to be stationary by converging to $\boldsymbol{f}^{\mathrm{opt}}$, which in turn entails $\boldsymbol{q}^{\star}(t)$ being converged to a stationary $\boldsymbol{q}^{\star}$ too.

Given $\boldsymbol{q}^{\star}$, the equality (76) entails that $\boldsymbol{f}^{\mathrm{opt}}$ must maximize $D$, which implies $f^{\mathrm{opt}}_{ij}=0$ if $q_{ij}^{\star}\!\leqslant 0$, and $f^{\mathrm{opt}}_{ij}\!=\,\overline{\!\phi_{ij}\!}\;q_{ij}^{\mathrm{opt}}$ otherwise. In a matrix form, this is equivalent to \boldsymbol{f}^{\mathrm{opt}}=\overline{\boldsymbol{\Phi}}\,\max\bigl{\{}\boldsymbol{0},\>\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{q}_{\circ}^{\star}\bigr{\}}. Putting the latter against (50) leads to

[TABLE]

Consider a directed edge $ad$ with its head at the destination node, which has zero queue on the wireless network and zero temperature on the reference thermal model. By (77), $q_{a}^{\star}$ and $\boldsymbol{q}_{a}^{\mathrm{opt}}$ must be equal. Repeating this argument eventually yields $(\boldsymbol{q}_{\circ}^{\star})^{+}\!=(\boldsymbol{q}_{\circ}^{\mathrm{opt}})^{+}$, as any node with positive queue (resp. positive temperature) on the wireless network (resp. on the reference thermal model) must be connected to the destination node $d$ through a directed path. Further, observe that $\boldsymbol{q}_{\circ}^{\star}\!\succcurlyeq 0$, as queues cannot be negative in a wireless network, and $\boldsymbol{q}_{\circ}^{\mathrm{opt}}\!\succcurlyeq 0$, as temperatures cannot fall below zero in a thermal system with no negative heat source. Thus, $\boldsymbol{q}_{\circ}^{\star}=\boldsymbol{q}_{\circ}^{\mathrm{opt}}$, which together with (76) conclude the proof. $\blacksquare$

Proof of Theorem 6 (Nonlinear Dirichlet Principle)

One can verify, by the $\boldsymbol{\vec{L}}_{\circ}$ structure in (37), that

[TABLE]

where each entry of $\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{q}_{\circ}$ represents temperature-difference along the corresponding edge. Let $\boldsymbol{q}_{\circ}^{*}$ be the $\vec{E}_{D}(\boldsymbol{q}_{\circ})$ minimizing solution and let us rearrange and partition $\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{q}_{\circ}^{*}$ into positive, zero and negative components. Accordingly, $\boldsymbol{B}_{\circ}$ gets partitioned into $\boldsymbol{B}_{\scriptscriptstyle\oplus}$, $\boldsymbol{B}_{\varnothing}$ and $\boldsymbol{B}_{\scriptscriptstyle\ominus}$, which respectively contain the incidence information of edges with positive, zero and negative values in $\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{q}_{\circ}^{*}$. Likewise, $\boldsymbol{\sigma}$ gets partitioned into $\boldsymbol{\sigma}_{\scriptscriptstyle\oplus}$, $\boldsymbol{\sigma}_{\varnothing}$ and $\boldsymbol{\sigma}_{\scriptscriptstyle\ominus}$. Then at $\boldsymbol{q}_{\circ}=\boldsymbol{q}_{\circ}^{*}$, we obtain

[TABLE]

Observe that (78a) is strongly zero due to the $(\cdot){{}^{+}}$ operation. On the other hand, (78b) vanishes since $\boldsymbol{B}_{\varnothing}^{\top}\boldsymbol{q}_{\circ}^{*}=\boldsymbol{0}$. In light of $(\boldsymbol{B}_{\scriptscriptstyle\oplus}^{\top}\boldsymbol{q}_{\circ}^{*}){{}^{+}}=\boldsymbol{B}_{\scriptscriptstyle\oplus}^{\top}\boldsymbol{q}_{\circ}^{*}$, we then obtain

[TABLE]

Since $\boldsymbol{a}_{\circ}$ is feasible, each nonzero heat source connects to the sink through at least one directed path. Thus, under any flow that keeps $\boldsymbol{q}_{\circ}$ entrywise finite, the edges with positive temperature-difference build a connected graph with the node $d$. On the other hand, $\boldsymbol{q}_{\circ}^{*}$ is entrywise finite as it minimizes $\vec{E}_{D}(\boldsymbol{q}_{\circ})$, and so the corresponding edges in $\boldsymbol{B}_{\scriptscriptstyle\oplus}$ build a connected graph with the node $d$. This implies that $\boldsymbol{B}_{\scriptscriptstyle\oplus}\,\mathrm{diag}(\boldsymbol{\sigma}_{\scriptscriptstyle\oplus})\boldsymbol{B}_{\scriptscriptstyle\oplus}^{\top}$ is a positive definite matrix. Thus, the functional $\frac{1}{2}\,\boldsymbol{q}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{B}_{\scriptscriptstyle\oplus}\,\mathrm{diag}(\boldsymbol{\sigma}_{\scriptscriptstyle\oplus})\boldsymbol{B}_{\scriptscriptstyle\oplus}^{\top}\boldsymbol{q}_{\circ}-\boldsymbol{q}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{a}_{\circ}$ is strictly convex in $\boldsymbol{q}_{\circ}$ and so finds its minimum at the critical point, where its first order variation with respect to $\boldsymbol{q}_{\circ}$ vanishes. Comparing this with (79), it turns out that the minimizing $\boldsymbol{q}_{\circ}^{*}$ must satisfy

[TABLE]

Utilizing $\boldsymbol{B}_{\scriptscriptstyle\oplus}^{\top}\boldsymbol{q}_{\circ}^{*}\!=(\boldsymbol{B}_{\scriptscriptstyle\oplus}^{\top}\boldsymbol{q}_{\circ}^{*}){{}^{+}}$ and $(\boldsymbol{B}_{\scriptscriptstyle\ominus}^{\top}\boldsymbol{q}_{\circ}^{*}){{}^{+}}\!=(\boldsymbol{B}_{\varnothing}^{\top}\boldsymbol{q}_{\circ}^{*}){{}^{+}}\!=\boldsymbol{0}$, one can rephrase (80) as

[TABLE]

that recovers the nonlinear Poisson equation (54) at $\boldsymbol{q}_{\circ}^{*}$. Further, $\boldsymbol{q}_{\circ}^{*}$ is unique as it minimizes the strictly convex functional $\frac{1}{2}\,\boldsymbol{q}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{B}_{\scriptscriptstyle\oplus}\,\mathrm{diag}(\boldsymbol{\sigma}_{\scriptscriptstyle\oplus})\boldsymbol{B}_{\scriptscriptstyle\oplus}^{\top}\boldsymbol{q}_{\circ}-\boldsymbol{q}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{a}_{\circ}$, which concludes the proof. $\blacksquare$

Proof of Theorem 7 (Nonlinear Thomson Principle)

Consider (56) as the primal optimization problem and let us construct its Lagrangian dual problem as

[TABLE]

where $\boldsymbol{\lambda}\succcurlyeq 0$ is the vector of Lagrange multipliers. From the first order condition $\nabla_{\!\boldsymbol{f}}\,{\cal L}=\boldsymbol{0}$, we get $\boldsymbol{f}^{\mathrm{opt}}\!=\mathrm{diag}(\boldsymbol{\sigma})\,\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{\lambda}$. Then enforcing the constraint $\boldsymbol{f}^{\mathrm{opt}}\!\succcurlyeq\boldsymbol{0}$ leads to $\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{\lambda}\succcurlyeq\boldsymbol{0}$, which is equivalent to $\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{\lambda}=(\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}\boldsymbol{\lambda}){{}^{+}}$. Thus, we obtain

[TABLE]

Plugging this $\boldsymbol{f}^{\mathrm{opt}}$ into the Lagrangian ${\cal L}$ and utilizing the structure of $\boldsymbol{\vec{L}}_{\circ}$ in (37), we obtain

[TABLE]

Then the dual problem reads $\max_{\boldsymbol{\lambda}}\,{\cal L}(\boldsymbol{\lambda})$, which is equivalent to the following minimization problem:

[TABLE]

Further, as $\boldsymbol{f}\succcurlyeq\boldsymbol{0}$ makes a convex set and ${\cal L}(\boldsymbol{\lambda},\boldsymbol{f})$ is a convex function, the duality gap is zero, and so both the primal and dual problems result in the same optimal solution.

Comparing (82) with the nonlinear Dirichlet equation (55), it remains to show that the Lagrangian multipliers $\boldsymbol{\lambda}$ are identical to the node temperatures $\boldsymbol{q}_{\circ}$. In (81), multiplying both sides by $\boldsymbol{B}_{\circ}$ and using the $\boldsymbol{\vec{L}}_{\circ}$ expression, we obtain

[TABLE]

where the second equality comes from the constraint in the primal problem (56). Further, by Th. 6, the nonlinear Poisson equation $\boldsymbol{\vec{L}}_{\circ}\boldsymbol{q}_{\circ}=\boldsymbol{a}_{\circ}$ has a unique solution. Putting this against (83) leads to $\boldsymbol{\lambda}=\boldsymbol{q}_{\circ}$, which concludes the proof. $\blacksquare$

Proof of Theorem 8 (HD Minimum Routing Cost)

It was shown by Th. 6 that if $\overline{\boldsymbol{a}_{\circ}\!}\,$ is feasible, then under the nonlinear heat equations (50)–(51), the stationary value of the nonlinear Dirichlet energy $\vec{E}_{D}(\boldsymbol{q}_{\circ})$ is strictly minimized. It was shown by Th. 7, on the other hand, that minimizing $\vec{E}_{D}(\boldsymbol{q}_{\circ})$ is equivalent to minimizing the stationary value of total energy dissipation $\vec{E}_{R}(\boldsymbol{f})$ on the graph. Then the proof immediately follows from Th. 5 which states that under a stabilizable arrival rate $\overline{\boldsymbol{a}_{\circ}\!}\,$, HD fluid model complies with the nonlinear heat equations (50)–(51). Note that if $\overline{\boldsymbol{a}_{\circ}\!}\,$ is stabilizable, i.e., it satisfies condition (17), then its feasibility is trivial in the sense of Def. 10. $\blacksquare$

Proof of Theorem 9 (HD Pareto Optimality)

Observe that HD policy minimizes $\overline{Q}$ at $\beta\!=\!0$, minimizes $\overline{R}$ at $\beta\!=\!1$, and changes weight on these two objectives by altering $\beta$ between 0 and 1. In fact, HD transforms the two objectives of minimizing $\overline{Q}$ and $\overline{R}$ into an aggregated objective function by multiplying each objective function by a weighting factor and summing up the two weighted objective functions. Further, the weighted sum is a convex combination of objectives as the sum of weighting factors $\beta$ and $1\!-\!\beta$ equals 1. Under the assumption that the region $(\overline{Q}\,,\overline{R})$ has a convex Pareto boundary, the proof then follows by the fact that the entire boundary can be reached using the weighted-sum method [40, 26, 25] that changes the weight on the weighted convex combination of the two objective functions. $\blacksquare$

Appendix B

Graph Laplacian

Consider a connected, weighted graph with set of nodes $\cal V$, set of undirected edges $\cal E$, node-edge incidence matrix $\boldsymbol{B}$, edge weight vector $\boldsymbol{\sigma}$ and the Laplacian $\boldsymbol{L}:=\boldsymbol{B}\,\mathrm{diag}(\boldsymbol{\sigma})\,\boldsymbol{B}^{\!\top}$. One can verify that $\boldsymbol{x}^{\!\top}\boldsymbol{L}\boldsymbol{x}=\sum_{ij\in\cal E}\sigma_{ij}(x_{i}-x_{j})^{2}\geqslant 0$, which means $\boldsymbol{L}$ is positive semi-definite. Observe that $\boldsymbol{L}\boldsymbol{1}=\boldsymbol{0}$, which implies that $\boldsymbol{1}$ is an eigenvector corresponding to the smallest eigenvalue $\lambda_{1}=0$. For the second smallest eigenvalue $\lambda_{2}$, let $\boldsymbol{\nu}$ be the eigenvector orthogonal to $\boldsymbol{1}$. Thus, ${\boldsymbol{\nu}^{\!\top}\boldsymbol{1}=0}$ and $\lambda_{2}=\boldsymbol{\nu}^{\!\top}\boldsymbol{L}\boldsymbol{\nu}=\sum_{ij\in\cal E}\sigma_{ij}(\nu_{i}-\nu_{j})^{2}$. Assume $\lambda_{2}=0$. Since the graph is connected, there exists a path between every two nodes, which enforces $\boldsymbol{\nu}=c\boldsymbol{1}$ for a constant $c$. This contradicts $\boldsymbol{\nu}^{\!\top}\boldsymbol{1}=0$, and so $\lambda_{2}$ must be positive.

Let $\boldsymbol{L}^{\dagger}$ be the Moore-Penrose pseudoinverse of $\boldsymbol{L}$. Both $\boldsymbol{L}$ and $\boldsymbol{L}^{\dagger}$ have the same eigenvectors, while two corresponding eigenvalues are reciprocals of each other, except that 1 replaces the zero eigenvalue of $\boldsymbol{L}$. One can verify that

[TABLE]

Also $\boldsymbol{L}^{\dagger}$ enjoys the structural property of $\boldsymbol{L}^{\dagger}\boldsymbol{1}=\boldsymbol{0}$.

The Dirichlet Laplacian $\boldsymbol{L}_{\circ}:=\boldsymbol{B}_{\circ}\,\mathrm{diag}(\boldsymbol{\sigma})\,\boldsymbol{B}_{\circ}{\!\!\!{}^{\top}}$ is made from $\boldsymbol{L}$ by discarding the entries corresponding to a reference node $d$. Let ${\cal E}={\cal E}_{1}+{\cal E}_{\circ}$, where ${\cal E}_{\circ}$ is the set of edges with one end connected to node $d$. One can verify that $\forall\,\boldsymbol{x}\neq\boldsymbol{0}$,

[TABLE]

which means $\boldsymbol{L}_{\circ}$ is positive definite and so invertible. Further, $\boldsymbol{L}_{\circ}\boldsymbol{1}=\boldsymbol{y}$ is a nonzero vector with non-negative coordinates, where $y_{i}=0$ if $id\in{\cal E}_{1}$ and $y_{i}>0$ if $id\in{\cal E}_{\circ}$.

Sometimes, it is misunderstood that $\boldsymbol{L}_{\circ}$ carries the same eigenvalues as $\boldsymbol{L}$ but the zero, which is not true.