Jonatha Anselmi
email: jonatha.anselmi@inria.fr

N S Walton
email: n.s.walton@uva.nl

Decentralized Proportional Load Balancing

Keywords: Multipath routing, decentralized load balancer, utility optimization, proportional fairness AMS subject classifications. 68M10, 68W15, 90B18

  L'archive ouverte pluridisciplinaire

1. Introduction. We focus on communication and computer networks whose performance (or behavior) at equilibrium can be modelled by the solution of an utility maximization problem of the form max r∈R n r U r Λ r n r subject to r:j∈r

Λ r ≤ C j , j ∈ J over Λ r ≥ 0, r ∈ R, (1.1)
where R denotes the set of network paths (or routes), J denotes the set of network congestible resources (or links, queues), C j denotes the service speed (or capacity) of link j, n r denotes the number of packets (or customers) using route r, U r is a strictly-concave function, and the decision variable Λ r , r ∈ R, is interpreted as the throughput given by the network to all packets moving on route r, i.e. the rate at which packets flow on route r. The value of function U r (•) is interpreted as the utility of a packet sent on route r. The constraints of such optimization can be interpreted as the capacity constraints or the stability region of the network and the form of the utility functions depends on which dynamics or communication protocols packets are forced to follow. Furthermore, within the time scale of interest, the numbers of packets circulating on each route n r , r ∈ R, is considered constant.

Knowing the throughput of packets on each route r, that is the unique optimizer of (1.1), one immediately obtains other performance metrics of interest such as the delay experienced by each packet on r, which follows by an application of Little's law [START_REF] Little | A proof of the queueing formula L = λW[END_REF]. In the following, we say that the performance (or behavior) of a network is known if the per-route throughputs of that network are known. We will also refer to (1.1) as the single-path utility optimization.

The utility maximization (1.1) is known to provide a general model able to approximate the performance achieved by a wide class of both communication and computer networks. In the networking community, it is well accepted that the throughputs achieved by packets in telecommunication networks that implement a form of congestion control, e.g. TCP (Transmission Control Protocol), can be approximated by the unique optimizer of an optimization of the form (1.1) [START_REF] Kelly | Rate control in communication networks: shadow prices, proportional fairness and stability[END_REF][START_REF] Mo | Fair end-to-end window-based congestion control[END_REF][START_REF] Yi | Stochastic network utility maximisation[END_REF]. Similarly, in computer science, the throughputs achieved by customers in data-centers are often modelled by the throughputs of a multiclass BCMP or Kelly queueing network [START_REF] Urgaonkar | An analytical model for multi-tier internet services and its applications[END_REF][START_REF] Anselmi | Heavy-traffic revenue maximization in parallel multiclass queues[END_REF], and it is known that the throughputs of such queueing network can be approximated again by the unique optimizer of an optimization of the form (1.1) [START_REF] Schweitzer | Bottleneck determination in networks of queues[END_REF][START_REF] Walton | Proportional fairness and its relationship with multi-class queueing networks[END_REF][START_REF] Anselmi | Closed queueing networks under congestion: Nonbottleneck independence and bottleneck convergence[END_REF].

Keeping (1.1) as a general framework to model the behavior of communication and computer networks, we are interested in the problem of injecting a form of decentralized load-balancing control for the packets. This can be done by allowing each network entry point to send packets among all of its available routes. Such approach is motivated by the natural fact that communication and computer network resources can send and receive data from multiple paths simultaneously, and that developing such control is currently a desirable feature of strong interest, as it aims at improving performance, fairness and robustness of the network by better exploiting its intrinsic path diversity. There is a vast literature on this subject and we point the reader to Section 1.1 for an overview of the related work. In several cases of practical interest, it is desirable to design a decentralized load-balancing control with the property that the resulting throughputs converge, at equilibrium, to a solution of an utility maximization problem of the form where S denotes the set of network source-destination pairs, w s = r∈s n r denotes the total number of packets using all routes belonging to s, and the remainder of the notation is as above. Thus, the load-balancer is distributed among all the sourcedestination pairs s, and on each s it aims at maximizing the utility of the total throughput among all the available routes. This is often referred to as resource pooling. Here the individual pathwise constraints of (1.1) are combined and throughput ia optimized over different source-destination pairs. There are several desirable properties in designing a decentralized load-balancer achieving the optimization (1.2); see [START_REF] Key | Path selection and multipath congestion control[END_REF][START_REF] Raiciu | Coupled congestion control for multipath transport protocols[END_REF]. For instance, the optimizers of (1.2) can be interpreted as a Nash or Wardrop equilibrium where the delays on the routes of a given source-destination are all equal and minimum. Due to its utilitarian objective, throughputs achieved by (1.2) are fair. We will refer to (1.2) as the multipath utility optimization.

1.1. Related work. During the last decades, decentralized load balancing has been a central theme for different research areas. The objective is to improve the performance, robustness and fairness of a network by better exploiting its intrinsic path diversity. In data-centers, for instance, web or application servers are usually replicated several times [START_REF] Urgaonkar | An analytical model for multi-tier internet services and its applications[END_REF], are often multi-homed [START_REF] Raiciu | Improving datacenter performance and robustness with multipath tcp[END_REF] and, due to their large sizes, the network structure has several entry points. In telecommunication networks, smartphones can simultaneously use 3G and WiFi [START_REF] Paasch | Exploring mobile/wifi handover with multipath tcp[END_REF], and end hosts of peer-to-peer net-works can usually choose among different download sources, e.g. as in BitTorrent, [START_REF] Hacker | The end-to-end performance effects of parallel tcp sockets on a lossy wide-area network[END_REF]. This recurrent network structure is ideal for developing decentralized load-balancing schemes.

