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  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 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 . (ii) HD minimizes average total queue congestion , 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 and the average quadratic routing cost and can be made to move along this boundary by changing a control parameter that compromises between the two objectives and (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 -weighted schedulers using -exponent of queue lengths . As another extension in packet switches,  introduced Projective Cone Schedulers (PCS) to allow scheduling with non diagonal weight assignments. The work in  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 . Replacing queue-length by packet-age,  introduced a delay-based BP policy. To improve BP delay performance,  proposed place-holders with Last-In-First-Out (LIFO) forwarding. Adaptive redundancy was used in  to reduce light traffic delay in intermittently connected mobile networks. Using graph embedding,  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 .
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 , 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 , original BP policy , 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 , very different indeed from what is called HD in this paper. The results on minimum network delay are extended to multiclass wireless networks in . 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 . 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 .
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.
We denote vectors by bold lowercase and matrices by bold capital letters. Bywe denote the vector of all zeros, by the vector of all ones, and by
the identity matrix. On arrays:and are taken entrywise; and express entrywise comparisons; and denotes the Schur product. For as a vector, denotes its transpose, its diagonal matrix expansion, its Euclidean norm, and . For as a set, denotes its cardinality. We use as the scalar indicator function, and as the vector indicator function that its entry takes the value 1 if , and 0 otherwise. By we denote the time derivative of . For a variable related to a directed edge from node to node , we use notations and interchangeably. We use and to denote the sets of nodes neighbor to node with respectively incoming links to and outgoing links from node .
Note: To keep continuity and enhance readability of the manuscript, proofs are all placed in Appendix.
Consider a uniclass wireless network that operates in slotted time with normalized slots . The network is described by a simple, directed connectivity graph with set of nodes and directed edges . New packets, all with the same destination at node , 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 and 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.
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 where if channel is included, and otherwise.
Given a connectivity graph , scheduling set is the collection of all maximal scheduling vectors.
With denoting expectation, the expected time average of a discrete-time stochastic process is defined as
A queuing network is stable if queue at each node and at each slot , denoted as , has a bounded time average expectation, viz., .
Given a wireless network, an arrival vector is stabilizable if there exists a routing policy that can make the network stable under .
For a link , its capacity , which is frequently called transmission rate in literature, counts the maximum number of packets the link can transmit at slot . The link actual-transmission , on the other hand, counts the number of packets genuinely sent over the link at slot . Each link is also endowed with a cost factor that represents the cost of transmitting one packet over the link at slot ; for example, , with as defined in , or a cost associated with greedy embedding .
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:
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 . 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:
The loss function, 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 . 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:
where 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 region, with requiring no knowledge of traffic and topology.
2.2 Back-Pressure (BP) Policy
At each slot , the original BP  observes queue backlogs at network layer and estimates channel capacities to make a routing decision as follows.
BP weighting: For every link find link queue-differential and weight the link with
BP scheduling: Find a scheduling vector such that
where ties are broken arbitrarily.
BP forwarding: Over each activated link with transmit packets at full capacity . If there is no enough packets at node , 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 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
where trades queue occupancy for routing penalty, while recovers the original BP.
The V-parameter BP yields a Dirichlet routing cost within of its minimum feasible value to the detriment of growing average delay of relative to that of original BP . 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 and the network, with two negative consequences: (i) The same value of 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 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 , HD policy observes link queue-differentials at network layer and estimates channel capacities and channel cost factors to make a routing decision as follows.
HD weighting: For every link first calculate the number of packets it would transmit if it were activated as
where if node is the final destination, i.e., , and otherwise. The Lagrange control parameter 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
HD scheduling: Find a scheduling vector, in the same way as BP, using the max-weight scheduling, such that
where ties are broken arbitrarily.
HD forwarding: Over each activated link transmit number of packets, viz.,
where represents the number of packets genuinely sent over link at slot .
It is critical to discriminate among actual link transmissions , link transmission predictions and link capacities . Also notice that 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 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 or .
Table 1 compares HD and V-parameter BP algorithms, which emphasizes the same algorithmic structure, computational complexity and overhead signaling.
(i) Since by assumption, we get for all . (ii) If , we get due to (6) and from (7); in this case, even if the link were scheduled by (8), still no packet would be transmitted over it. (iii) If , we get and since due to (6), the link weight (7) still remains positive. (iv) In light of and , the value of never exceeds the number of packets in the transmitting node .
For , HD policy reduces to the initial adiabatic-based HD policy proposed in , where the packet forwarding follows a thermally adiabatic, and so insulated, heat process on each link.
In a special case that all links are of the same capacity, i.e., , and all link queue-differentials remain less than it, i.e., , HD policy with and -weighted policy of  with turn to be equivalent. Packet switches are well suited to this special case. It was suggested in  that a smaller may lead to a lower network delay, with a non-proven conjecture that heavy traffic delay is minimized when . A discussion of this was given in  along with some counterexamples. Even if the conjecture were true, note that for a multihop routing problem, the requirement of 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 , HD emphasizes on actual number of transmittable packets , 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 , where for is simplified to . This contrasts with BP weighting which is linear in . The quadratic weight is central to HD key property (Th. 1) which is fundamental to other HD qualities.
H3: Varying the penalty factor 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 (Th. 2).
At , the average total queue , 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 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 , the average routing cost 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 with changing between 0 and 1 as a function of , and . 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 is decreasing, and so as to have a higher impact, by increasing 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 and V-parameter BP for . The performance region is restricted to the set of all 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 (packets/slot) and link 2 of time-varying capacity . Assume one packet to arrive for each user at every timeslot. It is then easy to verify that for , the given arrival goes beyond the network capacity region.
For , Fig. 2 compares the performance of HD with and original BP. The left side panel depicts timeslot evolution of for . The right side panel shows the steady-state average of total queue length as a function of . For , both HD and BP perform the same. For , however, average total queue length increases linearly in 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 . Assume that at every timeslot a single packet arrives at node 1 destined for node . Following , let us take .
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 , trying with different values of indicates the weakness of V-parameter BP in aptly supporting the arrival. This simplistically shows one of the impacts of entering link cost factor 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 has a noise intensity which is randomly assigned at first and keeps fixed during the simulation. For each link, we adopt Shannon capacity with as power transmission and as bandwidth. At every timeslot, two packets arrive at nodes 1, 2, 3 and 4, destined for node . The aim is to minimize total with , which implicitly minimizes total power consumption in the network. For simplicity, let us fix and 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 and for the original BP (). 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 and for V-parameter BP with . 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 and in V-parameter BP as a function of . 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 and , and in Fig. 5 for total power consumption at and . 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 . As before, let be the number of existing packets at node at slot . State variables of the system can then be represented by the following vector:
Note that is discarded from state variables.
We use subscript to denote reduced vectors or matrices obtained by discarding the entries corresponding to the destination node .
Let a stochastic process represent the number of exogenous packets arriving into node at slot . Discard and compose the vector of node arrivals as
Also compose the vector of link actual transmissions as
where, as before, represents the number of packets actually sent over link at slot .
Given a directed graph , let denote the node-edge incidence matrix, in which is 1 if node is the tail of directed edge , if is the head, and 0 otherwise.111 In combinatorial geometry, one can view graph as a 1-complex, where is its 1-incidence matrix that describes the correlation between all oriented 1-cells (edges) and 0-cells (nodes) in the complex. Then denotes a reduction of that discards the row related to the destination node , which is referred to as reduced incidence matrix. One can verify that is a node vector, in which the entry corresponding to node reads the net outflow as
Using the above notation, the -controlled, stochastic state dynamics of a uniclass queuing network is captured by
Note that the link capacities vary by channel states, while the link actual transmissions are assigned by a routing policy subject to . This difference explains why despite traditional notation in literature, there is no need for 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
where represents the vector composed of as defined in (6).
Theorem 1 (HD Key Property).
Consider a uniclass wireless network constraint by capacity, directionality and interference. At every timeslot and for all , HD policy maximizes the -controlled functional
Consider the long-term average of functional defined as
Next assumption is being used in the analytical proofs of HD properties, stating that the greedy maximization of 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 . It also implicitly means no overlapping among slot-based substructures of maximization problem.
Consider a uniclass wireless network constraint by capacity, directionality and interference. Given a combination of network topology and traffic rates, timeslot maximization of is an optimal substructure for global maximization of .
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 , 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 to at each timeslot. To maximize , however, a portion of traffic must be forwarded through other paths that connect node to the destination.
5 HD Throughput Optimality
Let the stochastic process represent channel states at slot , describing all uncontrollable factors that affect wireless link capacities and cost factors. We assume that evolves according to an ergodic stationary process and takes values in a finite set . Thus, by Birkhoff’s ergodic theorem, each state has a probability of
where . Then the expected link capacities and cost factors are obtained as
where and represent the vectors composed of link capacities and link cost factors , 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
along with the law of large numbers imply
meaning that the expectations converge to the long-term averages. Thus, a routing policy could estimate and by observing timeslot variables and 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 , a destination node , and an ergodic stationary channel state process .
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.
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 .
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 is in the network capacity region, i.e., stabilizable, if and only if there exists a set of link actual transmissions that satisfy
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 .
Link actual transmissions are not fixed, but depend on routing policy. Further, there could potentially exist infinite number of routing policies that meet (17) for any stabilizable . 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
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:
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 makes the following lemmas possible.
Given a connected uniclass wireless network, is pseudo positive definite in the sense that all of its eigenvalues are positive and
is pseudo positive definite in the sense that all of its eigenvalues are positive andfor any vector , with equality if and only if .
Given a connected uniclass wireless network, for any vector , the following identity holds:
Given a connected uniclass wireless network, there exists such a scalar that for any vectors , the following inequality holds:
To analyze the HD throughput optimality, consider the Lyapunov candidate
Though is indeed an energy function in light of Lem. 1, due to the nonsymmetric weighting matrix , 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 or for the case that all links are of the same cost factor, is simplified to a scaled identity matrix that leads to , which in turn reduces to the sum of squares of queue lengths – a familiar Lyapunov function in most of previous results in literature.
Let be the Lyapunov drift. Substituting for from (10) leads to
Let us drop timeslot variable for ease of notation. Applying Lem. 3 to the first line of the above drift equation yields
Taking conditional expectation from the latter given the current queue backlogs and knowing that the term has a zero lower bound lead to
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 is a function only of control parameter and link cost factors . Since arrivals are independent of both and , we get
At the same time, both and are independent of , so is , which means . On the other hand, since the network layer routing controller has no impact on arrivals, turns to be an independent system variable that is not influenced by anything, which implies . Putting these results together yields
Given the current queue backlogs , let be the link actual transmissions provided by HD policy. As compared to any alternative transmission decision , Th. 1 secures for all and at each slot . Considering this with the equality (19) of Lem. 2 implies
Taking conditional expectation given current queues yields
As one alternative transmission decision to be compared with the provided by HD policy, consider the case where is produced by a routing algorithm which makes independent, stationary and randomized transmission decisions at each slot based only on arrivals and link capacities and so independent of both queue backlogs and link cost factors . Let us fix for such an algorithm and refer to it as . Using equality and considering that is independent from and , we obtain
, 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 bothand have finite upper bounds, and (iii) is a pseudo positive definite matrix in the sense of Lem. 1 with finite entries (recall ). Thus, the expected value of is finite at each slot , and so there exists a finite positive scalar such that . Utilizing this in the Lyapunov drift inequality yields
In the wake of (23), the next theorem is proven by showing that is always negative for all . (Refer to the Appendix for the end of the proof.)
Theorem 2 (HD Throughput Optimality).
Over any uniclass wireless network, HD policy with any is throughput optimal, meaning that it guarantees network stability under all stabilizable arrival vectors.
6 HD Minimum Delay at
Pareto optimal performance of HD policy stands on two pillars: minimization of the average queue congestion with , and minimization of the average routing cost with . 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 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.
At and under the K-hop interference model, timeslot maximization of the functional in (12) is equivalent to timeslot maximization of
It is critical to understand that Lem. 4 does not claim about the same maximum values for functionals and , which is obviously not true, but about the same maximizing control action at each slot . Another point is that while at each timeslot, HD maximizes for all , it maximizes for only .
Consider a uniclass wireless network under an arrival rate that is stabilized by a routing policy, resulting in average queue occupancies and average link actual transmissions . Then the following identity holds:
where for two random variables
where for two random variablesand , and .
To gain an insight into this lemma, consider a constant arrival vector which makes the right-hand side of (25) vanished. In light of , 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).
We specify -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, -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  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 -class routing policies and under the K-hop interference model, HD with minimizes the average total queue congestion 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  and references therein.
7.1 Heat Equations on Manifolds
On a smooth manifold charted in local coordinates , consider as spatial distribution of temperature, as heat flux, and as scalar field of heat sources (with minus for sinks). The law of heat conservation entails
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:
To have a unique solution, besides time initial condition, one must prescribe conditions on a boundary .
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 is replaced by a 0-chain vector representing the discrete domain, the pointwise functions and are respectively replaced by 0-cochain vectors and (node variables), the line integral is replaced by 1-cochain vector (edge variable), and the thermal diffusivity is replaced by a vector of edge weights .
As a 1-complex, the graph structure is fully described by its node-edge incidence matrix . 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 as the heat sink, combinatorial analogue of the classical heat equations (26)–(28) are obtained as
Note that the boundary on the manifold collapses to the single node on the graph at the fixed zero temperature, which absorbs all the heat generated by the heat sources .
where as before, subscript denotes a reduced vector or matrix that discards the entries corresponding to the destination node . The linear operator is called the Dirichlet Laplacian with respect to the node , which is a symmetric and diagonally dominant matrix. Further, it can be shown that for any connected graph, 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 from equations by fixing as boundary condition. Then we get the reduced set of continuous-time heat equations on an uncapacitated directed graph as
We refer to as nonlinear Dirichlet Laplacian that acts on a directed graph and, unlike , is an operand-dependent operator that retains neither linearity nor symmetry.
For the first time, heat diffusion on directed graphs is formulated via a nonlinear Laplacian. This is in agreement with the recent work in  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 (35) 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.
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)–(35) 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 be a realization of a continuous-time stochastic process along a sample path . Define the scaled process for any . A deterministic function is a fluid limit if there exist a sequence and a sample path such that 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 .
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