A centre–free approach for resource allocation with lower bounds

ABSTRACT Since complexity and scale of systems are continuously increasing, there is a growing interest in developing distributed algorithms that are capable to address information constraints, specially for solving optimisation and decision-making problems. In this paper, we propose a novel method to solve distributed resource allocation problems that include lower bound constraints. The optimisation process is carried out by a set of agents that use a communication network to coordinate their decisions. Convergence and optimality of the method are guaranteed under some mild assumptions related to the convexity of the problem and the connectivity of the underlying graph. Finally, we compare our approach with other techniques reported in the literature, and we present some engineering applications.


Introduction
The increasing scale and complexity of systems have motivated the development of distributed methods to deal with situations where optimisation and decisionmaking are required. An important issue within this field is the optimal resource allocation over networks of agents, a problem that is closely related to network utility maximisation (NUM) problems (Palomar & Chiang, 2006;Tan, Zhu, Ge, & Xiong, 2015). Resource allocation arises when there is a limited amount of a certain resource (e.g. electric power, computing capacity or execution time), and it is necessary to establish an optimal distribution policy between some entities (e.g. loads, processors or controllers) that are connected by a communication network. This kind of problems has a large number of applications in economics (Ayesta, Erausquin, Ferreira, & Jacko, 2016;Conrad, 2010), smart energy systems (Hansen, Roche, Suryanarayanan, Maciejewski, & Siegel, 2015;Pantoja & Quijano, 2012), cloud computing (Pietrabissa et al., 2016;Pillai & Rao, 2016), and communications (H. Lee, K.J. Lee, Kim, Clerckx, & I. Lee, 2016;Tan et al., 2015).
Although there exists an extensive literature regarding distributed methods for solving resource allocation problems, this field still attracts considerable research attention (Cherukuri & Cortés, 2015;Obando, Pantoja, & Quijano, 2014;Pantoja, Quijano, & Passino, 2014;Poveda & Quijano, 2015;Ramirez-Llanos & Martinez, 2015;Tan, Yang, & Xu, 2013). Most of the solution methods are CONTACT Germán Obando ge-oband@uniandes.edu.co based on multi-agent systems (e.g. a survey that deals with the general class of NUM problems can be found in Palomar & Chiang, 2006), where the agents make decisions based on local information in order to obtain a desirable global behaviour. Appropriate coordination of agents is crucial because it avoids converging to suboptimal solutions. In order to ensure this coordination, a large number of methods require either the inclusion of a centralised agent or the use of restrictive information structures (as it is pointed out in Mosk-Aoyama, Roughgarden, & Shah, 2010). For instance, in classic decomposition techniques (Bemporad, Heemels, & Johansson, 2010;Boyd, Parikh, Chu, Peleato, & Eckstein, 2010;Palomar & Chiang, 2006), the Lagrange multiplier related to the 'price' of the resource is centrally adjusted to reach the optimum. By contrast, other methods are fully decentralised (e.g. Barreiro-Gomez, Obando, & Quijano, 2016;Xiao & Boyd, 2006;Zhu & Martinez, 2012). These methods exploit the communication capabilities of the agents to coordinate their decisions based on the information received from their neighbours. Fully decentralised methodologies have important advantages, among which we highlight the increase of the autonomy and resilience of the whole system since the dependence on a central authority is avoided.
In this paper, we propose a distributed resource allocation algorithm that does not require a central coordinator. An important characteristic of our method is the capability of handling lower bounds on the decision variables. This feature is crucial in a large number of practical applications, e.g. in Conrad (2010), Pantoja and Quijano (2012), and Lee et al. (2016), where it is required to capture the non-negativity of the resource allocated to each entity. We use a Lyapunov-based analysis in order to prove that the proposed algorithm asymptotically converges to the optimal solution under some mild assumptions related to the convexity of the cost function, and the connectivity of the graph that represents the communication topology. In order to illustrate our theoretical results, we perform some simulations and compare our method with other techniques reported in the literature. Finally, we present two engineering applications of the proposed algorithm. The first one seeks to improve the energy efficiency in large-scale air-conditioning systems. The second one is related to the distributed computation of the Euclidean projection onto a given set.
Our approach is based on a continuous time version of the centre-free algorithm presented in Xiao and Boyd (2006). The key difference is that the method in Xiao and Boyd (2006) does not allow the explicit inclusion of lower bounds on the decision variables, unless they are added by means of barrier functions (either logarithmic or exact; Cherukuri & Cortés, 2015). The problem of using barrier functions is that they can adversely affect the convergence time (in the case of using exact barrier functions) and the accuracy of the solution (in the case of using classic logarithmic barrier functions), especially for large-scale problems (Jensen & Bard, 2003). There are other methods that consider lower bound constraints in the problem formulation. For instance, Dominguez-Garcia, Cady, and Hadjicostis (2012) and Tan et al. (2013) have developed a decentralised technique based on broadcasting and consensus to optimally distribute a resource considering capacity constraints on each entity in the network. Nonetheless, compared to our algorithm, the approach in Dominguez-Garcia et al. (2012) and Tan et al. (2013) is only applicable to quadratic cost functions. On the other hand, Pantoja and Quijano (2012) propose a novel methodology based on population dynamics. The main drawback of this technique is that its performance is seriously degraded when the number of communication links decreases. We point out the fact that other distributed optimisation algorithms can be applied to solve resource allocation problems, as those presented in Nedic, Ozdaglar, and Parrilo (2010), Yi, Hong, and Liu (2015), and Johansson and Johansson (2009). Nevertheless, the underlying idea in these methods is different from the one used in our work, i.e. Nedic et al. (2010), Yi et al. (2015), and Johansson and Johansson (2009) use consensus steps to refine an estimation of the system state, while in our approach, consensus is used to equalise a quantity that depends on both the marginal cost perceived by each agent in the network and the Karush-Kuhn-Tucker (KKT) multiplier related to the corresponding resource's lower bound. In this regard, it is worth noting that the method studied in this paper requires less computational capability than the methods mentioned above. Finally, there are other techniques based on game theory and mechanism design (Kakhbod & Teneketzis, 2012;Sharma & Teneketzis, 2009) that decompose and solve resource allocation problems. Nonetheless, those techniques need that each agent broadcasts a variable to all the other agents, i.e. a communication topology given by a complete graph is required. In contrast, the method developed in this paper only uses a communication topology given by a connected graph, which generally requires lower infrastructure.
The remainder of this paper is organised as follows. Section 2 shows preliminary concepts related to graph theory. In Section 3, the resource allocation problem is stated. Then, in Section 4, we present our distributed algorithm and the main results on convergence and optimality. A comparison with other techniques reported in the literature is performed in Section 5. In Section 6, we describe two applications of the proposed method: (i) the optimal chiller loading problem in large-scale airconditioning systems, and (ii) the distributed computation of Euclidean projections. Finally, in Sections 7 and 8, arguments and conclusions of the developed work are presented.