In the networking community, researchers have recently implemented a multipath version of the TCP (Transmission Control Protocol) [START_REF] Wischik | The resource pooling principle[END_REF][START_REF] Ford | Tcp extensions for multipath operation with multiple addresses[END_REF], and a myriad of (decentralized) load-balancing congestion-control algorithms have been developed during the last decade; see, e.g., [START_REF] Kelly | Stability of end-to-end algorithms for joint routing and rate control[END_REF][START_REF] Han | Multi-path tcp: A joint congestion control and routing scheme to exploit path diversity on the internet[END_REF][START_REF] Wischik | Design, implementation and evaluation of congestion control for multipath tcp[END_REF][START_REF] Honda | Multipath congestion control for shared bottleneck[END_REF][START_REF] Wang | Optimal flow control and routing in multi-path networks[END_REF][START_REF] Key | Path selection and multipath congestion control[END_REF][START_REF] Cao | Delay-based congestion control for multipath tcp[END_REF][START_REF] Peng | Multipath tcp algorithms: theory and design[END_REF][START_REF] Khalili | Mptcp is not pareto-optimal: Performance issues and a possible solution[END_REF][START_REF] Le | Improving the performance of multipath congestion control over wireless networks[END_REF]. Here an objective is to allocate bandwidth that solves a version of the multipath utility optimization (1.2). Necessary and sufficient conditions for achieving the multipath utility optimization are developed in [START_REF] Peng | Multipath tcp algorithms: theory and design[END_REF]. Convergence to the solution of the multipath utility optimization is usually proven through models based on ordinary differential equations.

1.2. Our Contribution. On top of the general model for network performance (1.1), we propose a decentralized algorithm that balances the amount of packets associated to each network route. On each entry point of the network, our load balancer only uses local information and works as follows:

Every time a packet completes on a route: i) increase the number of packets on that route by one; ii) remove one packet from a route that is selected with probability proportional to the number of packets on that route. Underlying this algorithm, there is the heuristic argument that the dynamics of the number of packets on each route should mimic the dynamics of a Processor-Sharing (PS) queue1 [START_REF] Kleinrock | Time-shared systems: a theoretical treatment[END_REF] (from the standpoint of each network resource that implements such control). This analogy with the PS queue is intuitively driven by the following coupling argument. When our load balancer increases by one the number of packets in the route where the completion occurred, this corresponds to an arrival event at the processor-sharing queue. Consequently, a packet removal from its route is interpreted as a service event at the queue and, as in the processor-sharing queue, this event occurs at a rate that is proportional to the number of packets in each route.

Using differential equations, we develop a mathematical model describing the evolution of the number of packets on each route that is induced by our load balancer. As we describe in Section 3.2, our model assumes that between any two consecutive loadbalancing updates, i.e. steps i) and ii) above, performed on any source-destination s, the network achieved the (equilibrium) throughputs given by the single-path utility optimization (1.1). Though we argue that our qualitative results hold even without it, this time-scale separation argument helps us in reducing the complexity of the underlying mathematical analysis and more importantly allows for an interpretation of the resulting equations. In fact, we will exhibit a Lyapunov function L(t) for the proposed dynamical system satisfying the following property:

dL dt = - s∈S Γ s D Λ r Γ s r∈s n r w s r∈s - s∈S Γ s D n r w s r∈s Λ r Γ s r∈s ,
where D(•||•) denotes the Kulback-Liebler divergence or relative entropy. From this, we understand the effect of our processor-sharing heuristic: our load-balancer can be interpreted as to be equilibrating the relative entropy between the distribution of packets along each route and the distribution of throughputs along these routes. This entropy approach contrasts previous studies of multipath routing, which tend to focus on placing quadratic bound on Lyapunov functions [START_REF] Kelly | Stability of end-to-end algorithms for joint routing and rate control[END_REF][START_REF] Key | Path selection and multipath congestion control[END_REF] and/or applying Nyquist's stability criterion [START_REF] Han | Multi-path tcp: A joint congestion control and routing scheme to exploit path diversity on the internet[END_REF]. Under proportional fairness [START_REF] Kelly | Rate control in communication networks: shadow prices, proportional fairness and stability[END_REF], i.e. U r (Λ r /n r ) = log(Λ r /n r ), we prove in our main result that the throughputs achieved by the packets converge to the solutions of the multipath utility optimization (1.2). Our proof is essentially based on duality arguments connecting the structure of the Lyapunov function L(t) to the objective of the multipath utility optimization. Then, we adapt our approach to show that it is possible to design an other load-balancer able to achieve the multipath utility optimization under a general class of utility functions. This follows through an optimization decomposition argument that is similar to the one used in [START_REF] Kelly | Charging and rate control for elastic traffic[END_REF].

The remainder of the paper is organized as follows. Section 2 introduces the modeling framework. Section 3 presents our decentralized load-balancer, describes application areas and develops a mathematical model for the dynamics of network packets. Section 4 is devoted to our main results (Theorems 4.2 and 4.4), which are proven in Section 5. A number of technical lemmas are deferred to the Appendix for simplicity of exposition.

2. Multipath Utility Optimization. Let the set J , indexed by j, be the set of links. These are the network resources affected by congestion and, depending on the context, can be also interpreted as queues, Internet routers, or web-servers. Packets circulate in the network following given routes. A route is an ordered set of links, r = (j 1 , . . . , j kr ) with j 1 , . . . , j kr ∈ J . Thus, k r is the number of links on route r. The set R is the set of routes. A source-destination, s, is a set of routes, thus s ⊂ R. Each source-destination contains the set of routes available to a host. We think of hosts as the entry points of the network and can be also interpreted as dispatchers or end users. The set S gives the set of source-destinations. We assume that there are no source-destinations that share a common route, i.e., s ∩ s = ∅ when s = s . Thus, we can define s(r) ∈ S as the source-destination to which route r belongs. We assume that there is one host for each source-destination.

Let n r denote the number of packets on route r ∈ R. Correspondingly, n = (n r : r ∈ R) is the number of packets on each route. Let w s denote the number of packets on source-destination s ∈ S. Accordingly, w = (w s : s ∈ S) is the number of packets on each source-destination. Thus, for s ∈ S (2.1)

w s = r∈s n r .
The strictly positive vector C = (C j : j ∈ J ) gives the capacity (or service rate) of each link, i.e. the number of packets served per unit of time in a fully utilized link. Let Λ r denote the throughput (or flow rate or bandwidth) allocated by the network to the packets on route r ∈ R, and let Γ s denote the throughput on source-destination s ∈ S. Thus, for s ∈ S and j ∈ J

