Load-balancing for multi-skilled servers with Bernoulli routing

Annals of Operations Research (2022) 312:949–971 https://doi.org/10.1007/s10479-022-04532-7 ORIGINAL RESEARCH Load-balancing for multi-skilled servers with Bernoulli routing Fernando Miguelez1·Josu Doncel2·Balakrishna J. Prabhu3 Accepted: 3 January 2022 / Published online: 15 January 2022 © The Author(s) 2022 Abstract We study the optimal Bernoulli routing in a multiclass queueing system with a dedicated server for each class as well as a common (or multi-skilled) server that can serve jobs of all classes. Jobs of each class arrive according to a Poisson process. Each server has a holding cost per customer and use the processor sharing discipline for service. The objective is to minimize the weighted mean holding cost. First, we provide conditions under which classes send their traffic only to their dedicated server, only to the common server, or to both. A fixed point algorithm is given for the computation of the optimal solution. We then specialize to two classes and give explicit expressions for the optimal loads. Finally, we compare the cost of multi-skilled server with that of only dedicated or all common servers. The theoretical results are complemented by numerical examples that illustrate the various structural results as well as the convergence of the fixed point algorithm. Keywords Bernoulli routing ·Parallel-servers ·Multi-skilled servers 1 Introduction 1.1 Motivation We investigate the performance of a multi-skilled queueing system formed by parallel servers with Processor Sharing (PS) queues. Jobs of different classes of customers arrive to the system following a Poisson process. There is one dedicated server for each class of customer and one multi-skilled server that can execute jobs of all classes. Furthermore, we assume that jobs are assigned to the servers according to the Bernoulli policy. Our goal is to find the optimal load balancing so as to minimize the weighted mean number of jobs in the system. The main application of our model comes from wireless networks. Consider a region divided in different subregions. Each dedicated server models an antenna that provides serBJosu Doncel [email protected] 1Universidad Publica de Navarra, Navarra, Spain 2University of the Basque Country, UPV/EHU, Leioa, Spain 3LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France 123 950 Annals of Operations Research (2022) 312:949–971 vice to a unique subregion and the multi-skilled server models a central antenna that provides service to all the subregions. Using the results of this article, one can determine how the traffic of each subregion must be shared between the antenna of that region and the central one in order to minimize the performance of the system. This architecture has been previously considered by Taboada et al. (2017) in a different context where dedicated servers (or microcells in their model) can be switched on and off so as to minimize the weighted sum of the mean delay and the mean power consumption in the system. 1.2 Related work Load balancing has been widely investigated in different contexts. In data centers, for example, various policies depending upon the information available to the dispatcher have been proposed. In general, optimal policies for the typical performance measures such as mean processing times are not easy to determine albeit in some specific cases. For example, when no information on the state of the servers is available, the optimal Bernoulli routing policy was determined in Altman et al. (2011) for mono-skilled servers only. For FCFS servers, a policy based on Sturm sequences (Gaujal et al. 2006) are known to be optimal. With more information on the server state, a number of heuristics such as Join the Shorter of dqueues (Mitzenmacher 2001; Vvedenskaya et al. 1996) and Join the Shortest Queue (Graham 2000) have been analyzed in the large server asymptotic case. In addition, there are various pullbased policies such as Join the Idle Queue (Lu et al. 2011) that are known to work well in practice. Another important routing policy is the Size Interval Task Assignment (HarcholBalter et al. 1999) where jobs of different sizes are executed in different servers and, therefore, the service requirement of incoming tasks need to be known. This policy has been further studied in Feng et al. (2005) and the author in Harchol-Balter (2000) presented a variation of this policy in which the size of jobs does not need to be known. Load balancing has also been investigated for balancing energy costs in data centers using Energy Packet Networks model (Fourneau 2020), whereas in Liu et al. (2015) it is considered that data centers are located in different geographical zones. The above works are mostly concerned with mono-skilled or homogeneous servers. In networks with multi-skilled agents or servers, skill-based routing policies have been proposed and investigated (Koole et al. 2003; Wallace and Whitt 2005). These works are mainly oriented towards call-center architectures with Erlang-B or Erlang-C type of queues. An illustrative example is an overflow-type policy, where each incoming call has a list of agents ordered by priority, with highest priority given to mono-skilled ones, and is routed to the first available agent of this list. If no agent is available, the call can be queued or blocked depending on the architecture. These routing policies are usually difficult to analyze and the cited works are interested in approximations for the various performance measures for a given policy. In these models, obtaining the optimal policy analytically is not easy. We refer to Chen et al. (2020) for a recent survey on multi-skilled systems. Multi-skilled queues appear also in the analysis of redundancy systems (Gardner et al. 2015; Bonald et al. 2017) in which incoming requests can be sent simultaneously to a subset of queues. We do not investigate the redundancy aspect. The network topology we consider makes our model different from Altman et al. (2011) in which all servers can execute all type of tasks. 123 Annals of Operations Research (2022) 312:949–971 951 1.3 Contributions The main contributions of the article are summarized as follows: •We provide a necessary and sufficient condition for the stability of the system. •We show that the optimization problem in terms of probabilities can be reformulated in terms of the loads of the servers, which is a convex problem with linear constraints. We also show that if there are routing probabilities that satisfy the original problem (that we will call (PROB-OPT)), then it is possible to find loads that satisfy the reformulated problem (that we will call (LOAD-OPT)). •We fully characterize the optimal loads on the servers for two classes of customers. For more than two customers, we provide in Proposition 2conditions under which each class of traffic satisfies one of the following: (i) it sends all its traffic to its dedicated server, (ii) it sends all its traffic to the multi-skilled server and (iii) it shares its traffic among the multi-skilled and its dedicated server. •Using the result of Proposition 2, we present a fixed-point algorithm whose convergence ensures that the optimal loads on the servers are achieved. This algorithm starts with an initial condition of the set of servers (according to one of the three possible traffic sharing policies of Proposition 2) and its fixed point is given by the partition of the set of servers. Providing an analytical proof of this convergence on the partition of the set of servers seems to be an extremely difficult task. However, we illustrate the convergence of this algorithm using numerical experiments. •We compare the performance of our model with the performance of two models. The first model consists of a system where all the servers are multi-skilled and we show the existence of a switching curve, i.e., when the arrival rate of one of the traffic increases, the model whose performance is better changes. The second model consists of a system with no sharing, that is, all the servers are dedicated or mono-skilled, and we provide conditions on the arrival rates such that the performance of the no-sharing model is larger than the performance of our model. •We delve into the comparison of the aforementioned models using numerical experiments. First, we show the uniqueness of the switching curve when we compare our model with a system where all the servers are multi-skilled. We also observe that, in a system formed by servers with equal capacity and different (but not extremely large) holding costs, the region where our model outperforms the all-sharing system is very large. 1.4 Organization In the next section, we describe the network model and define the optimization problem. Section 3gives the stability condition and presents an equivalent problem in terms of loads on the servers. In Sect. 4, the main results on the structure of the optimal policy are provided for the model under consideration. We compare in Sect. 5the performance of our model with the performance of models with other network topologies. We present our numerical experiments in Sect. 6. Finally, we discuss our main conclusions in Sect. 7. 123 952 Annals of Operations Research (2022) 312:949–971 Fig. 1 The model under study in this article 2 Model description 2.1 Notation We consider a server farm with Processor Sharing (PS) queues and an input traffic of different classes. Let K={1,2,...,C}be the set of classes. We assume that jobs of class i∈K arrive to the system according to a Poisson process and have generally distributed service times1.Letηibe the traffic intensity of jobs of class i. The class of a job defines the set of servers that can be assigned to this job (Fig. 1). We consider a system with C+1 servers. Let S={0,1,...,C}be the set of servers. For a server j∈S,wedenotebyrjthe capacity or speed of Server j(i.e., the amount of traffic that can be served per unit of time) and by cjits holding cost. We denote by Sithe set of servers that can execute jobs of class iand, for A⊂K,SA=∪ i∈ASi.For j=1,...,C, Server jexecutes jobs of class j, i.e., they are dedicated servers. On the other hand, Server 0 executes jobs of all the classes, i.e., it is a multi-skilled server. For i=1,...,C,wedenotebypithe probability that a class ijob is executed in its dedicated server, i.e., Server i.For j=1,...,C, the load of Server jis defined as follows ρj(p)=ηjpj rj ,(1) whereas for Server 0 as ρ0(p)=j∈Kηj(1−pj) r0 .(2) 1Since our goal is to analyze the mean number of jobs and, in a M/G/1-PS queue, the mean number of jobs depends on the arrival rate and on the service time requirements only through the intensity, we do not specify the arrival rate of each class 123 Annals of Operations Research (2022) 312:949–971 953 2.2 Problem formulation For a given routing strategy p=(pi), the mean number of jobs of server jis denoted by E[Nj(p)]. In this article, we aim to find the routing matrix that minimizes the total cost of the system. More specifically, we analyze the following optimization problem: min p j∈S cjE[Nj(p)](PROB-OPT) 0≤pi≤1,for all i∈K;(3) ηipi<ri,for all i∈K;(4)  i∈K ηi(1−pi)<r0.(5) The first constraint ensures that pi’s are probabilities. The second and third constraints ensure that all the servers are stable, that is, that the total incoming traffic into a server is smaller than its service capacity. 3 Preliminary results We first study the existence of a feasible solution of (PROB-OPT). This is the same as characterizing the conditions under which the system can be stabilized. In the following proposition, we provide this result. Proposition 1 (Stability) The system under consideration can be stabilized if and only if r0+ i∈A ri> i∈A ηi,∀A⊂K.(6) Proof See “Appendix A”.  Since servers are M/G/1-PS queues, we know from Thm 3.8 and Thm 3.9 of Kelly (1979) that the probability of being njobs in Server iis ρn i(1−ρi). Therefore, it follows directly that E[Nj(p)]= ρj 1−ρjand, as a consequence, we can reformulate (PROB-OPT) in terms of the loads on the servers as follows: min æ j∈S cj ρj 1−ρj (LOAD-OPT)  i∈K ηi=r0ρ0+ j∈K rjρj;(7) 0≤ρj<1,for all j∈S;(8)  i∈A ηi≤r0ρ0+ j∈A rjρj,∀A⊂K.(9) We now show that the optimization problems we have considered so far are related. More precisely, we show that, if there are routing probabilities that satisfy (PROB-OPT), then it is possible to find loads that satisfy (LOAD-OPT). 123 954 Annals of Operations Research (2022) 312:949–971 Lemma 1 Let pbe a routing strategy that satisfies (3)–(5). Then, for all j ∈S,ρj(p)also satisfies the constraints of (LOAD-OPT). Proof First, we observe that, if (3)-(5) are satisfied, using (1)and(2), it follows that 0 ≤ ρj<1forall j∈S. We now show that i∈Kηi=r0ρ0+j∈Krjρjin the following way: r0ρ0+ j∈K rjρj= i∈K ηi(1−pi)+ i∈K ηipi= i∈K ηi, where the first equality is given using (1)and(2). Finally, we focus on the constraint i∈Aηi≤r0ρ0+j∈Arjρj,∀A⊂K. Using again (1)and(2), we have for all A⊂Kthat r0ρ0+ j∈A rjρj= i∈K ηi(1−pi)+ i∈A ηipi≥ i∈A ηi(1−pi)+ i∈A ηipi= i∈A ηi. And the desired result follows.  Note that (LOAD-OPT) is a convex problem with linear constraints and has an unique solution as long as the stability condition in Proposition 1is verified. Moreover, from the above lemma, the solution of (PROB-OPT) can be obtained by optimizing directly over the loads. Then, the optimal routing probabilities can be determined later from (1), once the optimal load on each server is determined. Remark 1 Let us remark that we assume that servers are M/G/1-PS queues. However, the results of this article are also valid if we assume that servers are M/M/1 queues with any work-conserving queueing discipline (note that the mean number of customers of a server j in both cases is ρj/(1−ρj)). 4 Analysis of the solution of (LOAD-OPT) Let δj=cj/rj c0/r0for all j∈K.WedenotebyCbthe set of classes that route traffic to two servers, by C0the set of classes that routes all the traffic to Server 0 and by Cdthe set of classes that send all the traffic to its dedicated server. In the following proposition, we present the first result of this section. It gives the conditions under which a class of traffic belongs to Cb,C0or Cd. Proposition 2 Jobs of class i routes traffic to Server 0 if and only if δi>1−ηi ri 1−ρ∗ 0 , and all the traffic of class i is routed to Server 0 if and only if δi≥1 1−ρ∗ 0 , where ρ∗ 0is the optimal load at Server 0and is given by ρ∗ 0=1−r0+j∈Cbrj−j∈Cb∪C0ηj r0+j∈Cbδjrj .(10) 123 Annals of Operations Research (2022) 312:949–971 955 Besides, if j ∈Cdthe optimal load of Server j is ηj rj,if j ∈Cbthe optimal load of Server j is ρ∗ j=1−δj(1−ρ∗ 0)(11) and if j ∈C0the optimal load of Server j is zero. Proof See “Appendix B”.  The above result leads to this corollary which gives a simple sufficient condition to determine when a given class will not send all its traffic to the multi-skilled server. Corollary 1 Let j ∈S.Ifδj<1,then j /∈C0. Proof Since δj<1, we have that the condition δj<1 1−ρ∗ 0 is always satisfied and this implies that j/∈C0according to Proposition 2. From the above corollary, it follows another interesting property that says that, if δj<1 for all j∈K,thenC0=∅. The next result gives an ordering which can help identify classes that use both the dedicated and the multi-skilled server. This can be seen as a way to determine, for a given set of input parameters (arrival rate, server speeds, holding costs, etc.), the skills for which we need to train the multi-skilled servers in order for the system to be optimal. Proposition 3 Let δi≤δj. (a) If i ∈C0,then j ∈C0. (b) If i ∈Cb∪C0and ηi ri≤ηj rj,then j∈Cb∪C0. Proof We first show (a). We consider that i∈C0.Sinceδj≥δi, it follows that δj≥δi≥ 1 1−ρ∗ 0 . Therefore, from Proposition 2,j∈C0. We now show (b). We consider that i∈Cb∪C0.Since ηi ri≤ηj rjand δj≥δi, it follows that δj≥δi≥1−ηi ri 1−ρ∗ 0 ≥ 1−ηj rj 1−ρ∗ 0 . Therefore, from Proposition 2,j∈Cb∪C0. We note that (b)of the above result can be stated as follows: if class jroutes all the traffic to Server j,classiroutes all its traffic to Server iwhen ηi ri≤ηj rjand δi≤δj. In the following result, we show that, under similar conditions, the set of classes that send traffic to two servers can never be {i,j}. Proposition 4 If δi≤δj<1and ηi ri≤ηj r0+rj,thenC bcannot be {i,j}. Proof We assume that Cb={i,j}. For this case, it follows from (10)that 1 1−ρ∗ 0 ≥r0+rjδj+riδi r0+rj+ri−ηj−ηi , where the above inequality is an equality if C0=∅. 123 956 Annals of Operations Research (2022) 312:949–971 Since ηi ri≤ηj r0+rj,wehaveforclassithat δi>1−ηi ri 1−ρ∗ 0 ≥1−ηi rir0+rjδj+riδi r0+rj+ri−ηj−ηi ⇐⇒ δi(r0+rj+ri−ηj−ηi) >1−ηi ri(r0+rjδj+riδi) ⇐⇒ δi(r0+rj−ηj)>1−ηi ri(r0+rjδj) ⇐⇒ δi>1−ηi ri 1−ηj r0+rj r0+rjδj r0+rj ≥r0+rjδj r0+rj ⇐⇒ δi>δ j+r0(1−δj) r0+rj >δ j, which is in contradiction with δi≤δj. Let Tidenote the sojourn time of jobs of class i.Wenowprovideaninterestingresult related to the sojourn time of jobs. Proposition 5 If δi≤δj. Then, ciE[Ti]≤cjE[Tj]. Proof We know that the sojourn time of jobs of class iand of class jfollow an exponential distribution with rate 1 ri(1−ρ∗ i)and 1 rj(1−ρ∗ j)respectively. Therefore, ciE(Ti)=ci ri(1−ρ∗ i) =ci riδi 1 1−ρ∗ 0 =c0 r0ci ri 1 1−ρ∗ 0 ≤c0 r0cj rj 1 1−ρ∗ 0 =cjE(Tj). And the desired result follows.  4.1 Charaterization of the solution of (LOAD-OPT)withC=2 We now focus on the case C=2. Throughout this article, we refer to this case as the M model. Without loss of generality, we assume that δ1≤δ2. The goal of this section is to fully characterize the solution of (LOAD-OPT) with C=2. We first note that, from Proposition 3, it can never be given the following cases: (i) C0={1}and Cd={2}and (ii) C0={1}and Cb={2}. For the remaining cases, we have the following options: 1. Cd={1,2}. In this case, each class sends all its traffic to the dedicated server. Therefore, ρ∗ i=ηi rifor i=1,2andρ∗ 0=0. According to Proposition 2this occurs when δi≤1−ηi ri,i=1,2. 123 Annals of Operations Research (2022) 312:949–971 957 2. Cd={1}and Cb={2}. In this case, all the traffic of class 1 is sent to Server 1 and the traffic of class 2 is sent to Server 0 and Server 2. As a result, ρ∗ 1=η1 r1and, from (10)and(11) we obtain that ρ∗ 2=1−δ2r0+r2−η2 r0+δ2r2and ρ∗ 0=1−r0+r2−η2 r0+δ2r2. According to Proposition 2, this case occurs when δ1≤1−η1 r1 1−ρ∗ 0 and 1 1−ρ∗ 0 >δ 2>1−η2 r2 1−ρ∗ 0 , i.e., δ1≤1−η1 r1r0+δ2r2 r0+r2−η2 and r0+δ2r2 r0+r2−η2 >δ 2>1−η2 r2r0+δ2r2 r0+r2−η2 , which simplifying gives δ1≤1−η1 r1r0+δ2r2 r0+r2−η2 and 1 1−η2 r0 >δ 2>1−η2 r2 . 3. Cd={1}and C0={2}. In this case, all the traffic of class 1 is sent to Server 1 and the traffic of class 2 is sent to Server 0. As a result, ρ∗ 1=η1 r1,ρ∗ 0=η2 r0and ρ∗ 2=0. According to Proposition 2, this case occurs when δ1≤1−η1 r1 1−ρ∗ 0 and δ2≥1 1−ρ∗ 0 ,i.e. δ1≤1−η1 r1 1−η2 r0 and δ2≥1 1−η2 r0 . 4. C0={1,2}. In this case, the traffic of both classes is sent to Server 0. Hence, ρ∗ i=0for i=1,2 and from (10)thatρ∗ 0=1−r0−η1−η2 r0=η1+η2 r0.According to Proposition 2and using that δ1≤δ2, this case occurs when δ1≥1 1−ρ∗ 0 ,i.e., δ1≥r0 r0−η1−η2 . 5. Cb={1,2}. In this case, the traffic of class iis sent to Server 0 and Server i,for i=1,2. From (10)and(11), it results that ρ∗ 0=1−r0+r1+r2−η1−η2 r0+δ1r1+δ2r2and, for i=1,2, ρ∗ i=1−δir0+r1+r2−η1−η2 r0+δ1r1+δ2r2. Moreover, we conclude from Proposition 2that this occurs when, for i=1,2, 1 1−ρ∗ 0 >δ i>1−ηi ri 1−ρ∗ 0 , which using that δ2≥δ1gives δ1>1−η1 r1r0+δ1r1+δ2r2 r0+r1+r2−η1−η2 and r0+δ1r1+δ2r2 r0+r1+r2−η1−η2 >δ 2>1−η2 r2r0+δ1r1+δ2r2 r0+r1+r2−η1−η2 . We simplify the above expressions and we obtain δ1>1−η1 r1r0+δ2r2 r0+r2−η2 and r0+δ1r1 r0+r1−η1−η2 >δ 2>1−η2 r2r0+δ1r1 r0+r1−η1 . 6. Cb={1}and Cd={2}. We observe that this case is symmetric to the case 2 (where Cb={2}and Cd={1}) and using the same arguments, we get the following conditions 1 1−η1 r0 >δ 1>1−η1 r1 and 1−η2 r2δ2≤r0+δ1r1 r0+r1−η1 . 123 964 Annals of Operations Research (2022) 312:949–971 Table 1 Parameters of the system considered in Sect. 6.3 j=0j=1j=2j=3j=4j=5 rj25 34 38 39 13 27 cj5 10 11 6 86 14 ηj28 33 20 12 24 Fig. 6 Convergence of the fixed-point algorithm presented in Sect. 4.2. The x-axis represents the iterations of the algorithm and the y-axis the load of each server whereas in the bottom line the loads of Server 3, Server 4 and Server 5. We observe that, when the algorithm converges, the load of Server 4 is zero, which means that, for this case, ρ∗ 4=0. We also see that, for Server 3, the initial load in the scenario that is represented by the solid line (that is, the scenario where all the classes send all the traffic to its dedicated server) equals to the load when the algorithm converges. This means that, for class 3, we have that ρ∗ 3=η3 r3. It is important to remark that, as we can also observe in Fig. 6, the algorithm converges to the same values for the three different initial partitions under consideration. We have also started the system with other initial partitions and the obtained results confirmed that the 123 Annals of Operations Research (2022) 312:949–971 965 Fig. 7 Minimum and maximum number of iterations until convergence for systems of different size algorithm always converges. Another interesting property of this algorithm is that the number of iterations required to reach the convergence is very small. Indeed, when the initial partition of Kis such that all the classes belong to Cd, the algorithm converges after 12 iterations. Moreover, for the rest of the cases, the algorithm converges for a less number of iterations. When the initial partition of Kis such that classes 1, 3 and 4 belong to Cdand classes 2 and 5toCb, it converges after 2 iterations and when the initial partition of Kis such that classes 1, 2, 3 and 4 belong to Cdandclass5toCb, it converges after 3 iterations. We now present further numerical work we have performed to analyze the convergence of Algorithm 1for larger systems. For this set of experiments, we consider that the number of dedicated servers, C, varies from 10 to 200 with step 10. For each case we run our algorithm 10 times where, in each run, the parameters of the system are randomly chosen (but satisfying the stability condition); the results are depicted in Fig. 7, where the blue bars represent the minimum number of iterations required for convergence and the yellow bars the difference up to the maximum. The main conclusions of these experiments are twofold: first, we observe that the algorithm converges in all the cases; and second, that the number of iterations required to converge varies between 2 and 6 in all the cases. This means that the convergence of this algorithm is very fast even for large systems with 200 dedicated servers. 7 Conclusions We study the optimal Bernoulli routing in a system with Cdedicated servers and a single multi-skilled server. We first provide a necessary and sufficient condition for the stability of the system. We then reformulate this problem as a optimization problem in terms of the loads of the system and we show the equivalence of both problems. We provide structural properties 123 966 Annals of Operations Research (2022) 312:949–971 on the solution of the derived problem, which allows us to fully characterize the optimal loads of the system when C=2 and also to present a fixed point algorithm whose convergence ensure that the optimal loads are obtained. We compare the performance of this system with optimal loads with a system where all the servers are multi-skilled and also with a system where all the servers are dedicated. Finally, we explore numerically the convergence of the fixed point algorithm and show that, in all the considered cases, the algorithm converges in a very few number of steps. For future work, we are interested in generalizing the results of this article to systems with a more complex topology. Besides, we think that an interesting extension of the performance analysis of this work would be to consider other popular load balancing policies such as Power of Two and Join the Shortest Queue. Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. A Proof of Proposition 1 We first show that if there exists a subset A⊂Ksuch that i∈Aηi>r0+i∈Ari,then the system is not stable. We know from (1)and(2)that r0ρ0+ i∈A riρi= i∈K ηi(1−pi)+ i∈A ηipi = i∈A ηi+ i∈K\A ηi(1−pi) ≥ i∈A ηi >r0+ i∈A ri. Therefore, we have obtained that r0ρ0+i∈Ariρi>r0+i∈Ari, which requires that, at least, the load of one server is larger than one, i.e., that the system is not stable. Let >0 small. We now show that, if (6) holds, then the system is stable. For this purpose, we define the following routing strategy: for all i∈Ksuch that ηi<ri,pi=1−and for all i∈Ksuch that ηi≥ri,pi=ri ηi(1−). For this choice, it is clear that ρj<1forall j∈K. We now focus on Server 0 and we aim to show that  i∈K ηi(1−pi)<r0. 123 Annals of Operations Research (2022) 312:949–971 967 We denote by K∗the set of classes such that ηi≥ri. Hence, the above expression is satisfied if and only if  i∈K\K∗ ηi+ i∈K∗ ηi(1−pi)<r0 Hence,  i∈K\K∗ ηi+ i∈K∗ ηi1−ri ηi (1−)<r0⇐⇒  i∈K\K∗ ηi+ i∈K∗ (ηi−ri(1−))<r0⇐⇒  i∈K\K∗ ηi+ i∈K∗ (ηi−ri+ri)<r0 We know from (6)thati∈K∗(ηi−ri)<r0and, therefore, the above inequality is satisfied if and only if ⎛ ⎝ i∈K\K∗ ηi+ i∈K∗ ri⎞ ⎠<r0− i∈K∗ (ηi−ri). In other words, the desired result follows if we choose >0 such that < r0−i∈K∗(ηi−ri) i∈K\K∗ηi+i∈K∗ri. B Proof of Proposition 2 In the following result, we provide a property that will be useful to show the result of Proposition 2. Lemma 3 Let ˜ C⊆K. Then,  i∈˜ C ηi=r0ρ0+ j∈˜ C rjρj⇐⇒ Cb∪C0⊆˜ C. Proof To simplify the notation, we write D=Cb∪C0.IfD=∅,thenρ0=0andρj=ηj/rj for j=1,2,...,C, which implies clearly that i∈Aηi=j∈SAρjrjfor all A⊂K. We now focus on the case D=∅. We know that ρj<η j/rj,∀j∈Dand ρj=ηj/rj,∀j∈K\D.(14) We now observe that K=D(K\D)and therefore from (7)  i∈D ηi+ i∈K\D ηi=r0ρ0+ i∈D riρi+ i∈K\D riρi. From ρj=ηj/rj,∀j∈K\D, it follows that  i∈D ηi+ i∈A ηi=r0ρ0+ j∈D rjρj+ j∈A rjρj,∀A⊆K\D.(15) 123 968 Annals of Operations Research (2022) 312:949–971 Therefore, for any ˜ C⊆Ksuch that D⊆˜ Cthe constraint (9) is satisfied as an equality. Besides, we now show that for any subset that does not contain D, the constraint (9)is satisfied as an inequality. For all B⊂D,(15) can be written as follows:  i∈B ηi+ i∈D\B ηi+ i∈A ηi=r0ρ0+ j∈B rjρj+ j∈D\B rjρj+ j∈A rjρj,∀A⊆K\D, which, by ρj<η j/rj∀j∈D, gives that  i∈B ηi+ i∈A ηi=r0ρ0+ j∈B rjρj+ j∈D\B rjρj− i∈D\B ηi + j∈A rjρj <r0ρ0+ j∈B rjρj+ j∈A rjρj,∀A⊆K\D. And the desired result follows.  We now prove the result of Proposition 2. Proof The Lagrangian corresponding to (LOAD-OPT)is L(ρ,ν,ζ,γ,ξ)=c0ρ0 1−ρ0 + C  j=1 cjρj 1−ρj + C  j=0 νj(−ρj)+ C  j=0 ζj(ρj−1) +γC  i=1 ηi−r0ρ0− C  j=1 rjρj + A⊂K A=∅ ξA i∈A ηi−r0ρ0− j∈A rjρj Given that the optimization problem is convex, ρ∗,ν∗,ζ∗,γ∗,ξ∗is a solution of (LOAD- -OPT) if it satisfies Karush-Kuhn-Tucker conditions: 0≤ρ∗ j≤1,∀j=0,1,...,C(16) c0 (1−ρ∗ 0)2−ν∗ 0+ζ∗ 0−r0γ∗+ A⊂K A=∅ ξ∗ A=0 (17) cj (1−ρ∗ j)2−ν∗ j+ζ∗ j−rjγ∗+ A⊂K j∈A ξ∗ A=0,∀j=1,2,...,C(18) ν∗ j≥0,ζ ∗ j≥0,γ ∗∈R,ξ ∗ A≥0,∀j=0,1,...,C,∀A⊂K,A=∅ (19) ν∗ jρ∗ j=0,ζ ∗ j(ρ∗ j−1)=0,∀j=0,1,...,C(20) C  i=1 ηi=r0ρ∗ 0+ C  j=1 rjρ∗ j(21)  i∈A ηi≤r0ρ∗ 0+ j∈A rjρ∗ j,∀A⊂K,A=∅ (22) 123 Annals of Operations Research (2022) 312:949–971 969 ξ∗ A i∈A ηi−r0ρ∗ 0− j∈A rjρ∗ j=0,∀A⊂K,A=∅.(23) We observe that the objective function tends to infinity when ρj→1, which implies that ρ∗ j<1,∀j=0,1,...,Cand, as a consequence of this and from (20), ζ∗ j=0,∀j= 0,1,...,C. Furthermore, from Lemma 3and (23), we conclude that ∀A∈Kthat does not contain Cb∪C0, its multiplier verifies that ξ∗ A=0, because for those subsets the constraint (9) is satisfied as an inequality. For Server 0, we know that ρ∗ 0=0ifC0∪Cb=∅and ρ0>0 otherwise. This clearly implies that ν∗ 0≥0ifC0∪Cb=∅and ν∗ 0=0 otherwise. For all j∈K, we know that ρ∗ j=ηj/rjif j∈Cd, whereas ρ∗ j<η j/rjotherwise. This clearly implies that ν∗ j=0 if j∈Cd. We also know that ν∗ j=0if j∈Cbbecause, in this case, ρ∗ j>0, whereas if j∈C0,wehavethatν∗ j≥0. We first prove this result when Cb∪C0=∅. For this case, the load of Server 0 is zero and, thus, it is enough to show that δj<1−ηj/rj.From(17)and(18), we get that c0−ν∗ 0−r0γ∗+ A⊂K A=∅ ξ∗ A=0(24) cj (1−ηj/rj)2−rjγ∗+ A⊂K j∈A ξ∗ A=0,∀j=1,2,...,C.(25) From (24) and since ν∗ 0≥0, it results that ν∗ 0=c0−r0γ∗+ A⊂K A=∅ ξ∗ A≥0⇐⇒ c0 r0 ≥γ∗+ A⊂K A=∅ ξ∗ A, From (25), we obtain that γ∗+ A⊂K j∈A ξ∗ A=cj rj 1 (1−ηj/rj)2,∀j=1,2,...,C. Therefore, cj rj 1 (1−ηj/rj)2=γ∗+ A⊂K j∈A ξ∗ A≤γ∗+ A⊂K A=∅ ξ∗ A≤c0 r0 , which gives the desired condition, i.e., δj≤1−ηj/rj,∀j=1,2,...,C. We focus on the case C0=∅or Cb=∅. We note that (17)and(18) can be written as follows: c0 (1−ρ∗ 0)2−r0γ∗+ A⊂K Cb∪C0⊆A ξ∗ A=0 (26) cj−ν∗ j−rjγ∗+ A⊂K Cb∪C0⊆A ξ∗ A=0,∀j∈C0(27) cj (1−ηj/rj)2−rjγ∗+ A⊂K Cb∪C0⊆A j∈A ξ∗ A=0,∀j∈Cd.(28) 123 970 Annals of Operations Research (2022) 312:949–971 cj (1−ρ∗ j)2−rjγ∗+ A⊂K Cb∪C0⊆A j∈A ξ∗ A=0,∀j∈Cb.(29) We aim to show that, for all j∈C0,δj≥1 1−ρ∗ 0 ,andforall j∈Cb1 1−ρ∗ 0 >δ j>1−ηj/rj 1−ρ∗ 0 . For the first condition, we observe that from (26)and(27), it follows that, for all j∈C0, γ∗+ A⊂K Cb⊆A ξ∗ A=c0 r0 1 (1−ρ∗ 0)2=cj rj −ν∗ j rj , which, using that ν∗ j≥0, gives that δj≥1 1−ρ∗ 0 . We now show the second condition, i.e., 1 1−ρ∗ 0 >δ j>1−ηj/rj 1−ρ∗ 0 for all j∈Cb.From(26) and (29), it follows that, for all j∈Cb, γ∗+ A⊂K Cb∪C0⊆A ξ∗ A=c0 r0 1 (1−ρ∗ 0)2=cj rj 1 (1−ρ∗ j)2(30) ⇐⇒ 0=cj rj 1 (1−ρ∗ j)2−c0 r0 1 (1−ρ∗ 0)2(31) ⇐⇒ δj=1−ρ∗ j 1−ρ∗ 0 ,(32) which gives that ρ∗ j=1−δj(1−ρ∗ 0), as desired. Using the last expression and that, for j∈Cb,0<ρ ∗ j<η j/rj, the desired result follows, i.e., 1 1−ρ∗ 0 >δ j=1−ρ∗ j 1−ρ∗ 0 >1−ηj/rj 1−ρ∗ 0 . To finish, we compute the loads of all the servers. First, for j∈Cd,we have clearly that ρ∗ j=ηj rj. Besides, we use that for all j∈Cb,ρ∗ j=1−δj(1−ρ∗ 0), and from the expression (7), it follows that C  i=1 ηi=r0ρ∗ 0+ C  j=1 rjρ∗ j⇐⇒ C  i=1 ηi=r0ρ∗ 0+ j∈Cb rj(1−δj(1−ρ∗ 0)) + i∈Cd ηi. And rearranging both sides of the above expression, we obtain that ρ∗ 0=1−r0+j∈Cbrj−i∈Cb∪C0ηi r0+j∈Cbδjrj . 123 Annals of Operations Research (2022) 312:949–971 971 And the desired result follows.  References Altman, E., Ayesta, U., & Prabhu, B. J. (2011). Load balancing in processor sharing systems. Telecommunication Systems, 47(1–2), 35–48. Bonald, T., Comte, C., & Mathieu, F. (2017). Performance of balanced fairness in resource pools: A recursive approach. Proceedings of the ACM on Measurement and Analysis of Computing Systems.https://doi.org/ 10.1145/3154500. Chen, J., Do, J., & Shi, P. (2020). A survey on skill-based routing with applications to service operations management. Queueing Systems, 96, 53–82. Feng, H., Misra, V., & Rubenstein, D. (2005). Optimal state-free, size-aware dispatching for heterogeneous M/G/-type systems. Performance Evaluation, 62(1–4), 475–492. Fourneau, J.-M. (2020). Modeling green data-centers and jobs balancing with energy packet networks and interrupted Poisson energy arrivals. SN Computer Science, 1(1), 28. Gardner, K., Zbarsky, S., Doroudi, S., Harchol-Balter, M., & Hyytia, E. (2015). Reducing latency via redundant requests: exact analysis. In Proceedings of the 2015 ACM SIGMETRICS International conference on measurement and modeling of computer systems, SIGMETRICS ’15, Association for Computing Machinery, New York, NY, USA (pp. 347–360). Gaujal, B., Hyon, E., & Jean-Marie, A. (2006). Optimal routing in two parallel queues with exponential service times. Discrete Event Dynamic Systems, 16(1), 71–107. https://doi.org/10.1007/s10626-006-6179-3. Graham, C. (2000). Chaoticity on path space for a queueing network with selection of the shortest queue among several. Journal of Applied Probability, 37(1), 198–211. https://doi.org/10.1239/jap/1014842277. Harchol-Balter, M. (2000). Task assignment with unknown duration. In Proceedings 20th IEEE international conference on distributed computing systems (pp. 214–224). IEEE. Harchol-Balter, M., Crovella, M. E., & Murta, C. D. (1999). On choosing a task assignment policy for a distributed server system. Journal of Parallel and Distributed Computing, 59(2), 204–228. Kelly, F. (1979). Reversibility and stochastic networks. Technical report. Koole, G., Pot, A., & Talim, J. (2003). Routing heuristics for multi-skill call centers,2, 1813–1816. Liu, Z., Lin, M., Wierman, A., Low, S., & Andrew, L. L. H. (2015). Greening geographical load balancing. IEEE/ACM Transactions on Networking, 23(2), 657–671. Lu, Y., Xie, Q., Kliot, G., Geller, A., Larus, J. R., & Greenberg, A. (2011). Join-idle-queue: A novel load balancing algorithm for dynamically scalable web services. Performance Evaluation, 68(11), 1056– 1071. Mitzenmacher, M. (2001). The power of two choices in randomized load balancing. IEEE Transactions on Parallel and Distributed Systems, 12(10), 1094–1104. Taboada, I., Aalto, S., Lassila, P., & Liberal, F. (2017). Delay and energy-aware load balancing in ultra-dense heterogeneous 5G networks. Transactions on Emerging Telecommunications Technologies, 28(9), e3170. Vvedenskaya, N. D., Dobrushin, R. L., & Karpelevich, F. I. (1996). Queueing system with selection of the shortest of two queues: An asymptotic approach. Problems of Information Transmission, 32(1), 15–27. Wallace, R. B., & Whitt, W. (2005). A staffing algorithm for call centers with skill-based Routing. Manufacturing & Service Operations Management, 7(4), 276–294. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123