Preliminaries
First, we describe the notation used throughout the paper and presents some preliminary results on graph theory that are used in the proofs of our main contributions.
In the multi-agent framework considered in this article, we use a graph to model the communication network that allows the agents to coordinate their decisions. A graph is mathematically represented by the pair G = (V, E ), where V = {1, . . . , n} is the set of nodes, and E ⊆ V × V is the set of edges connecting the nodes. G is also characterised by its adjacency matrix A = [a i j ]. The adjacency matrix A is an n × n non-negative matrix that satisfies: a ij = 1 if and only if (i, j) ∈ E, and a ij = 0 if and only if (i, j) / ∈ E. Each node of the graph corresponds to an agent of the multi-agent system, and the edges represent the available communication channels (i.e. (i, j) ∈ E if and only if agents i and j can share information). We assume that there is no edges connecting a node with itself, i.e. a ii = 0, for all i ∈ V; and that the communication channels are bidirectional, i.e. a ij = a ji . The last assumption implies that G is undirected. Additionally, we denote by N i = { j ∈ V : (i, j) ∈ E}, the set of neighbours of node i, i.e. the set of nodes that are able to receive/send information from/to node i.
Let us define the n × n matrix L(G) = [l i j ], known as the graph Laplacian of G, as follows: Properties of L(G) are related to connectivity characteristics of G as shown in the following theorem. We remark that a graph G is said to be connected if there exists a path connecting any pair of nodes. From Equation (1), it can be verified that L(G)1 = 0, . . . , 0] . A consequence of this fact is that L(G) is a singular matrix. However, we can modify L(G) to obtain a nonsingular matrix as shown in the following lemma. Lemma 2.1: Let L k r (G) ∈ R (n−1)×n be the submatrix obtained by removing the kth row of the graph Laplacian L(G), and let L k (G) ∈ R (n−1)×(n−1) be the submatrix obtained by removing the kth column of L k r (G). If G is connected, then L k (G) is positive definite. Furthermore, the inverse matrix of L k (G) satisfies (L k (G)) −1 l k r k = −1, where l k r k is the kth column of the matrix L k r (G). Proof: First, notice that L(G) is a symmetric matrix because G is an undirected graph. Moreover, notice that according to Equation (1), L(G) is diagonally dominant with non-negative diagonal entries. The same holds for L k (G) since this is a sub-matrix obtained by removing the kth row and column of L(G). Thus, to show that L k (G) is positive definite, it is sufficient to prove that L k (G) is nonsingular.
According to Theorem 2.1, since G is connected, L(G) has exactly n − 1 linearly independent columns (resp. rows). Let us show that the kth column (resp. row) of L(G) can be obtained by a linear combination of the other columns (resp. rows), i.e. the kth column (resp. row) is not linearly independent of the rest of the columns (resp. rows).
Since L(G)1 = 0, notice that l ik = − j∈V, j =k l i j , for all i ∈ V, i.e. the kth column can be obtained by a linear combination of the rest of the columns. Furthermore, since L(G) is a symmetric matrix, the same occurs with the kth row. Therefore, the submatrix L k (G) is nonsingular since its n − 1 columns (resp. rows) are linearly independent. Now, let us prove that (L k (G)) −1 l k r k = −1. To do so, we use the fact that (L k (G)) −1 L k (G) = I, where I is the identity matrix. Hence, by the definition of matrix multiplication, we have that where l k i j andl k i j are the elements located in the ith row and jth column of the matrices L(G) and L k (G) −1 , respectively. Thus, Let l k r k m be the mth entry of the vector l k r k . Notice that, according to the definition of L k (G) and since L(G)1 = 0, l k mi = − n−1 j=1, j =i l k m j − l k r k m . Replacing this value in Equation (3), we obtain According to Equation (2), k im l k r k m = −1, for all i = 1, … , n − 1. Therefore, L k (G) −1 l k r k = −1. Theorem 2.1 and Lemma 2.1 will be used in the analysis of the method proposed in this paper.