Γ s = r∈s Λ r , r:j∈r Λ r ≤ C j .
2.1. Proportional Fairness. When a host can distribute packets among multiple routes, one wishes to design a load-balancer maximizing the utility of the total throughput between each source-destination pair, and thus achieving the multi-path utility maximization (1.2). We will first assume that both the optimizations (1.1) and (1.2) have logarithmic utilities, i.e.

U r (Λ r /n r ) = log(Λ r /n r ), U s (Γ s /w s ) = log(Γ s /w s ), (2.2)
for r ∈ R and s ∈ S respectively. The throughputs achieved with this utility function are called proportionally fair. Several economic and game-theoretic interpretations hold in this case [START_REF] Kelly | Charging and rate control for elastic traffic[END_REF][START_REF] Kelly | Rate control in communication networks: shadow prices, proportional fairness and stability[END_REF][START_REF] Low | Optimization flow control, i: Basic algorithm and convergence[END_REF][START_REF] Srikant | The Mathematics of Internet Congestion Control[END_REF][START_REF] Lin | Utility maximization for communication networks with multipath routing[END_REF][START_REF] Yi | Stochastic network utility maximisation[END_REF]. Communication protocols achieving this optimization objective have been designed and implemented [START_REF] Low | Understanding tcp vegas: a duality model[END_REF][START_REF] Viswanath | Opportunistic beamforming using dumb antennas. Information Theory[END_REF][START_REF] Kushner | Convergence of proportional-fair sharing algorithms under general conditions[END_REF]. Furthermore, proportionally fair throughputs approximate the performance of closed multiclass BCMP or Kelly queueing networks [START_REF] Schweitzer | Approximate analysis of multiclass closed networks of queues[END_REF][START_REF] Walton | Proportional fairness and its relationship with multi-class queueing networks[END_REF][START_REF] Anselmi | Closed queueing networks under congestion: Nonbottleneck independence and bottleneck convergence[END_REF], which have been used to model the performance of data centers [START_REF] Urgaonkar | An analytical model for multi-tier internet services and its applications[END_REF].

Summarizing so far, our primary goal is to design a decentralized load-balancer ensuring that the resulting throughputs solve max Γ,Λ≥0 s∈S

w s log Γ s s.t. Γ s = r∈s Λ r , s ∈ S, r:j∈r Λ r ≤ C j , j ∈ J . (2.3a)
We refer to problem (2.3) as the multipath proportionally fair optimization.

We will achieve this optimization by considering a simple load-balancing mechanism which modifies the single-path proportionally fair optimization:

max Λ≥0 r∈R n r log Λ r s.t. r:j∈r Λ r ≤ C j , j ∈ J . (2.4a)
3. Processor Sharing Load Balancing. We develop a simple decentralized algorithm that balances the amount of packets among multiple routes of a network. Operated at each host, it sends incoming packets along a given route with a rate proportional to the number of packets in transfer on that route. This modus operandi is similar to how a Processor-Sharing (PS) queue serves different packet classes [START_REF] Kleinrock | Time-shared systems: a theoretical treatment[END_REF]. The proportion of service devoted to a class in a PS queue is proportional to the number of jobs in that class. Given this analogy, we refer to our algorithm as "Processor Sharing Load Balancer" (PSLB), which is as follows:

Every time a packet completes on a route: i) increase the number of packets on that route by one; ii) then, remove one packet from a route that is selected with probability proportional to the number of packets on that route.

It is clear that PSLB is easy to implement on each host, which only needs to store the amount of packet on each of its routes. Note that the probability of decreasing the number of packets on each route is proportional to the number of packets on each route, as in a PS queue.

Application areas.

There exist at least four scenarios where our algorithm can be implemented: multipath-TCP, data-centers, wireless networks and peer-to-peer networks. We discuss each of these. i) Multipath-TCP : Multipath-TCP (MPTCP) refers to a set of extensions of TCP (Transmission Control Protocol) that allow each connection to use multiple paths simultaneously. MPTCP allows smartphones to simultaneously use 3G and WiFi [START_REF] Paasch | Exploring mobile/wifi handover with multipath tcp[END_REF]; it reacts to link failures on a time-scale that is much quicker than that obtained by updating the routing tables of network routers; it increases resource pooling -that is where multiple network resources can collectively be seen as a single pooled resource [START_REF] Laws | Resource pooling in queueing networks with dynamic routing[END_REF][START_REF] Gibbens | Dynamic alternative routing. Routing in communications networks[END_REF]. Due to these potential benefits, the Internet Engineering Task Force (IETF) MPTCP working group has recently released three implementations of MPTCP [START_REF] Raiciu | How hard can it be? designing and implementing a deployable multipath tcp[END_REF]. The first large-scale commercial deployment of MPTCP appeared on September 18th, 2013, when Apple released iOS 7 [START_REF] Bonaventure | Apple seems to also believe in multipath tcp[END_REF]. In this setting, our approach aims at balancing the magnitude of the congestion window associated to each network route available to host s given the total congestion window w s .

ii) Data Centers: In a data-center, packets are interpreted as http requests, links as Internet servers and each source-destination as the set of routes available to execute a web-application. In practice, end-users submit an http request, read on a browser the result that has been produced by the data-center, and then submit again an other request. Http requests can be assumed to come from a constant number of end-users, and the network performance can be modelled by a closed multiclass BCMP/Kelly queueing network; see, e.g., [START_REF] Urgaonkar | An analytical model for multi-tier internet services and its applications[END_REF][START_REF] Anselmi | Heavy-traffic revenue maximization in parallel multiclass queues[END_REF] and references therein. Given their large size, datacenters usually admit multiple entry points, called dispatchers, which can implement a load-balancing mechanism. By doing so and assuming that throughputs without load-balancing can be modelled by an optimization of the form (1.1) (as in, e.g., [START_REF] Anselmi | Heavy-traffic revenue maximization in parallel multiclass queues[END_REF]), we show that our load-balancer will make throughputs converge to a solution of an optimization of the form (1.2).

