Energy-aware optimization of scalable load balancing strategies by Rajita Chandak B01090236 Faculty Advisor: Dr. Kavita Ramanan Post-doctoral Advisor: Dr. Debankur Mukherjee Second Faculty Reader: Dr. Matt Harrison Department of Applied Mathematics Brown University April 2019 Abstract Acknowledgements I would like to thank Dr. Kavita Ramanan for her support throughout this process. I would also like to thank Dr. Debankur Mukherjee for his guidance and mentorship. Abstract Queuing systems are one of the most prominent technological evolution in today’s digital world. Data centers and cloud networks operated by large companies like Microsoft, Google and Amazon rely on large-scale queuing systems for their operations. The goal of minimiz- ing energy consumption while optimizing performance has been the focus of recent literature around the subject. In this thesis, we consider the TABS scheme and investigate through simulations and a theoretical proof the implications of various parameters of the scheme on a variant of the decentralized parallel queuing system. Simulation results explore how changing the assumption of exponential service time distributions affects various performance metrics. Our theoretical result builds on previous results in order to identify patterns in fluctuations around mean behaviour of the system over long time intervals. Contents 1 Introduction 1 1.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3 Key Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Outline of Paper . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Details of Model and Algorithm 5 2.1 TABS Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Notation and Previous Results . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3 Detailed Outline of Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3 Simulation Results 7 3.1 Queue Lengths with Varying Number of Servers . . . . . . . . . . . . . . . . 7 3.2 Fraction of Active Servers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.3 Energy and Wait-Time Graphs . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.4 Folded Normal Service Time Distribution . . . . . . . . . . . . . . . . . . . . 13 3.5 Hyperexponential Service Time Distribution . . . . . . . . . . . . . . . . . . 14 3.6 Comparing Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.6.1 Comparing Different Setup Times . . . . . . . . . . . . . . . . . . . . 16 3.6.2 Comparing Number of Servers . . . . . . . . . . . . . . . . . . . . . . 18 3.6.3 Comparing Variance Effect . . . . . . . . . . . . . . . . . . . . . . . . 19 4 Limit Theorem Proof 21 4.1 System Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.2 Martingale Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.3 Martingale Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.4 Martingale Behaviour . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.4.1 Quadratic Variation . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.4.2 Quadratic Covariation . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.5 Martingale Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.6 Martingale Limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 CONTENTS CONTENTS 5 Conclusion 34 5.1 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.2 Open Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 References 36 1 INTRODUCTION 1 Introduction 1.1 Background and Motivation As we move towards an increasingly digitized world, concerns of energy consumption have become more prominent. Modern day data centers and cloud networks are just one such technological evolution that present financial and environmental challenges. It is estimated that US data centers consume about 70 million MegaWatt hours annually, the equivalent of 6 million homes [7]. A crucial consideration for data centers is to achieve operational efficiency in terms of user-perceived performance while managing energy consumption and efficient use of servers throughout an extended period of time. A big challenge with this set up is that client demand is uncertain and often changes drastically with time [2]. Thus we may dynamically scale the amount of resources that are actively utilized with the actual observed demand conditions in order to minimize energy consumption and operation costs, while achieving targeted performance goals. This methodology of scaling resources is referred to as service elasticity. Achieving ideal service elasticity is challenging because by turning servers off during periods of low demand, data centers incur additional additional setup costs in terms of the time required to turn on servers, producing a time lag in the system. Typically the setup time for servers is orders of magnitude higher than the time required to service tasks [5]. This time lag can be counteracted by maintaining a significant number of idle servers on in case of sudden demand spikes. Another scheme to balance energy consumption with data center performance is auto-scaling. Auto-scaling adjusts capacity as a response to changing demand. It is currently a widely employed practice in industry by companies like Facebook, Google and Microsoft [8]. Many auto-scaling approaches use historical information to forecast load predictions and manage capacity, information that may not always be available in data centers. These schemes often look at centralized queueing systems. In practice, data centers do not maintain a centralized queue. Instead, they distribute demand across parallel servers, each with its own queue. Operating within a parallel queuing network system provides less global information. Since parallel queuing systems allow for the chance that some servers remain idle while tasks are waiting at other servers, expected performance is not as easily analyzable as for the centralized queuing system and thus re- quires that algorithms to manage energy consumption operate independent of historical load prediction or global queuing information. 1.2 Literature Review Load-balancing schemes for large-scale systems found in cloud networks and data centers have been studied extensively in the past. In particular, stability and energy consumption of systems with centralized queues is well understood [4, 10]. Many papers consider the idea of putting servers to sleep in the centralized system because servers in sleep mode consume significantly less energy than idle servers. It is estimated that idle servers consume 60-70% of the energy used during peak demand [2, 5]. However, as the number of servers becomes large, 1 1.2 Literature Review 1 INTRODUCTION maintaining a centralized queue becomes complicated to handle and presents inefficiency in operation speed [14]. Stankovic [14] shows in his paper that three of the decentralized algorithms he tested for queuing systems show stability and improved performance with only a modest increase in costs. A decentralized system maintains servers with individual, parallel queues. This is an example of a load balancing algorithm that is used to improve performance. That is, the system distributes demand (the network load) efficiently across servers. Such systems have been the focus of more recent research on cloud networks [12, 11]. Furthermore, in typical data center architectures and cloud environments no centralized queue is maintained, and load balancing algorithms immediately distribute incoming tasks among parallel queues [1, 12]. It is estimated that servers account for about 45% of a data center’s amortized costs and utilities (electrical) make up about 15% of operational costs [2]. Thus, improving perfor- mance and lowering energy consumption of servers is of interest to companies that manage these centers. Since data centers, operating with tens of thousands of servers, face large amortized costs, operational efficiency is another prime concern for businesses. Unfortu- nately data centers tend to operate at low operational efficiency-usually less than 10%, either due to server memory constraints or safeguarding against unexpected demand shifts [2, 8]. That is, data centers tend to use significantly fewer number of servers than available and as a result a large portion of energy costs are due to idle servers (otherwise referred to as servers in standby) maintained through the day. One main reason for this, as stated by Greenberg et al. in [2] is uncertainty in expected demand. Since the demand for data centers can spike unexpectedly, businesses tend to err on the side of caution by operating with more than the necessary number of servers. This allows data centers to maintain higher capacity in case of a spike in demand and reduce the risk of system failure during times when profits can be maximized. There are many ways that costs could be reduced, as explored in [2]. One of the most impactful ones is reduced power consumption. Turning idle servers off to reduce energy consumption has been one of the most popular suggestions in literature. However, this comes at a cost. Data suggests the average setup time (the amount of time required for a server to turn on, warm up and be prepared to service new tasks) for a server is about 200 seconds whereas the time required to service a task is typically less than 1 second [4], meaning unexpected increases in demand could lead to system failure despite data centers having the capacity to keep up with the demand. Additionally, data centers could implement auto-scaling techniques. Auto-scaling refers to an automatic system of adjusting capacity of data centers to meet demand while maintaining performance goals. In much of the literature relevant to this thesis [12, 11, 4], auto-scaling is implemented in the form of turning servers off during periods of low demand. Thus, investigation into decentralized queuing systems presents great potential for reduced energy consumption, lower costs and higher operational efficiency, if managed carefully [2]. A particularly relevant idealized model for centralized load management is considered by Gandhi et al. in [4]. Specifically, they consider a variant of a standard queueing system, referred to as M/M/N/setup/delayedoff. In this scheme, servers can be turned on and off as needed. More specifically, in the proposed algorithm when a server finishes servicing a task, and finds no additional tasks in its queue, it waits for an exponential distributed amount of time with expectation µ. If a task arrives during this time then it is immediately 2 1.2 Literature Review 1 INTRODUCTION assigned to the server (or one of the other servers on standby at random), otherwise the server is turned off. When a task arrives, if there is no server in standby, then the task waits in the centralized queue. The server takes an exponentially distributed amount of time, with expectation ν to setup and become idle once the time available to service tasks in the queue has elapsed. . Gandhi et al.[4] present in their paper a closed form solution for the M/M/N/setup/delayedoff system, through what they term as a Recursive Renewal Reward method. The methodology combines features of Renewal Reward theory and analysis of high demand intervals to find closed form solutions for the system’s power consumption. The method introduce in this paper is robust in its application to Markov Chains. The idea of turning idle servers off is adapted to a decentralized system and introduced by Mukherjee et. al. in [12] as the TABS scheme. The TABS scheme was first proposed by Mukherjee et. al. [12] as a scheme that aims to increase efficiency and reduce energy consumption on average for large scale distributed systems with parallel queues. The main goal of the paper was to dynamically scale the number of servers that are actively utilized with the actual observed task load conditions to curtail cost and energy consumption, while still satisfying the target performance criteria. By looking at a joint auto-scaling and load balancing scheme (in the form of turning servers off when idle for an exponentially distributed time and turning them back on as needed) which does not require any knowledge of the global queue length or explicit information of system parameters the scheme is still able to achieve near-optimal service elasticity. The TABS scheme allows for servers to turn themselves off after staying idle for an exponentially distributed amount of time (detailed in Section 2.1 of this thesis). This scheme is shown to reduce average energy consumption that would have otherwise persisted were servers to remain in standby mode for long periods of time. If all servers are busy when a new task arrives, the task joins the queue of a randomly selected busy server and the dispatcher selects one of the off servers to begin setup. Mukherjee et. al. [12] also considers how not operating at full server capacity may affect average wait time of the tasks in the system. The results show that the TABS scheme provides similar waiting times on average when compared to the M/M/N/setup/delayedoff schem despite the former being a centralized scheme. Mukherjee et al. prove that both the waiting time of tasks and the relative energy portion consumed by idle servers vanish in the limit as the standby time µ tends to infinity. While [12] assumes finite buffer capacity, the question of system stability for an infinite buffer capacity system under the TABS algorithm is addressed in [11] by Mukherjee and Stolyar through a proof using induction and weak-monotonicity. Mukherjee and Stolyar [11] examine the asymptotic behaviour of the system as the number of servers becomes large and use mean-field analysis to prove stability and convergence of the TABS scheme (in Sections 2 and 3 of [11]). Important parameters in both [12] and [11] are the expected setup time, denoted by ν, and the expected standby time, represented by µ. From the results of [2], we may expect that the average amount of time for which servers remain idle before turning off or the expected amount of time it takes for servers to turn on and warm up before servicing tasks may affect the performance of the system. However, results of both papers [11] and [12] show that as the number of servers increases towards infinity, system performance metrics like the average waiting time and energy wastage do not depend on ν and µ values by virtue of the fluid-limit scaling techniques used. The parameters ν and µ are scaled according to 3 1.3 Key Contributions 1 INTRODUCTION the number of servers and thus the steady-state of the queuing system is independent of their values. In this thesis, we look at finer details of the system and their dependence on ν and µ. Specifically, we will look at a diffusion scaling of the system. By looking at the behaviour of the TABS scheme in this asymptotic regime, we try to understand the steady state of the system under a different scaling method. The goal of looking at a diffusive scaling for the system is to eventually understand how ν and µ affect performance metrics like average task wait time and energy wastage. This paper consists of two different analyses of the system to better understand it’s behaviour in the limit. First, we simulate the system with a large number of servers (ranging from 100-1000 servers) under different ν and µ values, with varying service time distributions. The simulation results show how the energy consumption and average task wait time varies with the service time distribution parameters and ν and µ in the long run. In this these, we also derive a proof for the limiting behaviour of this system. The proof extends ideas presented by Eschenfeldt and Gamarnik in [3] done on a similar queuing system. By understanding the steady state convergence of the diffusion scaling of the system under the TABS scheme, we aim to provide finer guidelines for building a more efficient system. 1.3 Key Contributions In this thesis, we provide provide a more refined study of the performance of the TABS algorithm. The first half of this thesis presents a simulation-based analysis of the effect of varying service time distributions on the stability and behaviour of the system. By looking at the effect of several distributions centered around expectation 1 with different variances on the system’s energy consumption and average task wait time, this paper explores robustness of the system with respect to the service time distributions. In the second half of this thesis, by looking at the diffusion scaling of the TABS system, we are able to understand the effect of changing expected standby time and setup time on the long term behaviour of the system. The limit theorem proof provides a finer scale analysis of the system than previously examined in literature pertaining to the TABS scheme. The proof sets up a martingale representation of the system, with appropriate scaling after which the limit of the system of equations converges to a known limit, proving stability and convergence of the system with respect to varying setup and standby time time for servers. 1.4 Outline of Paper The remainder of the paper is organized as follows. In section 2 we describe the system set up, the TABS Algorithm, and some key results from previous literature the provide support for following sections. In section 3 we present the simulations of the system, comparing various different system states. We also discuss the consequences of observed simulation results in this section. In Section 4 we present the limit theorem proof of convergence. Finally, in Section 5 we provide some concluding remarks and potential extensions of the work presented in this thesis. 4 2 DETAILS OF MODEL AND ALGORITHM 2 Details of Model and Algorithm 2.1 TABS Algorithm The queuing system we are looking at consists of 2 key actors. The first is the dispatcher that manages all incoming tasks. The second actor is the collection of servers. Servers operate based on task assignment from the dispatcher. The TABS algorithm uses the idea of tokens that are accessible to the dispatcher in order to assign new tasks to a server and turn servers on. Busy servers are identified by ‘yellow’ tokens, idle servers are identified by ‘green’ tokens, off servers are identified by ‘red’ tokens and servers in setup mode are identified by ‘orange’ tokens. Thus, in this system the only global knowledge available is the number of servers that are busy, in standby, in setup mode or off. No information about each server’s queue length is stored. The scheme is reproduced below is outlined in Section 2 of [12]: • When a server becomes idle, it sends a green token to the dispatcher, waits for an Exp(µ) time ( µ denotes expected standby time), and turns itself off by sending a red token to the dispatcher (this replaces the initial green token corresponding to the server) if it does not receive a task in the mean time. • When a task arrives, the dispatcher selects a green token at random, if there are any, and assigns the task to that server (the corresponding green token is replaced by a yellow token). Otherwise, the task is assigned to an arbitrary busy server, and if at that arrival epoch there is a red token at the dispatcher, the dispatcher selects one at random, and the setup procedure of the corresponding server is initiated (setup time distributed as Exp(ν)), replacing its red token by an orange token. If all servers are off when the task arrives, it is thrown away by the dispatcher. We consider this to mean the system has lost a task. • After the server in setup mode is ready to service tasks, it replaces the orange token with a green token to the dispatcher and waits for an Exp(µ) time for a possible assignment of a task, and turns itself off by replacing the green token with a red one if no task is assigned in the interim period. A key assumption in literature for previous results underlying data center simulation as well as the TABS algorithm is that tasks arrive with an exponential distribution with some expectation λ < 1 and tasks are serviced independently by servers with an exponential distribution with expectation 1. Previous literature alludes to additional constraints that can be imposed on the system such as buffer capacity for server queue lengths. In this case, the dispatcher will all be able to access information of whether each busy server has queue length equal to the buffer capacity, B. When the dispatcher receives a new task, the task will be assigned to a random server that has not reached capacity according to the same algorithm outlined above. In the event that all servers have queue length B, the dispatcher will throw the incoming task away. 5 2.2 Notation and Previous Results 2 DETAILS OF MODEL AND ALGORITHM 2.2 Notation and Previous Results Notation: For a system with n servers, Qni (t) denotes the number of servers with queue length greater than or equal to i at time t, including the task being serviced. ∆n0 (t) denotes the number of off servers at time t and ∆n1 (t) denotes the number of servers in setup mode at time t. Thus, (Qni (t), ∆n0 (t), ∆n1 (t)) provides a Markovian representation of the system since the servers are independent and identically distributed. Furthermore, for any fixed n, this process is an irreducible countable-state Markov chain [12]. Thus, we know the existence of a unique stationary distribution. Let µ denote the average standby time for an idle-on server before it turns off and let ν denote the average setup time required for an off server to turn on and be ready to start servicing tasks. Previous convergence results for the TABS algorithm from [11] look at fluid-scaled lim- iting behaviour of the system, that is, (Qni (t)/n, ∆n0 (t)/n, ∆n1 (t)/n) and its convergence to a limit. Specifically, for fixed µ, ν > 0 and λ < 1, the system converges to the following limit  n Q1 (t) ∆n0 (t) ∆n1 (t)  , , ⇒ (λ, 1 − λ, 0) n n n Qni (t) ⇒ 0 for any i > 1, (2.1) From these results we see that the limiting behaviour of the system under fluid scaling is independent of µ and ν. A finer analysis of the system should provide further insight into the extent to which µ and ν and the specific form of the service time distribution affect the system in the long run. 2.3 Detailed Outline of Goals In this section, we return to the goals outlined in Section 1.3 in the context of the notation and known results established in the previous sections. Thefirst mathematical result  we present Qn1 (t) ∆n0 (t) ∆n1 (t) in Section 3 through simulations is the behaviour of , , in the limit n n n given alternative service time distributions. We maintain the constraints established on the exogenous parameters introduced in [11] and [12]. That is, we set µ, ν > 0 and λ < 1, while also constraining our various service time distributions to have expectation 1. We explore the effects of generalized Pareto, Folded Normal and Hyperexponential distributions on the system’s performance in comparison with the Exponential distribution. The performance metrics used to understand the behaviour of the TABS scheme under these service time distributions will be the same as those in [11]-average task wait time and average energy consumption. In Section 4, we explore the behaviour of the TABS system under diffusion scaling and over long time intervals. That is, we√ center each of the system variables around their fluid limit expectations and then scale by n. This allows for us to understand the behaviour of the system on a finer scale and thus understand the impact of various µ and ν values. In 6 3 SIMULATION RESULTS mathematical notation, we are looking at the long run behaviour of  n Q1 (t) − λn ∆n0 (t) − (1 − λ)n ∆n1 (t)  √ , √ , √ and n n n Qni (t) √ , for all i ≥ 2. n 3 Simulation Results In this section we will explore the effects of alternate service time distributions on the queuing system under the TABS algorithm. We can simulate various features such as queue lengths, fraction of busy, standby and off servers, energy consumption and average task wait time to see how it compares to the TABS results from [11]. All graphs that follow were produced using MATLAB. We will start by first adjusting the service time distribution to be a generalized Pareto distribution with expectation (approximately) 1 to test how this affects long term behaviour of the system. The pdf of the generalized Pareto distribution is given by f  1  −1− 1 (x − θ) k  f (x|k, σ, θ) =  1+k  σ  σ  For the simulations in this thesis we set the parameters to be k = −0.6, σ = 0.7, θ = 0.56. The expectation of the Generalized Pareto distribution is 1 and the variance is 0.0870. The following sections will show similar simulations for Folded Normal and Hyperexponential service time distributions. We will also explore the effect of service time distribution variance on system performance. 3.1 Queue Lengths with Varying Number of Servers Queues with varying arrival rates: Before looking at the performance metrics shown in [12], we will look at the behaviour of the basic system with varying parameters. First, we start with a small number of servers and observe how the queue lengths vary over time based on the arrival rate. For this, we set the standby time to be large enough such that the likelihood of servers turning off during the simulation is very low (µ = 500s). We start the system with all servers idle and on. Figure 1 shows how a 5 server system evolves over 1000 seconds. 7 3.1 Queue Lengths with Varying Number of Servers 3 SIMULATION RESULTS (a) Arrival Rate: 1.11 (b) Arrival Rate: 0.67 Figure 1: Queue lengths, multiple servers As we can see in Figure 1(a), with a high arrival rate, the queue lengths for all servers are growing to infinity. This tells us that when more tasks arrive than are serviced in a given time period, on average, the queue lengths will tend to infinity. This fits into our model’s expectations. Since the service time takes 1 second on average for any given server, if tasks arrive faster than any single server can process them, the queue lengths will grow rapidly and tend to infinity. In Figure 1(b) where the arrival rate is smaller than the expected servicing time, however, we see the system is stable. Servers have queue lengths mostly between 1-3 tasks. The queue lengths do not grow rapidly because the servers are, on average, able to service tasks faster than they arrive. Thus, we can note in systems with arrival rates greater than service rate, the system will be unstable. Queues for reduced standby time: We can now reduce the expected standby time to see how this affects the system. The simulations below use an average standby time of 5 seconds. The average setup time is kept constant at 1s to minimize the effect this will have on the system. All other parameters of the system are unchanged from previous simulations. Thus, the resulting change in queuing behaviour in this new system should reflect the ef- fect of changing the expected standby time for servers. We set up the simulation such that if all servers are off, the task is thrown away by the dispatcher. As mentioned in Section 2.1, in this case, we consider the task to be ‘lost’ by the system. For the sake of visual representation, we maintain a small number of servers in the system. 8 3.1 Queue Lengths with Varying Number of Servers 3 SIMULATION RESULTS (a) Arrival Rate=0.2 (b) Arrival Rate=0.67 (c) Arrival Rate=0.83 Figure 2: Queue lengths, 3 servers, varying arrival rates Note that there are a large number of tasks lost in these simulations. This happens because once a server becomes idle, it will turn off on average after 5 seconds if it does not receive a new task to service. If all servers are off, the system looses the task. Since tasks are lost in the system and servers have an infinite buffer capacity, we know that there are many instances in these simulations where all servers have turned off. We also notice that as the arrival rate increases, the number of tasks lost decreases. This implies that the average standby time must change based on the arrival rate in order to minimize the number of tasks that are lost. If the arrival rate is very small, and we do not want to loose any tasks, servers should have a larger expected standby time. A system with small average standby time for the servers will loose fewer tasks if the arrival rate is relatively high. The smaller the inter-arrival time between new tasks, the smaller chance that an idle server will turn off before the arrival of the next task and thus the greater the chance there will be at least one server available to service tasks. Varying expected setup time: Finally, we can observe how increasing the average server setup time will affect queue lengths and number of tasks lost by the system. Once again, all other parameters of the system are left unchanged. In the simulations shown in Figure 3, the arrival rate is fixed at 0.67 and the average standby time is constant at 5s. As expected, the longer it takes for servers to set up, the more tasks are likely to be lost given all other 9 3.2 Fraction of Active Servers 3 SIMULATION RESULTS parameters of the system remain unchanged. (a) ν = 2s (b) ν= 1s (c) ν = 0.3s Figure 3: Queue lengths, 3 servers, varying average setup times 3.2 Fraction of Active Servers Now that we understand the evolution of the system over time and how each parameter affects system behaviour (in terms of number of tasks lost in the system), We can look at systems with a large number of servers, starting with 100 servers. Plotting queue lengths for every server, as was done earlier, would be inefficient and not very informative. Instead, we look at the proportion of servers that are busy, standby, off or in setup mode. This gives a better picture of how the resources, in terms of capacity, are being used in the system. From the results of [11], we expect that with a large number of servers, the system to converge to the limit given in Equation 2.1. That is, with a large number of servers, in a stable system, the fraction of busy servers will converge to λ, the fraction of idle-off servers converges to 1 − λ and the fraction of idle servers and servers in setup mode converges to 0, where λ is the arrival rate for tasks in the system. 10 3.3 Energy and Wait-Time Graphs 3 SIMULATION RESULTS (a) λ=0.14, ν=20s, (b) λ=0.11, ν=20s, µ=10s µ=20s (c) λ=0.14, ν=20s, µ=10s Figure 4: Server Proportions The graphs in Figure 4 show that the simulations with a generalized Pareto distribution match the theoretical result from [11] for varying arrival rates, setup time and standby times. Thus, the fact that the graphs in Figure 4 show similar patterns to theoretical expectations for exponential distribution, provides a starting point to analyse the TABS scheme under various other distributions and test for robustness of the results in literature [11, 12]. 3.3 Energy and Wait-Time Graphs The efficiency of the queuing system can be determined in two different ways. The first is to look at the average energy consumption of the system for varying standby times. This concerns the operational costs of data centers and aims to minimize them. The second is to 11 3.3 Energy and Wait-Time Graphs 3 SIMULATION RESULTS look at efficiency in terms of how fast tasks can be serviced. This concerns user experience and affects total business engagement by data centers. Data centers want to service tasks as soon as possible—delayed servicing can lead to loss in number of customers and total revenue. The goal is that both the average energy consumption and the average wait time for a new task in the system be as low as possible for quality performance and high operational efficiency. The values in following graphs are obtained using the methodology of calculating average wait time and energy consumption from [12].Average wait time is determined by observing the average of number of servers with queue lengths greater than 1, proportional to arrival rate. Average energy consumption is determined by taking the weighted average of the number of busy servers, servers in setup mode and idle-on servers. The graphs in Figure 5 show average energy consumption and wait-times corresponding to the system evolution plotted in Figure 4. (a) Arrival Rate=0.14, setup time=10s (b) Arrival Rate=0.11, setup time=10s (c) Arrival Rate=0.14, setup time=20s Figure 5: Average Energy Consumption and Wait-time From Figure 4 we see that as the expected standby time increases, the average wait time for a task approaches zero. Also note that as µ increases, average energy consumption levels out. We also observe that as setup time increases, the limit that average energy consump- tion approaches, increases. These results match theoretical predictions and simulation results from [11]. This provides the first step in establishing the hypothesis that TABS convergence results may be independent of the service time distribution given that the expectation of the 12 3.4 Folded Normal Service Time Distribution 3 SIMULATION RESULTS distribution is 1. In the next sections we simulate the behaviour of the TABS scheme under the Folded Normal and Hyperexponential distributions with expectation 1. 3.4 Folded Normal Service Time Distribution Just as was done with the generalized Pareto distribution, we simulate the TABS scheme with a folded normal distribution with expectation 1. The service time samples for each server are taken from the Folded Normal pdf, given by 1 −(x−α)2 f (x) = √ exp( ) σ 2π 2σ 2 Where α denotes the expectation of the distribution. The mean and variance of this pdf are set to equal 1 for all of the results presented below. We find the proportion of servers plots and average energy consumption and wait time. One such simulation result is shown in Figure 6. The parameters for this simulated run are λ = 0.14, ν = 20s, µ = 20s, σ 2 = 1. The folded normal distribution is simulated first from the normal distribution with expectation 1 and variance 1 and then the absolute value of the sample is taken to be the service time for a given server when it begins servicing a task. (a) Server Proportions (b) Average energy consumption and wait time Figure 6: Simulation with folded normal service distribution As we see from Figure 6, the proportion of active servers behaves similarly to the graphs from the generalized Pareto distribution and exponential service distribution. Furthermore, the average energy consumption approaches a limit of about 0.45, which is very similar to the limit observed with the generalized Pareto distribution. The average wait time also vanishes as the expected standby time increases, aligning with the results from [12]. 13 3.5 Hyperexponential Service Time Distribution 3 SIMULATION RESULTS 3.5 Hyperexponential Service Time Distribution Finally, we simulate the system under TABS scheme with a Hyperexponential service time distribution. The Hyperexponential distribution is defined as the weighted average of n independent exponential distributions with expectation λi ∀i = 1, ..., n. The probability density function of a Hyperexponential distribution with n independent exponential random variables (X1 , X2 , ..., Xn ), where each Xi has the pdf of an exponential distribution with mean from the corresponding entry in the mean vector λ = (λ1 , λ2 , ..., λn ). The Pn probability of each Xi is determined by the probability vector p = (p1 , p2 , ..., pn ) s.t. i=1 pi = 1 is given by n X fY (y) = pi fXi (x) i=1 To simulate from the Hyperexponential with expected value of 1 we use the following pseu- docode: • Set n = 3, that is, (X1 , X2 , X3 ) with a mean vector of (1.2, .9, 1). Assign probabilities of (0.1, 0.3, 0.6) to for the respective Xi . (This distribution has expectation 1 and variance approximately equal to 1 ) • Generate a sample, U, from the uniform distribution on [0, 1] • Given, U sample for the Hyperexponential from fY (y|U = u)  1  X 1 ∼ exp(1.2) if u ≤ 0.1 1.2        fY (y|U = u) = X2 ∼ 1 exp(0.9) if 0.1 ≤ u ≤ 0.4    0.9      X ∼ exp(1) 3 if 0.4 ≤ u The final sample generated from the exponential random variable is the service time re- quired to service the next task for the server that just started processing a new task. Figure 7, below, shows a sample evolution of the fraction of active servers and the corresponding expected wait time and average energy consumption. We notice that the graphs appear similar to those seen in previous sections as well as the previous TABS results. The conver- gence limits for average energy consumption and average wait time appear to be identical to the limits observed with generalized Pareto, Folded Normal and exponential service time distributions. 14 3.6 Comparing Distributions 3 SIMULATION RESULTS (a) Server Proportions (b) Average energy consumption and wait time Figure 7: Simulation with Hyperexponential service distribution 3.6 Comparing Distributions So far, we have analysed the affects of various service time distribution on system behaviour individually. These results provide evidence that long term convergence behaviour of TABS remains consistent regardless of the service time distribution. To best understand how the system depends on the various service time distributions, we can plot the average energy consumption and average wait time from simulations with all of the service time distributions together and compare the results. Figure 8 shows the comparisons between the generalized Pareto, Folded Normal, Hyper- exponential and Exponential distributions. The key difference in the plots is the variance of the each of the distributions. The variance for the generalized Pareto distribution is 0.087, for the rest of the distributions the variance is approximately 1. Another key change made in the following simulation plots is that the simulation was run 10 times and each point on the average wait time and energy consumption plots is the average of the 10 simulated runs. The reason for taking the average of 10 simulations is to remove any randomness that may show up in a single run and present the limiting behaviour of the system under this modified TABS algorithm. 15 3.6 Comparing Distributions 3 SIMULATION RESULTS (a) Average Energy Consumption (b) Average wait time Figure 8: Comparing Multiple Service Time Distributions We can see clearly in Figure 8(b) that average wait time for a task is approximately the same for each distribution at each µ value as well as in the limit µ → ∞. The key difference in results appears in 8(a) where the Folded Normal distribution has a strictly higher energy consumption than any other distribution for all µ values. Although, we do note that the general trend of energy consumption remains the same across all distributions. It is also worth noting that although the folded normal seems to consume more energy than any of the other distributions, the difference between the amount of the energy consumption is at most 0.02, which may arguably be a statistically insignificant difference. 3.6.1 Comparing Different Setup Times To ensure this pattern is not unique to the set of parameters chosen for the simulation pre- sented in Figure 8, we run the simulation with various different setup times and compare the average wait time and energy consumption across the different service time distributions. Average Energy Consumption: We notice identical behaviour in energy consumption across distributions and different average setup time in Figure 9. It is noteworthy that the convergence limits for each distribution remain unchanged with varying setup times. 16 3.6 Comparing Distributions 3 SIMULATION RESULTS (a) ν=100s (b) ν=150s (c) ν=200s Figure 9: Average Energy consumption with varying setup times Average Wait Time: As with the average energy consumption trend, Figure 10 shows that the trend across varying distributions remains unchanged with respect to setup time. Furthermore, the limit of each distribution approaches zero as µ → ∞ for all setup times at about the same rate. 17 3.6 Comparing Distributions 3 SIMULATION RESULTS (a) ν=100s (b) ν=150s (c) ν=200s Figure 10: Average wait time with varying setup times 3.6.2 Comparing Number of Servers We can also run the simulation with varying number of servers to see whether the observed patterns hold. The results in Figure 11 are from a system with 50 servers with ν = 100s and λ = 0.4. We note that the observed patterns remain constant within these graphs as well, further contributing to the hypothesis that TABS is consistent independent of the service time distribution. 18 3.6 Comparing Distributions 3 SIMULATION RESULTS (a) Average Energy Consumption (b) Average wait time Figure 11: Comparing Multiple Service Time Distributions 3.6.3 Comparing Variance Effect Finally, we can take a closer look at the effect of changing the variance of the service time distribution. For simplicity of adjusting variance, we consider only the Folded Normal distri- bution. Each point on the graphs in Figures 12 and 13 is obtained by running the simulation 10 times and taking the average value to remove the randomness of any one simulation. We generate each average wait time and average energy consumption figure for various different arrival rates to observe the interaction between changing variance and arrival rates. The average setup time in all simulations is fixed at 150s. Average Energy Consumption: Figure 12 shows that as the variance of the distribution increases, the average energy consumed by the system increases across all µ values. (a) λ=0.15 (b) λ=0.2 19 3.6 Comparing Distributions 3 SIMULATION RESULTS (c) λ=0.4 Figure 12: Average energy consumption with varying arrival rates The general trend in Figure 12 is that energy consumption falls as µ increases holds across all combinations of distribution variance and arrival rates. Furthermore, the magni- tude of average energy consumption does not change drastically across varying arrival rates. Average Wait time: In Figure 13 we see once again that as variance increases, the average task wait time increases. Upon close inspection of the graphs, we note the main difference in these trends from the energy consumption trends in Figure 12 is that the increase in average wait time seems to persist only for lower µ values. After about µ = 70 for arrival rates λ = 0.15, 0.2 the difference in wait time across service distributions with different variances appear to vanish. For the graph with λ = 0.4 this difference seems negligible for µ ≥ 100 for all service distributions except the distribution with variance of 2. The distribution with variance of 2 in the system with λ = 0.4 seems to realize consistently higher wait times for all µ values. However, we maintain that the general trend across different variance values remains the same. (a) λ=0.15 (b) λ=0.2 20 4 LIMIT THEOREM PROOF (c) λ=0.4 Figure 13: Average wait time with varying arrival rates It is worth noting that this analysis looks only at first order behaviour of the system and that by conducting second order analysis on the average wait time and number of servers in setup mode, we may obtain a better understanding of the impact of different service time distributions on the performance of the TABS Scheme. We conclude this section with the hypothesis that TABS behaviour is independent of the service time distribution and, to an extent, its variance conditioned upon the expectation of the distribution being 1. Each of the figures in this section provide an argument in favor of this hypothesis. Further sensitivity analysis will allow for a rigorous proof of this hypothesis. 4 Limit Theorem Proof We now move to discuss the behaviour of the queuing system under TABS with respect to µ and ν. As noted in Section 2, previous results for the TABS scheme have proven its convergence limit and stationary distribution under fluid-limit scaling [11, 12]. The limits proved to be independent of both µ and ν parameters [11]. From the discussion in Section 3.1, we are motivated to find the dependence of mean behaviour of the system on µ and ν values for long time intervals. This result will allow large-scale data centers to optimize µ and ν with respect to system performance. In this section we present a proof for the convergence of the system representation under TABS with diffusion scaling to a multi-dimensional Brownian motion that captures the effect of µ and ν on the system. Eschenfeldt and Gamarnik [3] provide a proof for the diffusion scaling limit theorem for a queuing system with all servers always on. They consider a slightly modified version of the scheme. Specifically, they prevent any server’s queue from being greater √ than 2. They also restrict the number of servers with queue length 2 to be no more that B n. In the limit these restrictions prove to emulate the system without any restrictions on queue lengths. These 21 4.1 System Evolution 4 LIMIT THEOREM PROOF assumptions are helpful in simplifying the analysis and largely do not affect the validity and application of the result to the general model. Our goal is to find a similar limiting behaviour for the TABS scheme. We √ will also assume that queue lengths cannot be greater than 2 and that a maximum of B n servers can have a queue length of 2. Additionally, we will assume that the service time distribution is exponential with expectation 1 as was done in the initial proposal of the TABS scheme [12, 11]. In the following analysis Qni (t) represents the number of queues with at least i customers at time t ≥ 0. We know that Qn1 (t) ≤ n, ∀t ≥ 0 We start by introducing a centered and scaled version of the system, defined by (Qn1 (t) − λn) X1n (t) = √ (4.1) n Qn2 (t) X2n (t) = √ (4.2) n ∆n0 (t) − (1 − λ)n ∆n0 (t) = √ (4.3) n ∆n1 (t) ∆n1 (t) = √ (4.4) n In addition to the earlier bounds set on Qn1 (t) and Qn2 (t), we will require the following equation to hold for all t ≥ 0. This is the reflected diffusion constraint of the system. ∆n0 (t) + ∆n1 (t) + Qn1 (t) ≤ n (4.5) 4.1 System Evolution Given the scaled and centered representations of the 4 random variables that define the system, we can build equations to show how each variable evolves with time. Equations 4.6-4.9 show the evolution of the system with the random variables that model the system. Z t  n n Q1 (t) = Q1 (0) + A1 (λnt) − D1 (Q1 (s) − Q2 (s))ds − U1n (t) n n (4.6) 0 Z t  Qn2 (t) = Qn2 (0) + U1n (t) − D2 Qn2 (s)ds − U2n (t) (4.7) 0 Z t  ∆n0 (t) = ∆n0 (0) + A2 µ(n − Qn1 (s) − ∆n1 (s))ds 0 Z t  −D3 1 {Qn1 (s) + ∆n0 (s) + ∆n1 (s) = n)}1 {∆n0 (s) > 0}dA(λns) (4.8) 0 22 4.2 Martingale Equations 4 LIMIT THEOREM PROOF Z t ∆n1 (t) = ∆n1 (0) + 1{Qn1 (s) + ∆n0 (s) + ∆n1 (s) = n)}1{∆n0 (s) > 0}dA(λns) 0 Z t  −D4 ν∆n1 (s)ds (4.9) 0 In the equations above, Ai (·) and Dk (·) represent point processes that account for increases and decreases, respectively, for each variable over time. U1n (t) is the number of arrivals in [0,t] when every server has at least one task. U2n√(t) is the number of arrivals in [0,t] when every server has at √ least one task and exactly B n servers have two tasks. (Implicitly we require that B ≤ n). Below we represent how U1n (t) and U2n (t) evolve with time. We also rewrite the evolution in terms of the scaled and centered random variables from 4.1-4.4. This will prove useful for calculations in the following section. Z t n U1 (t) = 1{Qn1 (s) − n + ∆n0 + ∆n1 = 0}dA(λns) 0 Z t = 1{Qn1 (s) − λn + ∆n0 − (1 − λ)n + ∆n1 = 0}dA(λns) 0 t Qn1 (s) − λn + ∆n0 − (1 − λ)n + ∆n1 Z   = 1 √ = 0 dA(λns) 0 n Z t 1 X1n (s) + ∆n0 (s) + ∆n1 (s) = 0 dA(λns)  = (4.10) 0 Z t 1{Qn1 (s) − n + ∆n0 + ∆n1 = 0, Qn2 (s) − B n − ∆n0 − ∆n1 = 0}dA(λns) p U2n (t) = 0 ( ) p t n − ∆n0 − ∆n1 Qn2 (s) − B Z = 1 Qn1 (s) − λn + − (1 − λ)n + = 0, ∆n0 √ ∆n1 = 0 dA(λns) 0 n Z t ( n p ) Q1 (s) − λn + ∆n0 − (1 − λ)n + ∆n1 Qn2 (s) − B n − ∆n0 − ∆n1 = 1 √ = 0, √ = 0 dA(λns) 0 n n Z t ( r n n ) λn + (1 − λ)n − ∆ − ∆ = 1 X1 (s) + ∆n0 (s) + ∆n1 (s) = 0, X2 (s) − B n n 0 1 = 0 dA(λns) 0 n Z t  √ q  = 1 X1 + ∆0 + ∆1 = 0, X2 − B λ − (∆0 − ∆1 ) n = 0 dA(λns) n n n n n n (4.11) 0 4.2 Martingale Equations Now we can rewrite the time changes of Poisson processes from 4.6-4.9 as time changes of scaled Poisson processes. To do this, we first define scaled martingales that will simplify the 23 4.3 Martingale Decomposition 4 LIMIT THEOREM PROOF scaled and centered representation of the system’s evolution in the next section. 1 √ M1n,1 = √ A1 (λnt) − λ nt (4.12) n Z t  Z t 1 1 M1n,2 = √ A2 µ(n − Qn1 (s) − ∆n1 (s))ds −√ µ(n − Qn1 (s) − ∆n1 (s))ds (4.13) n 0 n 0 Z t Z t  1 1 M2n,1 = √ D1 − (Qn1 (s) −√ Qn2 (s))ds (Qn1 (s) − Qn2 (s))ds (4.14) n 0 n 0 Z t  Z t n,2 1 n 1 M2 = √ D2 Q2 (s)ds − √ Qn2 (s)ds (4.15) n 0 n 0 Z t  1 M2n,3 = √ D3 1{Q1 (s) + ∆0 (s) + ∆1 (s) = n)}1{∆0 (s) > 0}dA(λns) n n n n n 0 Z t 1 −√ 1{Qn1 (s) + ∆n0 (s) + ∆n1 (s) = n)}1{∆n0 (s) > 0}dA(λns) (4.16) n 0 Z t  Z t 1 1 M2n,4 = √ D4 ν∆n1 (s)ds −√ ν∆n1 (s)ds (4.17) n 0 n 0 For simplicity of calculations in the following sections, we also introduce scaled versions of U1n (t) and U2n (t) n U1n (t) n U2n (t) V1 (t) = √ and V2 (t) = √ (4.18) n n 4.3 Martingale Decomposition Given the martingales from section 4.2, we can center, scale and simplify each of the evo- lution equations (4.6-4.9) to get an appropriate martingale decomposition to represent the entire system. Starting with Qn1 (t), we subtract and add the expectation of each of the Poisson processes Qn1 (t) = Qn1 (0) + (A1 (λnt) − λnt) + λnt  Z t  Z t  n n n n − D1 (Q1 (s) − Q2 (s))ds − (Q1 (s) − Q2 (s))ds 0 0 Z t − (Qn1 (s) − Qn2 (s))ds − U1n (t) 0 24 4.3 Martingale Decomposition 4 LIMIT THEOREM PROOF 1 Next, we scale by √ and rearrange the variables. Note that we now have martingales n constructed in 4.12-4.17 appear on the RHS above. By substituting for these martingales, we get the following equation: t (Qn1 (t) − λn) Qn (0) − λn Z  1 √ = 1 √ + M1n,1 (t) − M2n,1 (t) − √ (Qn1 (s) − Qn2 (s))ds + U1n (t) n n n 0 Finally, we substitute in the scaled representations to get a final martingale decomposition for the first random variable of the system. Z t n n n,1 n,1 X1 (t) = X1 (0) + M1 (t) − M2 (t) − (X1n (s) − X2n (s))ds − V1n (t) 0 We repeat this process for each of the other three evolution equations from Section 4.1. For each of the remaining three equations from 4.6-4.9, first we center, then we scale and finally we substitute the martingale representations. Z t  n n n Q2 (t) = Q2 (0) + U1 (t) − D2 Q2 (s)ds − U2n (t) n  0Z t  Z t  Z t n n n n = Q2 (0) + U1 (t) − D2 Q2 (s)ds − Q2 (s)ds − Qn2 (s)ds − U2n (t) 0 0 0  R  R  R t n t n t Qn2 (t) Qn2 (0) + U1n (t) − D2 0 Q2 (s)ds − 0 Q2 (s)ds − 0 Qn2 (s)ds − U2n (t) √ = √ n n Z t X2n (t) = X2n (0) + V1n (t) − M2n,2 (t) − X2n (s)ds − V2n (t) 0 Z t  ∆n0 (t) = ∆n0 (0) + A2 µ(n − Qn1 (s) − ∆n1 (s))ds 0 Z t  −D3 1 {Qn1 (s) + ∆n0 (s) + ∆n1 (s) = n)}1 {∆n0 (s) > 0}dA(λns) 1 0 = (∆n0 (0) − (1 − λ)n) + (1 − λ)n  Z t  Z t  n n n n + A2 µ(n − Q1 (s) − ∆1 (s))ds − µ(n − Q1 (s) − ∆1 (s))ds 0 0  Z t  Z t  − D3 1{∆0 (s) > 0}dU1 (s) − 1{∆0 (s) > 0}dU1 (s) n n n n 0 Z t Z t 0 + µ(n − Qn1 (s) − ∆n1 (s))ds − 1{∆n0 (s) > 0}dU1n (s) 0 0 1 Note that by nature of the indicator function inside the integral in M2n,3 , we can simplify the integral to Z t Z t 1{Qn1 (s) + ∆n0 (s) + ∆n1 (s) = n)}1{∆n0 (s) > 0}dA(λns) = 1{∆n0 (s) > 0}dU1n (s) 0 0 25 4.3 Martingale Decomposition 4 LIMIT THEOREM PROOF ∆n0 (t) (∆n0 (0) − (1 − λ)n) + (1 − λ)n √ = √ n n  R  R  t t + A2 0 µ(n − Qn1 (s) − ∆n1 (s))ds − 0 µ(n − Qn1 (s) − ∆n1 (s))ds √ n  R  R  − D3 0 1{∆n0 (s)n > 0}dU1n (s) − 0 1{∆n0 (s) > 0}dU1n (s) t t √ n + 0 µ(n − Qn1 (s) − ∆n1 (s))ds − 0 1{∆n0 (s) > 0}dU1n (s) Rt Rt √ n t µ(n − Qn1 (s) − ∆n1 (s)) Z ∆n0 (t) = ∆n0 (0) + M1n,2 (t) − M2n,3 (t) + √ ds 0 n t 1{∆n0 (s) > 0} Z − √ dU1n (s) 0 n t µ(λn + (1 − λ)n − Qn1 (s) − ∆n1 (s)) Z = ∆n0 (0) + M1n,2 (t) − M2n,3 (t) √ + ds 0 n 1{∆n0 (s) − (1 − λ)n > −(1 − λ)n} n Z t − √ dU1 (s) 0 n Z t n,2 n,3 √ n = ∆0 (0) + M1 (t) − M2 (t) + µ( n(1 − λ) − X1n (s) − ∆n1 (s))ds 0 Z t √ − 1{−∆n0 (s) < n(1 − λ)}dU1n (s) 0 Z t √ ∆n0 (t) = ∆n0 (0) + M1n,2 (t) − M2n,3 (t) + µ( n(1 − λ) − X1n (s) − ∆n1 (s))ds 0 Z t √ − 1{−∆n0 (s) < n(1 − λ)}dU1n (s) 0 Z t ∆n1 (t) = ∆n1 (0) + 1{Qn1 (s) + ∆n0 (s) + ∆n1 (s) = n)}1{∆n0 (s) > 0}dA(λns) 0 Z t  −D4 ν∆n1 (s)ds 0 Z t = ∆n1 (0) + 1{∆n0 (s) > 0}dU1n (s)  Z0 t  Z t  Z t n n − D4 ν∆1 (s)ds − ν∆1 (s)ds − ν∆n1 (s)ds 0 0 0 26 4.4 Martingale Behaviour 4 LIMIT THEOREM PROOF ∆n1 (0) + 0 1{∆n0 (s) > 0}dU1n (s) Rt ∆n1 (t) √ = √ n n  R  R  R t t t − D4 0 ν∆n1 (s)ds − 0 ν∆n1 (s)ds − 0 ν∆n1 (s)ds √ n t 1{∆n0 (s) − (1 − λ)n > −(1 − λ)n} t ν∆n1 (s) Z Z ∆n1 (t) = ∆n1 (0) − M2n,4 (t) + √ dU1n (s) − √ ds 0 n 0 n Z t √ Z t = ∆n1 (0) − M2n,4 (t) + 1{−∆n0 (s) < n(1 − λ)}dU1n (s) − ν∆n1 (s)ds 0 0 Z t √ Z t ∆n1 (t) = ∆n1 (0) − M2n,4 (t) + 1 {−∆n0 (s) < n(1 − λ)}dU1n (s) − ν∆n1 (s)ds 0 0 Thus, we have the final four martingale decompositions: Z t n n n,1 n,1 X1 (t) = X1 (0) + M1 (t) − M2 (t) − (X1n (s) − X2n (s))ds − V1n (t) (4.19) 0 Z t X2n (t) = X2n (0) + V1n (t) − M2n,2 (t) − X2n (s)ds − V2n (t) (4.20) 0 Z t √ ∆n0 (t) = ∆n0 (0) + M1n,2 (t) − M2n,3 (t) + µ( n(1 − λ) − X1n (s) − ∆n1 (s))ds 0 Z t √ − 1{−∆n0 (s) < n(1 − λ)}dU1n (s) (4.21) 0 Z t √ Z t ∆n1 (t) = ∆n1 (0) − M2n,4 (t) + 1 {−∆n0 (s) < n(1 − λ)}dU1n (s) − ν∆n1 (s)ds (4.22) 0 0 4.4 Martingale Behaviour 4.4.1 Quadratic Variation We now turn to inspect each of the martingale equations in further detail. The first step is to identify the square-integrable martingales and their quadratic variations and covariations with each martingale. 27 4.4 Martingale Behaviour 4 LIMIT THEOREM PROOF Definition 4.1. A martingale, M (t) ≡ {M (t) : t ≥ 0}(with respect to some filtration), is square-integrable if E[M (t)2 ] ≤ ∞ ∀ t ≥ 0. [13] From Theorem 7.2 of [13], we can identify M1n,1 , M1n,2 , M2n,1 , M2n,2 , M2n,3 , M2n,4 as square integrable with respect to an appropriate filtration. We also obtain the following predictable quadratic variations(PQV) for each of the martingales. hM1n,1 i(t) = λt (4.23) Z t 1 hM1n,2 i(t) = µ(n − Qn1 (s) − ∆n1 (s))ds (4.24) n 0 Z t 1 hM2n,1 i(t) = (Qn1 (s) − Qn2 (s))ds (4.25) n 0 Z t 1 hM2n,2 i(t) = Qn2 (s)ds (4.26) n 0 Z t 1 hM2n,3 i(t) = 1{∆n0 (s) > 0}dU1n (s) (4.27) n 0 Z t 1 hM2n,4 i(t) = ν∆n1 (s)ds (4.28) n 0 4.4.2 Quadratic Covariation Now, we find the quadratic covariation for each pair of martingales. From [13] we have the formula to derive the quadratic covariation of two square-integrable martingales as 1 hM1 , M2 i = (hM1 + M2 i − hM1 i − hM2 i) 2 To determine the quadratic covariation, we first need to identify the instantaneous rate change of each of the point processes that make up the each of the martingales. Once we have the functions, we can compute the variance of each pair to identify the limiting behaviour and thus the quadratic covariation between each pair of martingales. The instantaneous rate change of each point process below is identified by the taking the derivative of the rate function with respect to the time variable t. A1 (λnt) ⇒ λn Z t  A2 µ(n − Q1 (s) − ∆1 (s))ds ⇒ µ(n − Qn1 (t) − ∆n1 (t)) n n 0 Z t  D1 (Qn1 (s) − Qn2 (s))ds ⇒ Qn1 (t) − Qn2 (t) 0 28 4.4 Martingale Behaviour 4 LIMIT THEOREM PROOF Z t  D2 Qn2 (s)ds ⇒ Qn2 (t) 0 Z t  D3 1{∆0 (s) > 0}dU1 (s) ⇒ 1{∆n0 (t) > 0} n n 0 Z t  D4 ν∆1 (s)ds ⇒ ν∆n1 (t) n 0 Given these functions, we can now determine the quadratic covariance between each pair of martingales. Note that for the M1n,1 the instantaneous rate change is a constant function, i.e. the function does not depend on t. This gives us the feature that the covariance of M1n,1 with respect to all other martingales will be 0. Additionally, for M2n,3 the instantaneous rate is an indicator function which is assumed to be 1 for the purpose of this system set up. This is because from previous results we know that ∆n0 (t) hovers around (1 − λ)n and so the probability that ∆n0 (t) reaches 0 is exponentially decreases as n increases. Thus the quadratic covariation of M2n,3 with all other martingales is 0. We calculate the quadratic covariation of the rest of the paired equations below: Cov(µ(n − Qn1 (t) − ∆n1 (t)), Qn1 (t) − Qn2 (t)) = −µVar(Qn1 (t)) + µCov(Qn1 (t), Qn2 (t)) −µCov(∆n1 (t), Qn1 (t)) + µCov(∆n1 (t), Qn2 (t)) √ √ √ √ = µ[−o( n) + o( n) + O( n − O( n)] √ If we scale this covariance by n, we can obtain the limit for the quadratic covariance between M1n,2 and M2n,1 : 1 √ Cov(µ(n − Qn1 (t) − ∆n1 (t)), Qn1 (t) − Qn2 (t)) ⇒ 0 n hM1n,2 , M2n,1 i = 0 (4.29) 1 1 √ Cov(µ(n − Qn1 (t) − ∆n1 (t)), Qn2 (t)) = √ µCov(Qn1 (t), Qn2 (t)) − µCov(∆n1 (t), Qn2 (t)) n n 1 √ √ = √ µ(o( n) − O( n)) n ⇒ −aµ for some constant a > 0 Thus, the quadratic covariation is given by: 1 t 1 t n  Z Z  n,2 n,2 1 n n hM1 , M2 i = −aµ − µ(n − Q1 (s) − ∆1 (s))ds − Q (s)ds (4.30) 2 n 0 n 0 2 Repeating this process for the rest of the pairs of martingales we get the following equations: 1 1 √ Cov(µ(n − Qn1 (t) − ∆n1 (t)), ν∆n1 (t)) = √ µνCov(Qn1 (t), ∆n1 (t)) − µνVar(∆n1 (t)) n n 1 √ p = √ µν(O( n) − o( n)) n ⇒ bµν for some constant b > 0 29 4.4 Martingale Behaviour 4 LIMIT THEOREM PROOF  Z t Z t  1 1 1 hM1n,2 , M2n,4 i = bµν − µ(n − Qn1 (s) − ∆n1 (s))ds − ν∆n1 (s)ds (4.31) 2 n 0 n 0 1 1 √ Cov(Qn1 (t) − Qn2 (t), Qn2 (t)) = √ Cov(Qn1 (t), Qn2 (t)) − Var(Qn2 (t)) n n 1 √ √ = √ o( n) − o( n) n ⇒0 hM2n,1 , M2n,2 i = 0 (4.32) 1 1 √ Cov(Qn1 (t) − Qn2 (t), ν∆n1 (t)) = √ νCov(Qn1 (t), ∆n1 (t)) − νCov(Qn2 (t), ∆n1 (t)) n n 1 √ √ = √ ν(O( n) − O( n)) n ⇒0 (4.33) hM2n,1 , M2n,4 i = 0 (4.34) 1 1 √ Cov(Qn2 (t), ν∆n1 (t)) = √ νCov(Qn2 (t), ∆n1 (t)) n n 1 √ = √ νO( n) n ⇒ cν for some constant c > 0 (4.35)  Z t Z t  1 1 1 hM2n,2 , M2n,4 i = cν − Qn2 (s)ds − ν∆n1 (s)ds (4.36) 2 n 0 n 0 Furthermore, based on Equations 4.10 and 4.6-4.9, we can derive some additional system constraints. By Equation 4.10, we have: Z ∞ 0= 1{Qn1 (s) ≤ n − ∆n0 − ∆n1 }dU1n (s) 0 Z ∞ = 1{Qn1 (s) + ∆n0 + ∆n1 − λn − (1 − λ)n < 0}dU1n (s) 0 Z ∞ = 1{X1n (s) + ∆n0 (s) + ∆n1 (s) < 0}dV1n (s) (4.37) 0 30 4.5 Martingale Convergence 4 LIMIT THEOREM PROOF √ Additionally, √ since U2n (s) can increase only when Qn1 (s) = n−∆n0 −∆n1 and Qn2 (s) = B n, if Qn2 (s) < B n, U2n (s) cannot increase. From this, we arrive at the following condition: Z ∞ √ 0= 1{Qn2 (s) < B n}dU2n (t) 0 Z ∞ = 1{X2n (s) < B}dV2n (t) (4.38) 0 The next step in our proof requires convergence of the martingales from equations 4.19-4.22. 4.5 Martingale Convergence Given the martingales are square integrable and we have the predictable quadratic variations of each martingale, we can identify the limits of each PQV. We first define the time changes as follows. ΦA1 ,n (t) = λt (4.39) µ t n 1 t Z Z ΦA2 ,n (t) = µt − Q (s)ds − µ∆n1 (s)ds (4.40) n 0 1 n 0 1 t n 1 t n Z Z ΦD1 ,n (t) = Q (s)ds − Q (s)ds (4.41) n 0 1 n 0 2 1 t n Z ΦD2 ,n (t) = Q (s)ds (4.42) n 0 2 1 t Z ΦD3 ,n (t) = 1{∆n0 (s) > 0}dU1n (s) (4.43) n 0 ν t n Z ΦD4 ,n (t) = ∆ (s)ds (4.44) n 0 1 Now we can rewrite the martingales as M1n,1 = MA1 ,n ◦ ΦA1 ,n , M1n,2 = MA2 ,n ◦ ΦA2 ,n , M2n,1 = MD1 ,n ◦ ΦD1 ,n M2n,2 = MD2 ,n ◦ ΦD2 ,n , M2n,3 = MD3 ,n ◦ ΦD3 ,n , M2n,4 = MD4 ,n ◦ ΦD4 ,n Next, we define the limits of the time change equations (4.39-4.44). Before identifying the limits, we define the topological space over which these processes are defined. Let D = D([0, ∞), R) be the space of right continuous functions with left limits mapping [0, ∞) in R. Where necessary, Dk = D([0, ∞), Rk ) for k ≥ 0 is the product space of these functions. Thus all the following functions and limits are with respect to the domain D. Starting with ΦA1 ,n , since λ ∈ (0, 1) is a constant, ΦA1 ,n ⇒ λt as n → ∞ (4.45) 31 4.5 Martingale Convergence 4 LIMIT THEOREM PROOF Next we look at ΦD1 ,n . Note that in ΦD1 ,n the second term is exactly ΦD2 ,n so ΦD1 ,n ⇒ f − g as n → ∞ where g is the limit of ΦD2 ,n and f is the limit of the first term of ΦD1 ,n . To find f and g, we need to show the fluid limits for Qn1 (s) and Qn2 (s). We find the fluid limit just as was done in Lemma 6 in [3], reproduced below. Lemma 4.2. Let Ψni for i = 1, 2 be defined as Qni (t) Ψni (t) = , t ≥ 0. n Then, Ψn1 ⇒ ω and Ψn2 ⇒ 0 as n → ∞ (4.46) where ω(t) = λ for t ≥ 0 since λ is a constant in our system. The proof of this lemma is given in Section 5.1 of [3]. From this we have that f is the identity function in D. That is, 1 t n Z Q (s)ds ⇒ λt as n → ∞. n 0 1 Rt Similarly, g(t) = 0 0ds = 0, so g = 0 on D. Putting the two limits together, we have ΦD1 ,n ⇒ λt as n → ∞ (4.47) and ΦD2 ,n ⇒ 0 as n → ∞ (4.48) ν Rt n Next we look at ΦD4 ,n = ∆ (s)ds. Since ν represents the average startup time for a n 0 1 server, we know ν is generally some fixed constant, c, independent of other parameters in R t ∆n (s) the system. Thus, we only need to understand the convergence of the integral 0 1 ds. n We find the fluid limit for ∆n1 (s) just as was done for Qn1 (s) and Qn1 (s) Lemma 4.3. Let Ψn3 be defined as ∆n1 (t) Ψn3 (t) = , t ≥ 0. n Then, Ψn3 ⇒ 0 as n → ∞ (4.49) The proof of this lemma is identical to that of Lemma 4.2. Thus, we have that given ν → c, ΦD4 ,n ⇒ 0 (4.50) 32 4.6 Martingale Limit 4 LIMIT THEOREM PROOF For ΦA2 ,n , note that the second term is µ times the first term of ΦD1 ,n and the second term is precisely ΦD4 ,n . Thus, we can write the limit as ΦA2 ,n ⇒ k − µf − h where k is the limit of ΦA2 ,n with ΦA2 ,n = µt and h is the limit of ΦD4 ,n . To find the limit of k we can use the fact that µ is a constant. ΦA2 ,n ⇒ µt as n → ∞ From previous calculations, we have that f ⇒ λt as n → ∞ and ΦD4 ,n ⇒ 0. Putting all the limits together, we have ΦA2 ,n ⇒ µ(1 − λ)t as n → ∞ (4.51) Finally, we look at ΦD3 ,n Z t 1 ΦD3 ,n (t) = 1{∆n0 (s) > 0}dU1n (s) (4.52) n 0 As before, we assume that the indicator function 1{∆n0 (s) > 0} is equal to 1 for all time intervals. This assumption is made from earlier results with respect to the TABS scheme limiting behaviour ∆n0 approaches (1−λ). The probability that ∆n0 will hit 0 is exponentially unlikely. Using this assumption we get the following convergence result: 1 ΦD3 ,n (t) = ⇒ 0 as n → ∞ (4.53) n 4.6 Martingale Limit We can now put together the martingales and show their convergence. That is, (M1n,1 , M1n,2 , M2n,1 , M2n,2 , M2n,3 , M2n,4 ) =⇒ (W1 , W2 , W3 , W4 , W5 , W6 ) (4.54) with a zero mean vector and Σ covariance matrix.   λt 0 0 0 0 0  0 µt 0 f (t) 0 g(t)   0 0 λt 0 0 0  Σ=  0 f (t) 0  (4.55)  0 0 h(t)  0 0 0 0 0 0  0 g(t) 0 h(t) 0 0 where 1 f (t) = (−aµ − 1 − λt) 2 33 5 CONCLUSION 1 g(t) = (bµν − µ − µλt) 2 cν h(t) = 2 for fixed a, b, c > 0 Thus we have identified the diffusion limit for the TABS system with exponential ser- vice time distribution under diffusion scaling. The system converges to Brownian motion, as expected with a specified covariance matrix. This completes the proof of martingale convergence. We now have a result that fits into the criteria of applying the Skorokhod map [9] to solve the stochastic differential equations governing the system and thus identify the dependence of mean behaviour of the system variables on µ and ν. The final result can then be understood as a Taylor expansion of mean behaviour of each of the random variables of interest. For example, the mean behaviour for Qn1 (t) can be broken down into: √ √ ˆ n (t) + n · X n (t) + o( n) Qn1 (t) = n · Q (4.56) 1 1 Here we see that the first part of the approximation for mean behaviour of Qn1 (t) was the result obtained by Mukherjee et. al. [12]. The completion of this proof provides the second term in Equation 4.56 (in red). 5 Conclusion 5.1 Summary of Results The TABS scheme proposed to optimize data center performance through load-balancing is known to be stable in the limit with respect to time and number of servers under fluid scaling. Previous results about the TABS scheme have assumed an exponential service time distribution. In this thesis, we investigated the behaviour of the TABS scheme under alternative service time distributions through simulations and presented a proof for TABS limit convergence under diffusion scaling conditional on certain system properties being satisfied. The first section of this thesis presented simulations of the TABS scheme on large-scale data centers with varying service time distributions. We considered three different service time distributions- generalized Pareto, Folded Normal and Hyperexponential service time distributions with expectation λ and variance 1. A comparison of performance metrics of average task wait time and average system energy consumption showed that varying service time distributions did not contribute to larges shifts in performance of the TABS scheme in comparison with the standard exponential distribution. We also simulated the TABS scheme under the Folded Normal service time distribution for a range of variance values. Within this range, the TABS scheme is consistent with respect to the service time distribution variance. We observed that higher variance typically led to higher average wait times for tasks in the system as as well as higher overall energy consumption of the data center. The 34 5.2 Open Questions 5 CONCLUSION simulated data center performance allows us to hypothesize the robustness of the TABS scheme’s performance under varying service time distribution parameters. The second half of this thesis presented a convergence proof for the diffusion-scaled TABS system under specific assumptions. We modified the typical setup of the TABS scheme from literature to scale each variable of interest around its long-term expected value. This diffusion scaling of the system variables allowed for a closer look at the randomness of the variables in the limit. We provided a Martingale representation of the four fundamental equations that characterize the evolution of the system under the TABS scheme and use information about the predictable quadratic variance and covariance of the martingales in the equations to identify limit for the system of martingales. The proof showed convergence of the set of Martingales to a multidimensional Brownian motion for which we explicitly calculated the mean vector and covariance matrix. 5.2 Open Questions The simulation results in this thesis appear to corroborate the hypothesis that the TABS scheme contributes to improved performance of large-scale data centers regardless of the service time distribution and to an extent, the variance of the distribution itself. One open question from the results of this thesis is whether the minor differences observed in simulation results with varying service time distributions is statistically significant. Future work could conduct a sensitivity analysis on the TABS scheme to investigate the robustness of the scheme to varying service time distributions. Furthermore, the results presented in this thesis looked at first order effects of √ changing service time distributions their variance. Second order analysis, done by looking at n scaled behaviour of number of servers in setup mode or average wait time in the system will provide a thorough analysis of how service time distributions and their variances affect the performance of the queuing system. Further work on the simulation results to take into account the varying tails of each of the distributions to observe their effect of the performance of the system. This could include analysis of TABS with additional service time distributions not considered in this thesis. The proof of the system convergence can be completed by the application of the Sko- rokhod map [9]. Additionally, the proof of the TABS system convergence under diffusion scaling in this thesis uses some assumptions on typical server behaviour, we also imposed constraints on number of servers with queue lengths greater than or equal to 1. The proof also assumes that the likelihood of extreme events that may impact the overall performance of the system to be exponentially small. To completely justify the conclusions of this thesis, one would need to prove that these assumptions are indeed satisfied by the system. 35 REFERENCES REFERENCES References [1] Aghajani, Reza, Xingjie Li, and Kavita Ramanan. ”The PDE method for the analysis of randomized load balancing networks.” Proceedings of the ACM on Measurement and Analysis of Computing Systems 1, no. 2 (2017): 38. [2] Barroso, Luiz Andr, and Urs Hlzle. ”The case for energy-proportional computing.” Com- puter 12 (2007): 33-37. [3] Eschenfeldt, Patrick, and David Gamarnik. ”Join the Shortest Queue with Many Servers. The Heavy-Traffic Asymptotics.” Mathematics of Operations Research (2018). [4] Gandhi, Anshul, et al. ”Exact analysis of the M/M/k/setup class of Markov chains via recursive renewal reward.” ACM SIGMETRICS Performance Evaluation Review. Vol. 41. No. 1. ACM, 2013. [5] Gandhi, Anshul, Mor Harchol-Balter, and Ivo Adan. ”Server farms with setup costs.” Performance Evaluation 67, no. 11 (2010): 1123-1138. [6] Greenberg, Albert, et al. ”The cost of a cloud: research problems in data center net- works.” ACM SIGCOMM computer communication review 39.1 (2008): 68-73. [7] Koomey, Jonathan G. ”Estimating total power consumption by servers in the US and the world.” (2007). [8] Koomey, Jonathan G. ”Worldwide electricity used in data centers.” Environmental re- search letters 3, no. 3 (2008): 034008. [9] Kruk, Lukasz, John Lehoczky, Kavita Ramanan, and Steven Shreve. ”An explicit formula for the Skorokhod map on [0, a].” The Annals of Probability 35, no. 5 (2007): 1740-1768. [10] Liu, Zhenhua, et al. ”Renewable and cooling aware workload management for sustain- able data centers.” ACM SIGMETRICS Performance Evaluation Review. Vol. 40. No. 1. ACM, 2012. [11] Mukherjee, Debankur, and Alexander Stolyar. ”Join-idle-queue with service elasticity: large-scale asymptotics of a non-monotone system.” arXiv preprint arXiv:1803.07689 (2018). [12] Mukherjee, Debankur, et al. ”Optimal service elasticity in large-scale distributed sys- tems.” Proceedings of the ACM on Measurement and Analysis of Computing Systems 1.1 (2017): 25. [13] Pang, Guodong, Rishi Talreja, and Ward Whitt. ”Martingale proofs of many-server heavy-traffic limits for Markovian queues.” Probability Surveys 4 (2007): 193-267. [14] Stankovic, John A. ”Simulations of three adaptive, decentralized controlled, job schedul- ing algorithms.” Computer Networks (1976) 8, no. 3 (1984): 199-217. 36