Problem statement
In general terms, a resource allocation problem can be formulated as follows (Patriksson, 2008;Patriksson & Strömberg, 2015): where x i ∈ R is the resource allocated to the ith zone; x = [x 1 , … , x n ] ; φ i : R → R is a strictly convex and differentiable cost function; X is the available resource; and x i , is the lower bound of x i , i.e. the minimum amount of resource that has to be allocated in the ith zone. Given the fact that we are interested in distributed algorithms to solve the problem stated in Equation (4), we consider a multi-agent network, where the ith agent is responsible for managing the resource allocated to the ith zone. Moreover, we assume that the agents have limited communication capabilities, so they can only share information with their neighbours. This constraint can be represented by a graph G = {V, E} as it was explained in Section 2.
Avoiding the individual inequality constraints (4c), KKT conditions establish that at the optimal solution Hence, a valid alternative to solve (4a-4b) is the use of consensus methods. For instance, we can adapt the algorithm presented in Xiao and Boyd (2006), which is described as follows:ẋ This algorithm has two main properties: if the nodes i and j are connected by a path; is the initial condition of x i . Therefore, if the graph G is connected and the initial condition is feasible (i.e. n i=1 x i (0) = X), x asymptotically reaches the optimal solution of (4a-4b) under (5). However, the same method cannot be applied to solve (4) (the problem that considers lower bounds in the resource allocated to each zone) since some feasibility issues related with the constraints (4c) arise.
In the following section, we propose a novel method that extends the algorithm in Equation (5) to deal with the individual inequality constraints given in Equation (4c).

Resource allocation among a subset of nodes in a graph
First, we consider the following subproblem: let G = {V, E} be a graph comprised by a subset of active nodes V a and a subset of passive nodes V p , such that V a V p = V. A certain amount of resource X has to be split among those nodes to minimise the cost function φ(x) subject to each passive node is allocated with its corresponding lower bound x i . Mathematically, we formulate this subproblem as: Feasibility of (6) is guaranteed by making the following assumption.
Assumption 4.1: At least one node is active, i.e. V a = ∅.
According to KKT conditions, the active nodes have to equalise their marginal costs at the optimal solution. Therefore, a consensus among the active nodes is required to solve (6). Nonetheless, classic consensus algorithms, as the one given in Equation (5), cannot be used directly. For instance, if all the nodes of G apply (5) and G is connected, the marginal costs of both passive and active nodes are driven to be equal in steady state. This implies that the resource allocated to passive nodes can violate the constraint (6c). Besides, if the resource allocated to passive nodes is forced to satisfy (6c) by setting x * i = x i , for all i ∈ V p , there is no guarantee that the new solution satisfies (6b). Another alternative, is to apply (5) to only active nodes (in this case, the neighbourhood of node i ∈ V a in Equation (5) has to be taken as { j ∈ V a : (i, j) ∈ E}, and the initial condition must satisfy i∈V a x i (0) = X − i∈V p x i ). However, the sub-graph formed by the active nodes is not necessarily connected although G is connected. Hence, marginal cost of active nodes are not necessarily equalised at equilibrium, which implies that the obtained solution is sub-optimal. In conclusion, modification of (5) to address (6) is not trivial. In order to deal with this problem, we propose the following algorithm: In the same way as in (5), the variables {x i , i ∈ V} in Equation (7) correspond to the resource allocated to both active and passive nodes. Notice that we have added auxiliary variables {x i , i ∈ V p } that allow the passive nodes to interact with their neighbours taking into account the constraint (6c). On the other hand, the term j∈N i (y j − y i ), in Equations (7a)-(7b), leads to a consensus among the elements of the vector y = [y 1 , … , y n ] , which are given in Equation (7c). For active nodes, y i only depends on the marginal cost φ i (x i ), while for passive nodes, y i depends on both the marginal cost and the state of the auxiliary variablex i . Therefore, if the ith node is passive, it has to compute both variables x i andx i . Furthermore, it can be seen that, if all the nodes are active, i.e. (V a = V), then the proposed algorithm becomes the one stated in Equation (5).
Notice that the ith node only needs to know y i and the values is a distributed map over the graph G (Cortés, 2008). This implies that the dynamics given in Equation (7) can be computed by each node using only local information. In fact, the message that the ith node must send to its neighbours is solely composed by the variable y i .