iii) Wireless Networks. In wireless systems, there are reasons to balance load amongst multiple resources. First, hosts may have access to multiple frequencies in the wireless spectrum and these are usually shared between users. Thus, load balancing is highly relevant in this dynamic spectrum access setting. Further, a wireless host may have access to multiple transmitters whose signal quality will, amongst other factors, depend on the distance between host and the transmitter. This access point selection problem represents a further model where load balancing is required. For an analogous proportional fair load-balancing system and discussion on these systems, see [START_REF] Yun | Distributed load balancing in heterogeneous systems[END_REF]. iv) Peer-to-peer networks. Although we refer to transfers between source destination pairs, it is clear from the above examples that a single-source is not required in our modeling setting. The case of peer-to-peer networks gives a further example. Here, a peer may download elements of a file from multiple peers. By exploiting the diversity of peers who share a specific file, users can receive high and fair download rates. Once again, the traffic must be balanced between multiple routes and thus load-balancing techniques are of importance, see [START_REF] Karger | Simple efficient load balancing algorithms for peer-to-peer systems[END_REF].

Differential Equation Model.

We model the time evolution of the number of packets induced by PSLB through differential equations. Between two consecutive load-balancing updates, we assume that the throughputs allocated to packets by the network are close to Λ(n), i.e., the equilibrium rates achieved by the single-path optimization (1.1). We stress that this is possible because rates converge to such single-path optimization. In other words, we assume that our load-balancer operates at a time-scale that allows the single-path equilibrium rates to be achieved.

A packet completion on a route r ∈ s will result in a increase in the number of packets on that route. In view of the time-scale separation argument above, the rate allocated by the network to such completions is Λ r (n). Each completion on a source-destination results in a decrease in the congestion of a route with probability proportional to the number of packets on that route. These source-destination completions occur at rate Γ s (n) and the probability of decrease on route r is nr ws . Thus, we model the time evolution of the number of packets in each route obtained when PSLB is used by the following system of ordinary differential equations: for n r (t) ≥ 0

dn r dt = Λ r (n) -Γ s (n) n r w s , (3.1)
for all r ∈ s and for all s ∈ S, where n(0) ≥ 0 is such that r∈s n r (0) = w s for all s. We refer to any solution n(t) of (3.1) as a fluid solution.

When each source-destination contains a single route and the total number of packets is kept constant, our model can be compared to the fluid-limit of the flow-level model studied in [START_REF] Bonald | Impact of fairness on internet performance[END_REF], which was proposed to model the time evolution of connections in a communication network with singlepath congestion control.

Before continuing, we observe that a special case of practical interest is a network of parallel servers, which gives a bipartite network structure between hosts and servers. The applications of wireless networks and data centers discussed above are two examples of a parallel server system. In this case, the solution to the proportional fair optimization takes a particularly simple form. Namely, the throughput at route r is proportional to the number of packets (or jobs) that have been sent there, i.e.,

Λ r (n) = n r r :j∈r n r C j ,
where j ∈ r. In this case, the above differential equation takes a very simple form.

Main Convergence Results.

Using the differential equation model developed in previous section, we prove our main results (Theorems 4.2 and 4.4). In Theorem 4.2, we prove that the throughputs achieved by PSLB converge to a solution of the multipath proportional fair optimization (2.3). In Theorem 4.4, we show that these properties can be also achieved under a wide class of utility functions.

4.1. Multipath Proportional Fair Optimization. We will demonstrate that our load balancer is able to modify the single-path flows, which solve the proportionally fair optimization, into a multipath flow control, which solves the multipath proportional fair optimization. For this purpose, we use the following function .

Letting L * def = min n≥0 L(n) subject to r∈s n r = w s for all s, we show that L(n) -L * is a Lyapunov function for the dynamical system (3.1).

The crucial observation about this function is given in the following proposition.

Theorem 4.1. For each point of differentiability of a fluid solution n(t), the following holds

dL(n(t)) dt = - s∈S Γ s (n) D (Λ r (n)) r∈s Γ s (n) (n r ) r∈s w s + D (n r ) r∈s w s (Λ r (n)) r∈s Γ s (n) , (4.2a)
where, for two probability distributions p s = (p r : r ∈ s) and q s = (q r : r ∈ s), D(p s ||q s ) defines the relative entropy

(4.3) D (p s ||q s ) = r∈s p r log p r q r .
Note that for two discrete probability distributions p and q on the same support, it is well known, using Jensen's inequality, that D (p||q) ≥ 0 and that D (p||q) = 0 if and only if p = q. This immediately implies that L(n(t)) is non-increasing.

Remark 1. The symmetric relative entropy

D(p s ||q s ) + D(q s ||p s ),
as first introduced by Kullback and Leibler [START_REF] Kullback | On information and sufficiency[END_REF], is critical to understanding the convergence of single-path connections to the multi-path optimum. Through the symmetric relative entropy (4.2), PSLB attempts to equate the distribution of packets within the network and the distribution of throughputs along different routes. If these are unequal then the Lyapunov function will decrease further towards its minimum. Through duality arguments presented shortly, we see that the minimum of the Lyapunov function L(n) is achieved only at solutions which minimize the multipath proportional fair optimization.

As a consequence of the last theorem, we obtain our main result, which says that dynamics of (3.1) converge to the multipath proportional fair optimization (2.3).

Theorem 4.2. For any fluid solution for which L(n(t)) is absolutely continuous, the source-destination throughputs converge to Γ * the optimal solution to the multipath proportionally fair optimization (2.3),

Γ(n(t)) ---→ t→∞ Γ * .
A proof of Theorem 4.1 and Theorem 4.2 can be found in Section 5.1.

Multipath Utility Optimization.

Our prior analysis achieved the multipath proportional fair objective. However, in this section we show that it can be generalized to a wide class of utility functions. Through an optimization decomposition argument, it is known that other utility optimizations can be achieved by altering the weights associated with the proportional fair objective [START_REF] Kelly | Charging and rate control for elastic traffic[END_REF]. Now, we apply a similar decomposition in order to modify proportional fair single-path flows towards a general multipath utility optimization.