.. Feasibility
Let us prove that, under the multi-agent system proposed in Equation (7), x(t) satisfies the first constraint of the problem given by Equation (6), for all t ࣙ 0, provided that Moreover, according to Equation (7), The above lemma does not guarantee that x(t) is always feasible because of the second constraint in Equation (6), i.e. x i = x i , for all i ∈ V p . However, it is possible to prove that, at equilibrium, this constraint is properly satisfied.

.. Equilibrium point
The next proposition characterises the equilibrium point of the multi-agent system given in Equation (7).

Proposition 4.1: If G is connected, the system in Equation
4.1: Proposition 4.1 states that, at the equilibrium point of (7), the active nodes equalise their marginal costs, while each passive node is allocated with an amount of resource equal to its corresponding lower bound. In conclusion, if n i=1 x * i = X, then it follows from Proposition 4.1, that x * minimises the optimisation problem given in Equation (6). Additionally, notice that the values {x * i , i ∈ V p } are equal to the KKT multipliers associated with the constraint (6c).

.. Convergence
Let us prove that the dynamics in Equation (7) converge and Assumption 4.1 holds,then x(t) converges to x * under Equation (7), where x * is the solution of the optimisation problem stated in Equation (6) Proof: According to Lemma 4.1, since n i=1 x i (0) = X, then x(t) satisfies the first constraint of the problem stated in Equation (6), for all t ࣙ 0. Therefore, it is sufficient to prove that the equilibrium point x * , {x * i , i ∈ V p } (which is given in Proposition 4.1) of the system proposed in Equation (7) is asymptotically stable (AS). In order to do that, let us express our multi-agent system in error coordinates, as follows: where ; e y = [e y 1 , . . . , e y n ] ; and L(G)e y i represents the ith element of the vector L(G)e y .
Since Assumption 4.1 holds, V a = ∅. Let k be an active node, i.e. k ∈ V a , and let e k , e k y be the vectors obtained by removing the kth element from vectors e and e y , respectively. We notice that, according to Lemma , for all t ࣙ 0. Therefore, Equation (8) can be expressed aṡ where L k (G) and l k r k are defined in Lemma 2.1. In order to prove that the origin of the above system is AS, let us define the following Lyapunov function (adapted from Obando, Quijano, & Rakoto-Ravalontsalama, 2014): The function V is positive definite since G is connected (the reason of this fact is that, according to Lemma 2.1, L k (G) and its inverse are positive definite matrices if G is connected). The derivative of V along the trajectories of the system stated in Equation (9) is given by, where φ i is strictly increasing given the fact that Given the fact that G is connected and V = V p (by Assumption 4.1), thenė = 0 iff e y = 0 (see Equation (8)). Therefore, the only solution that stays identically in S is the trivial solution, i.e. e i (t) = 0, for all i ∈ V, e i (t ) = 0, for all i ∈ V p . Hence, we can conclude that the origin is AS by applying the Lasalle's invariance principle.
In summary, we have shown that the algorithm described in Equation (7) asymptotically solves the subproblem in Equation (6), i.e. (7) guarantees that the resource allocated to each passive node is equal to its corresponding lower bound, while the remaining resource X − i∈V p x i is optimally allocated to active nodes.