First, we introduce some additional notation. For each source destination pair, let U s : R + → R be a strictly-increasing continuously-differentiable utility function. We further assume that the function Γ s → Γ s U s (Γ s ) is continuous, strictly increasing and with range (0, ∞). We note that this assumption holds for all the weighted α-fair allocations with α < 1 [START_REF] Mo | Fair end-to-end window-based congestion control[END_REF]. Given this assumption, we define the function g s (w s ) = Γ s such that (4.4)

g -1 s (Γ s ) = Γ s U s (Γ s ),
which is well-defined since Γ s U s (Γ s ) is increasing. We interpret g s (w s ) as the throughput source-destination pair s would expect given that it has utility function U s and given that its current number of packets from source to destination is w s . From this interpretation, we now consider a load balancer who obeys the following dynamics

(4.5) dn r dt = Λ r (n) -g s (w s ) n r w s ,
where Λ r (n) is the solution to the single-path proportionally fair optimization problem (2.4). In view of the fluid model (3.1), we note that these rates are a form of processorsharing load balancing: once again, we increment the number of packets for each packet arrival and decrement it proportionally to the number of packets present in each of the available routes. The only difference to our prior analysis is that the number of packets decreases at rate g s (w s ) rather than Γ s , so the number of packets between each source-destination pair is no longer conserved. Summing over r ∈ s, we see that w s must vary dynamically according to the equations

dw s dt = Γ s (n) -g s (w s ), where Γ s (n) = r∈s Λ s (n), w s = r∈s n r .
We want to control the number of packets between each source destination pair so that the logarithm in the proportionally fair optimization is replaced by the more general utility function U s . With this in mind, we consider the function

(4.6) G s (w s ) = max Γs≥0 {w s log Γ s -U s (Γ s )} .
With the addition of function G s and motivated by the structure of (4.1), we define the following function

(4.7) H(n) = r∈R n r log n r w s Λ r (n) + s∈S G s (w s ).
Letting H * def = min n≥0 H(n) subject to r∈s n r = w s for all s, we show that H(n) -H * is a Lyapunov function for the dynamical system (4.5).

Once again, the function H balances the relative entropy between different routes of the network. This is shown in the following proposition, which is proven in Section 5.2.

Theorem 4.3. For each point of differentiability of a fluid solution for (4.5), n(t), the following holds

dH(n(t)) dt = - s∈S Γ s D (Λ r ) r∈s Γ s (n r ) r∈s w s - s∈S g s (w s )D (n r ) r∈s w s (Λ r ) r∈s Γ s - s∈S Γ s -g s (w s ) log Γ s -log g s (w s ) . (4.8)
Using that the logarithm is increasing, the term (4.8) is negative. Furthermore, since D(•||•) ≥ 0, this means that the derivative of the entropy function H(n(t)) is negative.

As a consequence of Theorem 4.3 and adapting the arguments in the proof of Theorem 4.2, we can prove that the single-path flows converge to solutions of the multipath utility optimization (1.2).

Theorem 4.4. For any fluid solution for which H(n(t)) is absolutely continuous and w(t) is bounded from zero, the source-destination throughputs converge to Γ * the optimal solution to the multipath utility optimization

(1.2), Γ(n(t)) ---→ t→∞ Γ * .
A proof of Theorem 4.3 and 4.4 can be found in Section 5.2.

Proofs of Main Results

. We present proofs for our main results: Theorem 4.2 and Theorem 4.4. For simplicity, we often drop the dependency on n when referring to quantities Λ r (n) and Γ s (n).

5.1. Proof of Theorem 4.1 and Theorem 4.2. We require the following lemma which follows by somewhat standard calculus arguments, for instance see [START_REF] Walton | Concave switching in single and multihop networks[END_REF]Lemma 2] for a proof.

Lemma 5.1.

If n r > 0 then ∂L(n) ∂n r = log n r w s Λ r (n) + 1.
With this, we can prove Theorem 4.1.

Proof. [Proof of Theorem 4.1] For the following sequence of equalities, we use the shorthand that n = n(t), L = L(n(t)) and Γ = Γ(n(t)):

dL dt = r∈R dn r dt ∂L ∂n r = s∈S r∈s nr>0 Λ r - n r w s Γ s log n r w s Λ r + 1 = s∈S r∈s Λ r - n r w s Γ s log n r w s Λ r (5.1a) = - s∈S Γ s r∈s Λ r Γ s log w s Λ r n r Γ s - s∈S Γ s r∈s n r w s log n r Γ s w s Λ r (5.1b) = - s∈S Γ s D (Λ r ) r∈s Γ s (n r ) r∈s w s - s∈S Γ s D (n r ) r∈s w s (Λ r ) r∈s Γ s (5.1c)
In the first two equalities above, we apply the chain rule and the definition of the derivatives of n and L, from fluid equations (3.1) and Lemma 5.1. In the equality (5.1a), we cancel terms and observe that we can exclude terms where n r = 0 because dnr dt = 0 and Λ r (n) = 0 must hold at times when n r = 0. Then, in equalities (5.1b) and (5.1c) we add and subtract terms Γ s log Γ s and observe that the terms summed are relative entropy terms.

To prove Theorem 4.2, we will need the following technical lemma, which is straightforward to verify and gives first order properties for our Lyapunov function L.

Lemma 5.2. L(n) is continuous and convex and thus N * , the set of the optimal solutions to (5.1), is compact and convex.

A proof of Lemma 5.2 can be found in the Appendix. The following Lemma bounds terms Γ s (n).

Lemma 5.3. a) The function L(n) can be characterized through its dual as follows

L(n) = max q≥0 r∈R n r log   j∈r q j   - j∈J q j C j + r∈R n r - s∈S w s log w s . (5.

2)

b) There exists a constant Γ min > 0 not dependent on n such that for all n and each s ∈ S we have Γ s (n) ≥ Γ min . A proof of Lemma 5.3 can be found in the Appendix.

The following lemma finds a dual description for the Lyapunov function L and also equates the minimized Lyapunov function with the multipath proportionally fair optimization.

Lemma 5.4. The minimum of Lyapunov function L(n) relates to the multi-path proportionally fair optimization as follows

min n≥0 L(n) subject to r∈s n r = w s , s ∈ S = -max Γ,Λ≥0 s∈S w s log Γ s subject to Γ s = r∈s Λ r , s ∈ S.
Moreover, L(n) is minimal if and only if

(5.3) n r w s = Λ r (n) Γ s (n)
for each r ∈ R (with r ∈ s) such that n r > 0.

Proof. We have the following equalities

min n≥0 L(n) subject to r∈s n r = w s , s ∈ S = min n≥0 min Λ≥0 s∈S r∈s n r log n r Λ r w s s.t. r:j∈r Λ r ≤ C j , r∈s n r = w s , j ∈ J , s ∈ S = min Λ≥0 min n≥0 s∈S r∈s n r log n r Λ r w s s.t. r∈s n r = w s , r:j∈r Λ r ≤ C j , j ∈ J , s ∈ S = min Λ≥0 min n≥0 s∈S w s r∈s n r w s log n r w s Γ s Λ r - s∈S w s log Γ s subject to r∈s n r = w s , r:j∈r Λ r ≤ C j , r∈s Λ r = Γ s = min Λ≥0 - s∈S w s log Γ s s.t. r:j∈r Λ r ≤ C j , r∈s Λ r = Γ s
The first equality follows by the definition of L(n); in the second equality, we interchange the order of minimization; in the third, we add and substract the multipath proportionally-fair objective; in the final equality, we minimize over n and observe that the relative entropy term in the objective is minimized to zero when the vector n is proportional to Λ. Now we show that (5.3) is a condition for optimality of L(n). Consider the function

(5.4) L(n, Λ) = s∈S r∈s n r log n r Λ r w s .
It is clear that L(n, Λ) is convex. The minimum of L(n, Λ) over Λ (satisfying the network capacity constraints r:j∈r Λ r ≤ C j , j ∈ J ) is achieved at Λ(n) and the minimum attained at L(n). The minimum of L(n, Λ) over n (satisfying the packet constraints r∈s n r = w s , is achieved when n r ∝ Λ r . Since we are minimizing a convex function over a product space, its minimum is the intersection of the marginal minimality conditions, i.e. L(n, Λ) and thus L(n) is minimized iff Λ r = Λ r (n), for r ∈ R and, for all r ∈ s,

Λ r (n) Γ s (n) = n r w s .
The proof of Theorem 4.2 is now fairly straight-forward.

Proof. [Proof of Theorem 4.2] We let L * be the minimal value of L(n) over n with r∈s n r = w s , s ∈ S. The last lemma states that (5.5)

L(n) = L * iff n r w s = Λ r (n) Γ s (n) .
From Proposition 4.1, we have that

dL dt = - s∈S Γ s D (Λ r ) r∈s Γ s (n r ) r∈s w s - s∈S Γ s D (n r ) r∈s w s (Λ r ) r∈s Γ s .
Thus, as the relative entropy D(p||q) = 0 iff p = q, we can combine the last two displays to give that

(5.6) dL dt = 0 iff Λ r (n) Γ s (n) = n r w s , r ∈ s, s ∈ S iff L(n) = L * .
So L(n(t)) is decreasing unless n(t) achieves the minimum of L(n). Suppose that L(n(t)) decreases to some value L (which may or may not equal L * ). There must be a sequence of times {t k } k∈N such that dL(t k ) dt → 0 as t k → ∞. By Lemma 5.3b, it can not be that Γ s (n(t k )) → 0 as k → ∞, so for this sequence of times it must be that

(5.7) D (n r (t k )) r∈s w s (Λ r (t k )) r∈s Γ s (t k ) ----→ k→∞ 0.
Thus, by continuity of the relative entropy at 0, any limit point, n * , Λ * , Γ * of the sequences n(t k ), Λ(n(t k )), Γ(n(t k )) must satisfy

Λ * r Γ * s = n * r w s ,
for all r ∈ s, s ∈ S. Thus the minimum of L(n(t)) must be the optimal value L = L * . So L(n(t)) converges to L * , we now show that Γ(n(t)) converges to Γ * . Observe

L(n(t)) -L * = s∈S r∈s n r (t) log n r (t)Γ s (n(t)) w s Λ r (n(t)) + s∈S w s log Γ * s - s∈S w s log(Γ s (n(t))) ≥ s∈S w s log Γ * s - s∈S w s log(Γ s (n(t))) ≥0,
where in the first inequality, above, we use the positivity of the relative entropy. So since L(n(t)) -L * converges to 0, and since Γ → s w s log Γ s is continuous and Γ s (n(t)) is feasible for each t, it must be that, for s ∈ S,

(5.8) Γ s (n(t)) ---→ t→∞ Γ * s ,
as required.

Proof of Theorem 4.3 and Theorem 4.4.

Next lemma is a consequence of Fenchel duality and relates the function G s given in (4.6) with the utility function U s and g s in (4.4). The lemma also gives the derivative of our Lyapunov function.

Lemma 5.5. a) The function G s (w s ) is convex. The maximization for G s (w s ) given by (4.6) is optimized by g s (w s ). The function G s is differentiable with derivative (5.9)

G s (w s ) = log g s (w s ).

Further,

(5.10) -U s (Γ s ) = min ws≥0 {G s (w s ) -w s log Γ s } .
b) The partial derivatives of H are given by the expression