Optimal resource allocation with lower bounds
Now, let us consider our original problem stated in Equation (4), i.e. the resource allocation problem that includes lower bound constraints. Let x * = [x * 1 , . . . , x * n ] be the optimal solution of this problem. Notice that, if we know in advance which nodes will satisfy the constraint (4c) with strict equality after making the optimal resource allocation process, i.e. I := {i ∈ V : x * i = x i }, we can mark these nodes as passive and reformulate (4) as a subproblem of the form (6). Based on this idea, we propose a solution method for (4), which is divided in two stages: in the first one, the nodes that belong to I are identified and marked as passive; in the second one, the resulting subproblem of the form (6) is solved by using (7).
Protocol (7) can be also used in the first stage of the method as follows: in order to identify the nodes that will satisfy (4c) with strict equality at the optimal allocation, we start marking all nodes as active and apply the resource allocation process given by (7). The nodes that are allocated with an amount of resource below their lower bounds at equilibrium are marked as passive, and then (7) is newly applied (in this way, passive nodes are forced to meet (4c)). This iterative process is performed until all nodes satisfy their lower bound constraints. Notice that the last iteration of this procedure corresponds to solve a subproblem of the form (6) where the set of passive nodes is equal to the set I. Therefore, this last iteration is equivalent to the second stage of the proposed method.
Summarising, our method relies on an iterative process that uses the continuous-time protocol (7) as a subroutine. The main idea of this methodology is to identify in each step the nodes that have an allocated resource out of their lower bounds. These nodes are marked as passive, so they are forced to satisfy their constraints in subsequent iterations, while active nodes seek to equalise their marginal costs using the remaining resource. In the worst case scenario, the classification between active and passive nodes requires |V| iterations, where |V| is the number of nodes in the network. This fact arises when only one active node becomes passive at each iteration.
The proposed method is formally described in Algorithm 1. Notice that this algorithm is fully decentralised since Steps 4-6 can be computed by each agent using only local information.
Step 4 corresponds to solve Equation (7), while Steps 5 and 6 describe the conditions for converting an active node into passive. Let us note that Steps 4-6 have to be performed |V| times since we are considering the worst case scenario. Therefore, each agent needs to know the total number of nodes in the network. This requirement can be computed in a distributed way by using the method proposed in Garin and Schenato (2010, p. 90). We also notice the fact that the agents have to be synchronised (as usual in several distributed algorithms; Cortés, 2008;Garin & Schenato, 2010;Xiao & Boyd, 2006) in order to apply the Step 4 of Algorithm 1, i.e. all agents must start solving Equation (7) at the same time.

Input:
-Parameters of the problem in Equation (4).
According to the reasoning described at the beginning of this subsection, we ideally require to know the steady-state solution of Equation (7) at each iteration of Algorithm 1 (since we need to identify which nodes are allocated with an amount of resource below their lower bounds in steady state). This implies that the time t l in Step 4 of Algorithm 1 goes to infinity. Under this requirement, each iteration would demand infinite time and the algorithm would not be implementable. Hence, to relax the infinite time condition, we state the following assumption on the time t l .

Assumption 4.2: Let x *
i,l be the steady state of x i (t) under Equation (7), with initial conditions x(0) =x i,l−1 , V a = V a,l−1 , V p =Ṽ p,l−1 , and {x i (0) = 0, ∀i ∈ V p } 1 . For each l = 1, . . . , |V| − 1, the time t l satisfies the following condition: According to assumption 4.2, for the first |V| − 1 iterations, we only need a solution of (7) that is close enough to the steady-state solution. We point out the fact that, if the conditions of Proposition 4.2 are met in the lth iteration of Algorithm 1, then x i (t) asymptotically converges to x * i,l , for all i ∈ V, under Equation (7). Therefore, Assumption 4.2 is satisfied for large values of t 1 , . . . , t |V|−1 .
Taking into account all the previous considerations, the next theorem states our main result regarding the optimality of the output of Algorithm 1.

Theorem 4.1:
Assume that G is a connected graph. Moreover, assume that φ i is a strictly convex function for all i = 1, … , n. If t 1 , . . . , t |V|−1 satisfy Assumption 4.2, and the problem stated in Equation (4) is feasible, then the output of Algorithm 4 tends to the optimal solution of the problem given in Equation (4) as t |V| → ∞.
P8: i∈Ṽ a,l x * i,l ≥ i∈Ṽ a,l x * i,l+1 (we prove P8 as follows: using P1 and the result in Lemma 4.1, we know that i∈V x * i,l = i∈V x * i,l+1 = X. Moreover, according to P5, V can be expressed as V =Ṽ a,l Ṽ p,l , whereṼ p,l−1 ⊂Ṽ p,l (see P3).
, for some i ∈Ṽ a,l . According to Proposition 4.2, and since P1 and P9 hold, x * i,l has the characteristics given in Proposition 4.1, for all i ∈ V, and for all l = 1, . . . , |V|. Hence, φ i (x * i,l ) has the same value for all i ∈Ṽ a,l−1 , and φ i (x * i,l+1 ) has the same value for all i ∈Ṽ a,l . Moreover, sinceṼ a,l ⊂Ṽ a,l−1 (according to P3), we have that , for all i ∈ V a,l , because φ i is strictly increasing (this follows from the fact that φ i is strictly convex by assumption). Therefore, i∈Ṽ a,l x * i,l < i∈Ṽ a,l x * i,l+1 , which contradicts P8). Now, let us prove that {x * 1,|V| , . . . , x * n,|V| } solves the Problem in Equation (4). First, the solution {x * 1,|V| , . . . , x * n,|V| } is feasible according to P1 and P7. On the other hand, from P9, it is known that ∃k : k ∈Ṽ a,l , for all l = 1, . . . , = λ (this follows from the fact that φ j (x * j,l ) has the same value for all j ∈Ṽ a,l−1 , which in turn follows directly from step 4 of Algorithm 1,and Proposition 4.2).