(
- s∈S Γ s log Γ s + s∈S g s (w s ) log Γ s + s∈S Γ s G s (w s ) - s∈S g s (w s )G s (w s ) = - s∈S Γ s D (Λ r ) r∈s Γ s (n r ) r∈s w s - s∈S g s (w s )D (n r ) r∈s w s (Λ r ) r∈s Γ s (5.12b) - s∈S Γ s -g s (w s ) log Γ s -log g s (w s ) .
In the second equality, we used Lemma 5.5.b. In (5.12a), we have added terms Γ s log Γ s and multiplied/divided by Γ s to express summands in terms of the relative entropy in (5.12b). In (5.12b), we have also used (5.9).

We let H * be the minimal value of H(n) over n with r∈s n r = w s , s ∈ S. Next lemma extends Lemma 5.4.

Lemma 5.6. The minimum of H(n) relates to the multi-path utility optimization as follows min n≥0,w≥0

H(n) subject to r∈s n r = w s , s ∈ S = -max Γ,Λ≥0 s∈S U s (Γ s ) subject to Γ s = r∈s Λ r , s ∈ S, r:j∈r Λ r ≤ C j , j ∈ J .
Moreover, H(n) is minimal if and only if

(5.13) n r w s = Λ r (n) Γ s (n)
for each r ∈ R (with r ∈ s) such that n r > 0 and

(5.

= Λ r (n) Γ s (n) , n r > 0.
So H(n(t)) is decreasing unless n(t) achieves the minimum of H(n). Suppose that H(n(t)) decreases to some value H (which may or may not equal H * ). There must be a sequence of times {t k } k∈N such that dH(t k ) dt → 0 as t k → ∞. Since w(t) is bounded from zero, it can not be that Γ s (n(t k )) → 0 as k → ∞, so for this sequence of times it must be that

(5.18) D (n r (t k )) r∈s w s (Λ r (t k )) r∈s Γ s (t k ) ----→ k→∞ 0. and (5.19) Γ s (t k ) -g s (w s (t k )) log Γ s (t k ) -log g s (w s (t k )) ----→ k→∞ 0.
Thus, by continuity of relative entropy at 0, any limit point, n * , Λ * , Γ * of the sequences n(t k ), Λ(n(t k )), Γ(n(t k )) must satisfy

Λ * r Γ * s = n * r w s and g s (w s ) = Γ s (n)
for all r ∈ s, s ∈ S.

Thus the decreasing process H(t) must be decrease to the minimum it value H * . Finally, to consider the convergence of the processes n(t), w(t), Λ(t) and Γ(t), we note that

H(t) = H(n(t), w(t), Λ(t)).
for the continuous convex function H given by (7.10).

H(n, w, Λ) = s∈S r∈s n r log n r w s Λ r + s∈S G s (w s ).
It is clear that H has compact level sets and so must converge n(t), w(t) and Λ(t) must converge to the arguments minimizing H over r∈s n r = w s , for s ∈ S, r j Λ r ≤ C j , j ∈ J . As a consequence, Γ(t) must converge to the unique solution to the the optimization (1.2).

So H(n(t)) converges to H * , we now show that Γ(n(t)) converges to Γ * . Observe

H(n(t)) -H * = s∈S r∈s n r (t) log n r (t)Γ s (n(t)) w s Λ r (n(t)) + s∈S U s (Γ * s ) - s∈S w s log(Γ s (n(t))) ≥ s∈S U s (Γ * s ) - s∈S w s log(Γ s (n(t))) ≥0.
So since H(n(t)) -H * converges to 0, and since Γ → s w s log Γ s is continuous, it must be that, for s ∈ S,

(5.20) Γ s (n(t)) ---→ t→∞ Γ * s ,
as required.

7. Proof of Additional Lemmas. Proof. [Proof of Lemma 5.2] To prove that L(n) is convex, we use the dual formulation proven in Lemma 5.4. Let D(n, q) be the objective of the dual function (5.2). Observe D(n, q) is linear in n. For probabilities p = (p 0 , p 1 ) with p 0 + p 1 = 1 and for two vectors n 0 , n 1 ∈ R R + , we define n p = p 0 n 0 + p 1 n 1 . We have

L(n p ) = max q≥0 D(n p , q) = max q≥0 p 0 D(n 0 , q) + p 1 D(n 1 , q) ≤ p 0 max q≥0 D(n 0 , q) + p 1 max q≥0 D(n 1 , q) = p 0 L(n 0 ) + p 1 L(n 1 ).
Thus, we see that L(n) is a convex function.

We consider the continuity of L(n). Firstly, it is proven that n → Λ r (n) is a continuous function for n r > 0, see [START_REF] Kelly | Fluid model for a network operating under a fair bandwidthsharing policy[END_REF]Lemma A.3]. Thus each summand, r, in the above expression is continuous for n r > 0. However, the case where n r = 0 remains. Since Λ r (n) ≤ C max and by bound (7.7) proven in Lemma 5.4,

(7.1) n r log n r C max w s ≤ n r log n r Λ r (n)w s ≤ n r log w max C min w s .
Thus, lim nr→0 n r log nr Λr(n)ws = 0, which means that the summands of L(n) are also continuous at n r = 0. Therefore, L(n) is continuous because it is the sum of continuous functions.

Finally, we note that the set of feasible values of n belongs to a compact set, so the set of optimal solutions is bounded and since L(n) is convex continuous the set of optimal solutions is also closed and convex.

Proof. [Proof of Lemma 5.3] a) The function L(n) is the solution to the concave optimization

L(n) = min - r n r log Λ r w s n r s.t. j∈r Λ r ≤ C j , j ∈ J over Λ r ≥ 0, r ∈ R.
We calculate the dual of this optimization. After adding slack variables z j and Lagrange multipliers q j , the Lagrangian of this optimization is

(7.3) - r n r log Λ r w s n r - j q j C j - r:j∈r Λ r -z j .
We minimize such Lagrangian. As the slack variables are positive, it is clear that a finite solution q j must be positive for each j ∈ J , and by complementary slackness q j z j = 0. Differentiating with respect to Λ r , we see that the Lagrangian of this optimization is minimized when

(7.4) n r Λ r = j∈r q j .
Substituting this optimality condition back in to the Lagrangian (7.3), we see that the dual of the optimization is

max q≥0 r n r log j∈r q j - j q j C j + r n r - s w s log w s .
The primal problem, (7.2), satisfies Slater's condition and so is Strong Lagrangian. Hence, (5.2) holds.

b) We first place some bounds on the dual formulation proved above. The partial derivative of the objective, with respect to q j for any j ∈ r, is (7.5) n r j∈r q j -C j .

For an optimal choice of q this quantity must be negative or equal to zero if q j > 0. However, if we take q such that (7.6) r:j∈r qj > w max C min then, (7.5) is negative. Thus, it must be that