Early stopping criterion
Notice that, if the set of passive nodes does not change in the kth iteration of Algorithm 1 because all active nodes satisfy the lower bound constraints (see step 5), then the steady state solutions x * i,k and x * i,k+1 are the same, for all i ∈ V, which implies that the set of passive nodes also does not change in the (k + 1)th iteration. Following the same reasoning, we can conclude that Therefore, in this case, {x * 1,k , . . . , x * n,k } is the solution of our resource allocation problem. Practically speaking, this implies that Algorithm 1 does not need to perform more iterations after the kth one. Thus, it is possible to implement a flag z * i (in a distributed way) that alerts the agents if all active nodes satisfy the lower bound constraints after step 4 of Algorithm 1. A way to do that is by applying a min-consensus protocol (Cortés, 2008) with initial conditions z i (0) = 0 if the node i is active and does not satisfy its lower bound constraint, and z i (0) = 1 otherwise. Hence, notice that our flag z * i (i.e. the result of the min-consensus protocol) is equal to one, for all i ∈ V, only if all the active nodes satisfy the lower bound constraints, which corresponds to the early stopping criterion described above.

Simulation results and comparison
In this section, we compare the performance of our algorithm with other continuous-time distributed techniques found in the literature. We have selected three techniques that are capable to address nonlinear problems and can handle lower bound constraints: (i) a distributed interior point method (Xiao & Boyd, 2006), (ii) the local replicator equation (Pantoja & Quijano, 2012), and (iii) a distributed interior point method with exact barrier functions (Cherukuri & Cortés, 2015). The first one is a traditional methodology that uses barrier functions; the second one is a novel technique based on population dynamics; and the third one is a recently proposed method that follows the same ideas as the first one, but replaces classic logarithmic barrier functions by exact penalty functions. Below, we briefly describe the aforementioned algorithms.

Distributed interior point (DIP) method
This algorithm is a variation of the one presented in Equation (5) that includes strictly convex barrier functions to prevent the solution to flow outside the feasible region. The barrier functions b i (x i ) are added to the original cost function as follows: where φ b (x) is the new cost function, and ϵ > 0 is a constant that minimises the effect of the barrier function when the solution is far from the boundary of the feasible set. With this modification, the distributed algorithm is described by the following equation: is equal to the marginal cost plus a penalty term induced by the derivative of the corresponding barrier function.

Local replicator equation (LRE)
This methodology is based on the classical replicator dynamics from evolutionary game theory. In the LRE, the growth rate of a population that plays a certain strategy only depends on its own fitness function and on the fitness of its neighbours. Mathematically, the LRE is given byẋ where v i is the fitness perceived by the individuals that play the ith strategy. In this case, the strategies correspond to the nodes of the network, and the fitness functions to the negative marginal costs (the minus appears because replicator dynamics are used to maximise utilities instead of minimise costs). On the other hand, it can be shown that, if the initial condition x(0) is feasible for the problem given in Equation (4), then x(t) remains feasible for all t ࣙ 0, under the LRE.

Distributed interior point method with exact barrier functions (DIPe)
This technique follows the same reasoning of the DIP algorithm. The difference is that DIPe uses exact barrier functions (Bertsekas, 1975) to guarantee satisfaction of the lower bound constraints. The exact barrier function for the ith node is given by: is the feasible region of x for the problem (4). Using these exact barrier functions, the augmented cost function can be expressed as: The DIPe algorithm is given in terms of the augmented cost function and its generalised gradient ∂φ e b (x) = [∂ 1 φ e b (x), . . . , ∂ n φ e b (x)] as follows: where In Cherukuri and Cortés (2015), the authors show that the differential inclusion (13) converges to the optimal solution of the problem (4), provided that x(0) is feasible.

Comparison
In order to compare the performance of our algorithm with the three methods described above, we use the following simulation scenario: a set of n nodes connected as in Figure 1 (we use this topology to verify the behaviour of the different algorithms in the face of few communication channels since previous studies have shown that algorithms' performance decreases with the number of where a i and b i are random numbers that belong to the intervals [1,2] and [− 1 2 , 1 2 ], respectively; a resource constraint X = 1; and a set of lower bounds {x i = 0 : i ∈ V}. For each n, we generate 50 problems with the characteristics described above. The four distributed methods are implemented in Matlab employing the solver function ode23s. Moreover, we use the solution provided by a centralised technique as reference. The results on the average percentage decrease in the cost function reached with each algorithm and the average computation time (time taken by each algorithm for solving a problem 2 ) are summarised in Table 1. Results of DIPe for 100 and 200 nodes were not computed for practicality since the time required by this algorithm to solve a 100/200-nodes problem is very high.
We notice that the algorithm proposed in this paper always reaches the maximum reduction, regardless of the number of nodes that comprise the network. The same happens with the DIPe algorithm. This is an important advantage of our method compared to other techniques. In contrast, the algorithm based on the LRE performs far from the optimal solution. This unsatisfactory behaviour is due to the small number of links of the considered communication network. In Pantoja and Quijano (2012), the authors prove the optimality of the LRE in problems involving well connected networks; however, they also argue that this technique can converge to suboptimal solutions in other cases. On the other hand, the DIP method provides solutions close to the optimum. Nonetheless, its performance decreases when the number of nodes increases. This tendency is due to the influence of barrier functions on the original problem. Notice that, the larger the number of nodes, the bigger the effect of the barrier functions in Equation (11).
Regarding the computation time, although convergence of the proposed method is slower than the one shown by LRE and DIP, it is faster than the convergence of the method based on exact barrier functions, i.e. DIPe. Therefore, among the methods that guarantee optimality of the solution, our technique shows the best convergence speed. Computation time taken by DIPe is affected by the use of penalty terms that generate strong changes in the value of the cost function near to the boundaries of the feasible set. The drastic variations of the generalised gradient of exact barrier functions produces oscillations of numerical solvers around the lower bounds (a visual inspection of the results given in Figure 3 of Cherukuri and Cortés (2015) confirms this claim). These oscillations are the main responsible for the low convergence speed shown by DIPe. On the other hand, LRE and DIP exhibit the fastest convergence. Hence, LRE and DIP are appealing to be implemented in applications that require fast computation and tolerate suboptimal solutions.

Applications
This section describes the use of the approach developed in this paper to solve two engineering problems. First, we present an application for sharing load in multiple chillers plants. Although this is not a large-scale application (multi-chiller plants are typically comprised of less than ten chillers; Yu & Chan, 2007), it aims to illustrate the essence of the proposed method and shows algorithm's performance in small-size problems. One of the reasons to use a distributed approach in small-/medium-size systems is due to the need of enhancing systems resilience in the face of central failures (e.g. in multiple chiller plants, central failures can occur due to cyber-attacks (Manic, Wijayasekara, Amarasinghe, & Rodriguez-Andina, 2016) against building management systems (Yu & Chan, 2007)). The second application deals with the distributed computation of the Euclidean projection of a vector onto a given set. Particularly, we use the proposed algorithm as part of a distributed technique that computes optimal control inputs for plants composed of a large number of sub-systems. This application aims to illustrate the performance of the method proposed in this paper when coping with large-scale problems.

Optimal chiller loading
The optimal chiller loading problem in multiple chiller systems arises in decoupled chilled-water plants, which are widely used in large air-conditioning systems (Chang & Chen, 2009). The goal is to distribute the cooling load among the chillers that comprise the plant for minimising the total amount of power used by them. For a better understanding of the problem, below we present a brief description of the system. A decoupled chilled-water plant comprised by n chillers is depicted in Figure 2. The purpose of this plant is to provide a water flow f T at a certain temperature T s to the rest of the air-conditioning system. In order to do this task the plant needs to meet a cooling load C L that is given by the following expression: where m > 0 is the specific heat of the water, and T r is the temperature of the water returning to the chillers. Since there are multiple chillers, the total cooling load C L is split among them, i.e. C L = n i=1 Q i , where Q i is the cooling power provided by the ith chiller, which, in turn, is given by where f i > 0 and T i are, respectively, the flow rate of chilled water and the water supply temperature of the ith chiller. As it is shown in Figure 2, we have that f T = n i=1 f i . In order to meet the corresponding cooling load, the ith chiller consumes a power P i that can be calculated using the following expression (Chang & Chen, 2009): where k j, i , for j = 0, … , 6, are constants related to the ith chiller. If we assume that the flow rate f i of each chiller is constant, then P i is a quadratic function of the temperature T i . The optimal chiller loading problem involves the calculation of the chillers' water supply temperatures that meet the total cooling load given in Equation (14), and minimise the total amount of power consumed by the chillers, i.e. n i=1 P i . Moreover, given the fact that each chiller has a maximum cooling capacity, we have to consider the following additional constraints: where Q i is the maximum capacity (rated value) of the ith chiller. Summarising, the optimal chiller loading problem can be expressed as follows: Now, let us consider that we want to solve the aforementioned problem in a distributed way by using a multiagent system, in which each chiller is managed by an agent that decides the value of the water supply temperature. We assume that the ith agent knows (e.g. by measurements) the temperature of the water returning to the chillers, i.e. T r , and the flow rate of chilled water, i.e. f i . Moreover, agents can share their own information with their neighbours through a communication network with a topology given by the graph G. If each P i (T i ) is a convex function, then the problem can be solved by using the method proposed in Algorithm 1 (we take, in this case, x i = f i T i ). The main advantage of this approach is to increase the resilience of the whole system in the face of possible failures, due to the fact that the plant operation does not rely on a single control centre but on multiple individual controllers without the need for a centralised coordinator.

.. Illustrative example
We simulate a chilled-water plant comprised by 7 chillers. 3 The cooling capacity and the water flow rate of each chiller are, respectively, Q i = 1406.8 kW, and f i = 65 kg.s −1 , for i = 1, … , 7; the specific heat of the water is m = 4.19 kW.s.kg −1 . degC −1 ; the supply temperature of the system is T s = 11 degC; and the coefficients k j, i of Equation (16) are given in Table 2. We operate the system at two different cooling loads, the first one is 90% of the total capacity, i.e. C L = 0.9 n i=1 Q i , and the second one is 60% of the total capacity, i.e. C L = 0.6 n i=1 Q i . The P i -T i curves are shown in Figure 3(a) for both cases, it can be noticed that all functions are convex. In order to apply the distributed solution presented in Algorithm 1, we use an agent per chiller (i.e. the ith agent controls the supply temperature T i of the ith chiller) and the communication network shown in Figure 1. In all cases the initial conditions of the chillers' supply temperatures are T i (0) = T s , for i = 1, … , 7. The results for the first cooling load, i.e. C L = 8862.8 kW, are depicted in Figure 3(b), where it is shown that the cooling load is properly allocated among the chillers by adjusting the supply temperatures. More efficient chillers (i.e. chiller 3, chiller 6, and chiller 7 in Figure 3(a)) are more loaded than the less efficient ones (i.e. chiller 2 and chiller 5). This can be noticed from the fact that their supply temperatures, in steady state, reach the minimum value. Furthermore, the energy consumption is minimised and power saving reaches to 2.6%. The results for the second cooling load, i.e. C L = 5908.6 kW, are shown in Figure 3(c), where it can be noticed a similar performance to that obtained with the first cooling load. However, in this case, it is not necessary that the supply temperatures reach the minimum value to meet the required load. Newly, energy consumption is minimised and power saving reaches to 2.8%. As it is stated in Section 4, convergence and optimality of the method is guaranteed under the conditions given in Theorem 4.1.
In both cases we use the early stopping criterion given in Section 4. performance even considering that the communication graph is sparse and the optimal solution is not in the interior of the feasible domain. As shown in Section 5, this characteristic is an advantage of Algorithm 1 over population dynamics techniques as the one proposed in Barreiro-Gomez, Obando, Ocampo-Martinez, and Quijano (2015) to compute the Euclidean projection in a distributed way.

Discussion
The method developed in this paper solves the problem of resource allocation with lower bounds given in Equation (4). The main advantage of the proposed technique is its distributed nature; indeed, our approach does not need the implementation of a centralised coordinator. This characteristic is appealing, especially in applications where communications are strongly limited. Moreover, fully distributed methodologies increase the autonomy and resilience of the system in the face of possible failures. In Section 5, we show by means of simulations that the performance of the method presented in this paper does not decrease when the number of nodes (which are related to the decision variables of the optimisation problem) is large, or the communication network that allows the nodes to share information has few channels.
In these cases, the behaviour of our approach is better than the behaviour of other techniques found in the literature, such as the DIP method, or the LRE. Moreover, it is worth noting that our technique addresses the constraints as hard. This fact has two important consequences: (i) in all cases, the solution satisfies the imposed constraints, and (ii) the objective function (and therefore the optimum) is not modified (contrary to the DIP method that includes the constraints in the objective function decreasing the quality of the solution as shown in Section 5.4).
Other advantage of the method proposed in this paper is that it does not require an initial feasible solution of the resource allocation problem (4). Similarly to the DIPe technique, our method only requires that the starting point satisfies the resource constraint (4b), i.e. we need that n i=1 x i (0) = X. Notice that an initial solution x(0) that satisfies (4b) is not hard to obtain in a distributed manner. For instance, if we assume that only the kth node has the information of the available resource X, we can use (x k (0) = X, {x i (0) = 0 : i ∈ V, i = k}) as our starting point. Thus, an initialisation phase is not required. In contrast, other distributed methods, such as DIP and LRE needs an initial feasible solution of the problem (4), i.e. a solution that satisfies (4b) and (4c). Finding this starting point is not a trivial problem for systems involving a large number of variables. Therefore, for these methods, it is necessary to employ distributed constraint satisfaction algorithm (as the one described in Domınguez-Garcıa & Hadjicostis, 2011) as a first step.
On the other hand, we notice that to implement the early stopping criterion presented at the end of Section 4, it is required to perform an additional min-consensus step in each iteration. Despite this fact, if the number of nodes is large, this criterion saves computational time, because in most of the cases, all passive nodes are identified during the first iterations of Algorithm 1.

Conclusions
In this paper, we have developed a distributed method that solves a class of resource allocation problems with lower bound constraints. The proposed approach is based on a multi-agent system, where coordination among agents is done by using a consensus protocol. We have proved that convergence and optimality of the method is guaranteed under some mild assumptions, specifically, we require that the cost function is strictly convex and the graph related to the communication network that enables the agents to share information is connected. The main advantage of our technique is that it does not need a centralised coordinator, which makes the method appropriate to be applied in large-scale distributed systems, where the inclusion of centralised agents is undesirable or infeasible. As future work, we propose to use a switched approach in order to eliminate the iterations in Algorithm 1. Moreover, we plan to include upper bound constraints in our original formulation.