(7.7) n r Λ r (n) = j∈r q j ≤ w max C min
The equality holds by the optimality condition (7.6) holding. for the optimal q. The lowerbound is achieved by noting that Λ r (n) ≤ C max for all n and r. For each s, there must exists some route such that n r ≥ ws |s| and certainly Γ s (n) ≥ Λ r (n) for each r ∈ s. Applying this to the above upper-bound derived, we have that

w s |s|Γ s (n) ≤ n r Λ r (n) ≤ w max C min .
Rearranging gives, Γ s (n) ≥ wmax|s| Cminws =: Γ min , as required. Proof. [Proof of Lemma 5.5] a) By the definition of G s (w s ) (see (4.6)), we have G s (w s ) = max Γs≥0 {w s log Γ s -U s (Γ s )} = max γs∈R {w s γ s -U s (e γs )} (7.8) Thus, G s (w s ) is the Legendre-Fenchel transform of the convex function U s (e γs ) and thus it is convex. The objective w s γ s -U s (e γs ) is concave provided the function U s (e γs ) is convex, which is equivalent to our assumption that Γ s U s (Γ s ) is increasing. Differentiating this objective, we see that optimality is achieved by the Γ s satisfying w s = Γ s U s (Γ s ). In other words, the optimum for G s (w s ) is achieved by Γ s = g s (w s ). Next, we see that we can bound the derivative of G s from the right as follows G s (w s + h) -G s (w s ) = (w s + h) log g s (w s + h) -U s (g s (w s + h)) -w s log(g s (w s )) + U s (g s (w s )) ≥ (w s + h) log(g s (w s )) -U s (g s (w s )) -w s log(g s (w s )) + U s (g s (w s )) = h log(g s (w s )).

By a symmetric argument we have h log(g s (w s + h)) ≥ G s (w s + h) -G s (w s ). Dividing by h and by the continuity of the function g s , we have from that these bounds give (7.9)

G s (w s ) = log g s (w s ).

This gives the derivative of G s . Next, recall that G s (w s ) is the Legendre-Fenchel transform of the convex function U s (e γs ), see (7.8). Thus, Legendre-Fenchel duality holds and we can see that U s (e γs ) = max ws {w s γ s -G s (-w s )} .

And so -U s (Γ s ) = min ws {G s (w s ) -w s log Γ s }. This proves part a). b) We proceed in a similar manner to Lemma 5.1. We have Differentiating the first two terms using standard arguments, differentiating the third term as in Lemma 5.1, and using (7.9), we get (5.11).

H(n) =
Proof. [Proof of Lemma 5.6] We minimize H(n). For this we consider the following sequence of equalities min In the first equality above, we introduce the optimization taken over Λ; in the second equality, we reorder the minimization and add and subtract relevant terms; in the third, we minimize the above relative entropy term to zero; in the forth we take the Fenchel transform of G s (w s ) to obtain U s and the required optimization. The function is convex over (Λ r : r ∈ R) and (n r : r ∈ R) and (w s : s ∈ S) satisfying r∈s n r = w s , for s ∈ S. Minimizing over Λ satisfying network capacity constraints r j Λ r ≤ C j , j ∈ J , we arrive at the proportionally fair optimization, Λ(n), i.e. H(n) = H * (n, w, Λ(n)). Minimizing over n r (as in the third equality above) we obtain the condition (7.11)

n r w s = Λ r Γ s
for each r ∈ R (with r ∈ s) such that n r > 0. Further minimizing over w s we obtain the condition (7.12) g s (w s ) = Γ s (n).

Together we see that the minimality condition for H and thus for H are (5.13) and (5.14).

  maximize s∈S w s U s (Γ s /w s ) Λ r ≥ 0, s ∈ S, r ∈ R, (1.2d)

  n r w s Λ r (n)

r∈R n r

  log n r -r∈R n r log Λ r (n) -s∈S w s log w s + s∈S G s (w s ).

  s w s Λ r + s∈S {G s (w s ) -w s log Γ s } = min Λ,Γ≥0: r∈s Λs=Γs r j Λr≤Cj min w≥0 s∈S {G s (w s ) -w s log Γ s } = min Λ,Γ≥0 s∈S U s (Γ s ) subject to r j Λ r ≤ C j , j ∈ J , r∈s Λ s = Γ s , s ∈ S.

  [START_REF] Kelly | Stability of end-to-end algorithms for joint routing and rate control[END_REF] g s (w s ) = Γ s (n)for each s ∈ S.A proof of Lemma 5.6 can be found in the Appendix. With Lemma 5.6 we can prove Theorem 4.4.Proof. [Proof ofTheorem 4.4] We let H * be the minimal value of H(n) over n. -g s (w s ) log Γ s -log g s (w s ) .

	Last lemma says that
	(5.15)	H(n) = H * iff g s (w s ) = Γ s (n) and	n r w s	=	Λ r (n) Γ s (n)	, n r > 0.
	From Theorem 4.3 we have that
	(5.16a)
	dH(n(t)) dt	= -	s∈S	Γ s D	(Λ r ) r∈s Γ s	(n r ) r∈s w s	-	s∈S	g s (w s )D	(n r ) r∈s w s	(Λ r ) r∈s Γ s
	(5.16b)		-
	It is clear that dH dt is negative and combining the last two displays, (5.15) and (5.16),
	we have that
	(5.17)	dH dt	= 0	iff	g s (w s ) = Γ s (n) and		n r w s

s∈S

Γ s

In a PS queue, all jobs are processed simultaneously and each of them takes 1/n of the available capacity, provided that n denotes the number of jobs in the queue.

6. Acknowledgments. The work of Jonatha Anselmi has been supported by grant MTM2010-17405 (Ministerio de Ciencia e Innovación, Spain) and grant SA-2012/00331 of the Department of Industry, Innovation, Trade and Tourism (Basque Government), which also sponsored a one-month visit of Neil Walton to the Basque Center for Applied Mathematics (BCAM). The work of Neil Walton is fund by the VENI research programme, which is financed by the Netherlands Organisation for Scientific Research (NWO).