Eco-Friendly Powering and Delay-Aware Task Scheduling in Geo-Distributed Edge-Cloud System: A Two-Timescale Framework

Power consumption and task latency are two crucial issues in edge-cloud computing. This paper mainly aims to promote the use of clean power in geo-distributed data centers (DCs) in a deregulated electricity market where customers are allowed to buy power from multiple suppliers, combined with the guarantee of task latency. To alleviate the conflict between frequent switches of servers and the uncertainty of task arrivals in DCs, this paper proposes a two-timescale framework consisting of the long-term capacity planning of geo-distributed DCs and the real-time task dispatching from edge gateways (EGs) to DCs. First, DCs make long-term plans on the number of active servers aiming at the eco-friendly and delay-aware power cost minimization, which is formulated as problem $\mathcal {P}$ . Specifically, we introduce a convex pollution indicator function (PIF) to measure the pollution cost of the various types of powers sold by different suppliers, which can encourage the use of cleaner power and improve power savings. Second, in each sub-slot, each EG separately optimizes its individual mixed strategies of task dispatching to DCs with the knowledge of the planned capacities and the real-time queue backlogs of DCs, where a Lyapunov optimization framework is applied. Finally, we give the corresponding distributed algorithm design. Simulation results reveal that our method can realize the trade-off between the power cost and the delay cost of requests, and improve the clean power usage by up to 50%–60% of the total power usage in DCs. Additionally, comparisons with other schemes show that our approach can provide more stable guarantees of task latency in different situations of workload density, which benefits from the diverse-timescale optimizations of capacities of DCs and task routing from EGs to DCs.


I. INTRODUCTION
With the penetration of edge-cloud computing, the requests from users are explosively growing and the power consumptions in DCs are becoming enormous. On one hand, given the limitations of system's computation capacity and users' latency tolerance, it is of great significance to improve the efficiency of task scheduling inside the edge-cloud system The associate editor coordinating the review of this manuscript and approving it for publication was Ahmed F. Zobaa . in a situation of increasingly growing requests. On the other hand, the huge power consumptions of expanding DCs can not only increase the power cost of service providers, but also bring about massive carbon emissions along with severe air pollutions of sulfide and dust. According to [1], the total number of servers in the geo-distributed data centers of Google, Microsoft and Akamai were almost 1 million, 200,000 and 70,000 in 2010, respectively, and their corresponding power costs were on the order of millions of dollars per year. As reported in [2], the annual power consumption of DCs across China had reached up to 160 million MW.h, which was almost equivalent to the annual generation of the Three Gorges Dam. It is estimated that just a 1-MW data center powered by thermal power can cause over 10,000 metric tons of CO 2 emissions annually [3].
We consider a deregulated electricity market throughout this paper, which has been gradually popularized across the world with the development of energy Internet in recent years, e.g., in the south of China [4], the state of Texas in the US [5], several Nordic countries [6], etc. In a deregulated electricity market, local electricity companies (ECs) can buy wholesale electricity from different plants and then sell it to customers. At the same time, customers are allowed to freely choose and buy power from one or more ECs according to the price, cleanness and quality of the powers supplied by them. In this situation, geo-distributed DCs can take advantage of diversities of price and cleanness of electricity in different regions to optimize their power uses.
This paper mainly focuses on promoting the use of clean power in geo-distributed data centers (DCs) in a deregulated electricity market, combined with guarantees on task latency. A two-timescale framework is proposed, consisting of long-term capacity planning of geo-distributed DCs and real-time task dispatching from edge gateways (EGs) to DCs. To the best of the authors' knowledge, this is the first time that the number of active servers and the task scheduling strategies are optimized in diverse timescales, which can not only alleviate the conflict between frequent switches of servers and the uncertainty of task arrivals in DCs, but also provide more stable task latency guarantees in different situations of workload density, as shown by our simulation results.

A. RELATED WORK
First, we consider the literature about cost-aware task scheduling in edge-cloud systems. Reference [7] studied the minimization of power-cost via task scheduling among geo-distributed DCs by exploiting the spatial diversity of electricity price, where a boundary constraint of expected sojourn time in DCs was considered. Reference [8] investigated the temporal task scheduling with strict sojourn latency restrictions in a DC aiming to reduce the electricity bill by following the time-varying electricity price, where the considered DC should make decisions on task admissions and executions in each discrete time slot according to a multi-step ahead power cost. Reference [9] considered the spatial task scheduling among geo-distributed DCs on the basis of [8]. However, the transmission delay from sources to DCs was not considered in the literature mentioned above, which generally depended on the propagation distance, link bandwidth, traffic density, etc. Reference [10] modeled the transmission delay as an implicit function of the total task arrivals to a DC. The authors proposed an HBBF algorithm to search for the optimal solution by iteratively feeding back the value of transmission delay associated with the newly updated task scheduling strategies. Additionally, [11] investigated the offloading route selections and the task scheduling among fog nodes, and [12] discussed the joint optimization of link resources in a cloud radio access network (C-RAN), computation resources in mobile edge cloud (MEC) and task scheduling strategies between C-RAN and MEC, both of which applied Lyapunov optimization techniques to design online algorithms of power-delayed balanced online task scheduling. All of the above made great contributions on the cost-aware task scheduling in edge-cloud systems. However, they only considered electricity price in the power cost formulation, while clean power uses, power storage and multi-source power supplies were not mentioned. As stated in [13], the geo-distributed task scheduling in the case of static electricity price or the one decoupled from clean power generations was of nearly no social benefits. In addition, different from the majority of works on geographical task scheduling, [13] focused on the joint optimization of power cost and task latency, instead of the power cost minimization subject to task latency conditions. Then, in the area of eco-friendly power usage, [14] proposed a bid mechanism for the colocation of data centers where equipment, space, and bandwidth were available for rental, and operators could realize carbon-aware task scheduling by using economic incentives to reshape the tenants' demand. Reference [15] investigated the balance of power demands and supplies with the aid of task scheduling among geographically green DCs powered by fuel cells. References [16]- [18] assumed that DCs could harvest power from their private green micro-grids. They proposed online algorithms to minimize the carbon emission and power cost by controlling the number of active servers, and the amount of power bought from power grids and collected from microgrids. Reference [19] studied the joint scheduling of tasks and power storage over a time horizon to reduce the power cost of a DC, considering the temporal diversity of electricity. Reference [20] assumed that DCs could buy both clean power and brown power, and optimized the task scheduling among geo-distributed DCs besides the number of active servers and the purchase amounts of clean power and brown power, aiming to realize the trade-off among the request cost and carbon-inclusive power cost. However, it only simply modeled the cost of requesting servers as a linear function of the number of active servers, while ignoring the task latency experienced by users. Reference [21] studied the match between DCs and power suppliers in the background of the deregulated electricity market, which referred to the optimization of task arrivals to DCs from the perspective of demand response. All of the literature mentioned above only roughly considered the total expected task arrivals to DCs in the considered time slot, while ignoring the detailed process of task routing from different EGs to DCs, which on the other hand would be needed to characterize some crucial factors in practice, such as transmission delay and real time task arrivals to EGs.
To the best of the authors' knowledge, only a small number of papers in the previous literature jointly considered eco-friendly power usage and efficient task scheduling. Reference [22] focused on a carbon-aware power cost minimization problem with sojourn latency constraints of tasks in geo-distributed DCs. Reference [23] pointed out that large financial losses were inevitable once the failure of DCs occurred. To address this issue, it proposed a fault-tolerant task scheduling scheme from users to geo-distributed DCs, accounting for the trade-off between carbon-inclusive power costs and request costs. However, to avoid extreme complexities resulting from the integration of task scheduling and carbon-aware power usage, they ignored some crucial factors, such as the transmission delay, power storage, etc.
Finally, it is worth noting that all the above literature ignored the switch cost of servers in the sense that the operation timescale of switching servers on or off is the same with that of task scheduling. To address this issue, we concentrate our attention on a two-timescale framework design. Papers [24] and [25] were also interested in this. They used an extended Lyapunov optimization framework to investigate the trade-off between power costs and queue backlogs of DCs, which consisted of operations on two different timescales. Specifically, the number of active servers and the strategies for task routing were determined in large time slots, while the CPU frequency was optimized in each small time slot. However, neither transmission delay nor clean power usages were covered in [24] and [25]. In this paper, we determine the capacity provisioning of DCs in large time slots by jointly optimizing the task latency and power cost, as shown in problem P. Then, we apply a Lyapunov optimization framework to optimize EG's task dispatching strategies in each small time slot in a distributed manner. Different from references [24] and [25], we optimize the number of active servers and the task scheduling strategies on diverse timescales, and only use the Lyapunov optimization framework in small time slots to conduct delay-aware task routing in a distributed manner. In our proposed scheme, the uses of clean power, power storage and transmission delay are all considered.

B. OUR CONTRIBUTIONS
Our main objective is to ameliorate the situation of power consumption and task scheduling in geo-distributed edgecloud systems. To this aim, we propose a two-timescale scheme of eco-friendly powering and delay-aware task scheduling for edge-cloud systems in a deregulated electricity market. Specifically, (1) in large time slots, DCs optimize the number of active servers aiming to realize the trade-off between task latency and power costs, using a weighted sum of monetary cost and pollution cost. This power-latency trade-off is formulated as problem P and can be efficiently solved by our proposed Sequential Convex Programming (SCP) algorithm. We also introduce a convex pollution indicator function (PIF) to measure the pollution of the powers supplied by different ECs, which can encourage the uses of cleaner power and improve power savings. (2) In small time slots, called sub-slots, EGs separately optimize their strategies of task dispatching with the observation of capacities and queue backlogs of DCs. A Lyapunov optimization framework is applied for the design of a distributed task scheduling mechanism. To the best of the authors' knowledge, this is the first time the number of active servers and the task scheduling strategies are optimized on diverse timescales, which can efficiently alleviate the conflict between frequent switches of servers and the uncertainty of task arrivals in DCs. After that, we investigate the distributed algorithm design for capacity planning of DCs and task routing. It is worth noting that our scheme is a comprehensive approach, as shown by the illustrative comparisons of the proposed scheme with the current state-of-the-art summarized in Table 1.
Our main contributions can be summarized as follows.
• We propose a two-timescale framework for the capacity planning of geo-distributed DCs and task routing from edge-ends to cloud-ends. Benefiting from the diverse-timescale optimization of DC capacities and task routing, our scheme can provide stable and energy-efficient task latency guarantees.
• We formulate the long-term eco-friendly and delayaware capacity planning problem of DCs in a deregulated electricity market as problem P, and propose a low-complexity SCP algorithm to find the globally optimal non-integer solution of problem P.
• We introduce a convex PIF to measure the pollution cost of powers supplied by different ECs, by which our proposed strategies can improve the usage of clean power by up to 50%-60% of the total bought power in DCs, as well as encourage power savings and multiple-source power buying in a deregulated market. The rest of this paper is organized as follows. First, we introduce the system architecture and give a formulation of the two-timescale framework in Section II. Second, we discuss the solution and the distributed algorithm design in Section III. Third, we give simulation results about the performance of our model in Section IV. Finally, we conclude the paper in Section V. For clarity, abbreviations and notations frequently used in this paper are listed in Table 2 and Table 3, respectively.

A. OVERVIEW OF THE SYSTEM MODEL
We consider an edge-cloud system consisting of I geographically separate DCs and J EGs, as shown in Fig. 1, where J I . EGs are responsible for collecting requests from users and then dispatching tasks to local edge servers or distant DCs. Although the computation capacities of DCs are enormous, the remote distances between cloud-end and edge-end usually result in large transmission latency. Additionally, for the purpose of power conservation, we assume that only some of the servers in DCs are activated according to the instantaneous workload. We regard this edge-cloud system as a time-discrete system. To reduce the system power cost and task latency, we will develop a two-timescale framework, where DCs make decisions on the number of active servers at the beginning of large time slots, and then EGs can adjust their task dispatching strategies in each sub-slots according to the planned available capacities of DCs and the random request arrivals in real time. In this paper, we only concentrate on task routing from EGs VOLUME 8, 2020 to DCs rather than on the inner scheduling of various tasks in DCs (such as workflows [26], delay-sensitive tasks, etc.).
Let DC i , i ∈ I, represent the i-th DC, where I = {1, . . . , I }. All DCs are assumed to be equipped with smart energy controllers for the coordination of buying power from the market, charging or discharging batteries, and supplying power for loads. Specifically, we assume that DC i can buy power from multiple ECs in a deregulated electricity market. Let N i represent the set of candidate ECs for DC i , and q i,n denote the amount of power bought from the n-th candidate EC by DC i , where q i,n ≥ 0, n ∈ N i and |N i | = N i . Besides, we let Q bat i > 0 represent the amount of power consumed to charge batteries, Q bat i < 0 represent the power amount discharged from the batteries, and Q bat i = 0 represent that no power is charged or discharged. Then, the total power that DC i needs to buy can be calculated by where Q cons i is the power consumption of DC i and E{Q rew i } is the expected amount of renewable power generated by the green micro-grid of DC i . Given that the power consumption of DCs is very large in general, while the discharged power and generated renewable power are usually limited, we typically have Q buy i ≥ 0. Denote M i and m act i as the number of total servers and active servers of DC i in the considered time slot, respectively. According to [27], the power consumption of DC i can be approximately calculated as where T is the length of one large slot, s i and its index α are parameters related to the CPU frequency, 1 β i represents the basic power consumption of servers, and ϕ i > 0 is a constant to indicate the fraction of cooling power to servers' power in DC i . Given that the price and unit pollution cost of electricity usually vary in different ECs and in different regions, we need to optimize the number of active servers in geo-distributed DCs aiming at power cost savings in the whole system. On the other hand, let EG j , j ∈ J represent the j-th EG, J = {1, . . . , J }. Without loss of generality, we assume that the request arrivals to different EGs are independent stochastic progresses. Dividing the considered time slot into H subslots, each of equal length T /H , let the random variable l h j represent the request rate (per sub-slot) to EG j in sub-slot h ∈ {1, · · · , H }. We assume that EG j dispatches tasks to DC i with probability φ h ji on sub-slot h, where I i=1 φ h ji = 1, j ∈ {1, · · · , J } and h ∈ {1, · · · , H }. Denoting φ h j = (φ 1 ji , · · · , φ H ji ), we emphasize that EGs should make their optimal decisions on φ h j in each sub-slot according to the 1 Although we ignore the differences among servers in the same DC, this will not affect the generality of our model, which can be readily extended.
real-time information of task arrivals to EGs and queue backlogs of DCs in order to minimize the related task latency.

B. LARGE-SLOT LAYER: CAPACITY PROVISIONING OF DCs
Although less active servers contribute to more power savings, this has the risk of increasing the task latency, so that decisions on the number of active servers should take into account the trade-off between power cost and task latency.

1) POWER COST
It is noted that the power cost considered in this paper includes both pollution cost and monetary cost.
First, we use a pollution indicator function (PIF) f PI (·) to measure the pollution cost of the powers supplied by different ECs. In order to encourage eco-friendly power consumption, the PIF is assumed to meet the following conditions: (i) the more polluting the power, the larger the pollution cost, and (ii) the more the power bought, the larger the total and unit pollution cost. The simplest case of PIF is a quadratic function as where a i,n = γ i,n TP max i is a positive coefficient, γ i,n is the pollution factor of the power supplied by the n-th EC in the region of DC i , and P max i is the maximum power of DC i . We assume that pollution factors for different types of powers are proportional to their pollutant emissions. Then, the pollution factor of an electricity company depends on the proportions of different power types in its power structure. It is noted that a more polluting power always leads to a larger γ i,n . Introducing the divisor TP max i makes the pollution cost of DC i penalized by the proportion of the purchased amount q i,n with respect to its maximum power consumption TP max i . Second, we will address the monetary cost. Let p i,n andε i denote the electricity price published by the n-th EC and the estimated price of stored power in the region of DC i , respectively. The reason for considering storage price is that the stored power and bought power are homogeneous commodities. Present decisions on charging or discharging will influence the future purchases of power [28], so that the storage price is a kind of potential cost. Let i denote the variation of the energy stored in the batteries during the (dis)charging procedure. When i ≥ 0, batteries charge. Otherwise, batteries discharge. Then, the monetary cost of DC i can be given by When we consider the efficiency of charging and discharging, noted as η i ∈ (0, 1), the relationships between i and the (dis)charged power amount Q bat i can be expressed by If we make where η i > 0. Fig. 2a shows an example of η i [29] and the corresponding η i , where δ i = i TC i is the C rate of batteries' charging or discharging and C i is the battery capacity of DC i . We note that η i is a non-linear function of δ i . 2 Based on the above, the total power cost of DC i can be formulated as the sum of pollution cost and monetary cost where q i is the vector of q i,n , n = 1, · · · , N i . The required constraints are where c i is the current battery state of DC i . (8a) is obtained from (1) and (6). (8b) represents that batteries can discharge no more than the stored power and can charge no more than the remaining capacity. In addition, due to the fact that an excessively high speed of charge and discharge will cause severe damage to storage devices as well as exorbitant wastes of energy [30], we set lower and upper bounds for i in (8c) as lb i = −TC i and ub J j=1 l h j as the average task arrival rate of the whole system in the considered large time-slot, and let L σ represent the quantile satisfying Pr(l < L σ ) = σ , where σ is a predetermined probability threshold. As a latency guarantee, we consider the constraint of I i=1λ i ≥ L σ . We assume that the request arrival process to a DC can be modeled as a Poisson process, which is reasonable if the request arrivals are independent and stationary, and then we can use an M/M/n queuing model to estimate the queuing delay in DCs. Although some researches show that the dependency among request arrivals is likely to result in a deviation from the Poisson distribution, it is still useful in the area of queuing theory for its generality and tractability. In this part, our main purpose is to optimize the number of active servers, which can tolerate the estimation errors of queuing delay to some degree, so that we use the assumption of Poisson distribution here. In the sub-slot layer of task scheduling, a more precise approach will be implemented. According to the M/M/n queuing model [7], the task latency in DC i is given by where T /H is used for unit conversion.

3) JOINT OPTIMIZATION OF POWER COST AND TASK LATENCY
Incorporating the power cost in Eq. (7) and the latency cost in Eq. (9), the cost function of DC i , i ∈ I, can be formulated as where θ 1,i > 0 and θ 2,i > 0 are weight parameters. Then, the total cost function of the whole DC population can be given by = I i=1 i . Besides, considering that i is a monotonically increasing function ofλ i , ∀i ∈ I, the inequal- Based on all of the above, an eco-friendly and delay-aware power cost minimization problem in the large-slot layer, denoted as problem P, can be formulated as follows: where Q cons

C. SUB-SLOT LAYER: DISTRIBUTED TASK DISPATCHING
We focus on a mixed strategy of task dispatching in which EG j dispatches tasks to DCs 3 with probability φ h j in sub-slot h. A distributed mechanism of the mixed strategy generation based on the Lyapunov optimization framework is proposed in this part. Specifically, the task latency considered by EGs mainly includes the communication time from EGs to DCs and the sojourn time in DCs. We use the queue backlogs of DCs to measure the sojourn time in DCs and build the Lyapunov function associated with the queue backlogs. Then, we formulate the drift-plus-penalty function as the weighted sum of Lyapunov drift and communication latency. Each EG makes decisions on the mixed strategy of task dispatching separately by minimizing the upper bound of their own driftplus-penalty function on each sub-slot.

1) COMMUNICATION LATENCY
The communication latency mainly results from the transmission latency and propagation time which can not be ignored due to far distances between EGs and DCs in general. Let W h ji (in bits/s) and τ ji represent the effective bandwidth and the one-way propagation time between EG j and DC i , respectively. It is noted that W h ji changes over slots in the sense that the network congestion depends on the traffic state. Then, the communication latency between EG j and DC i in sub-slot h can be given by whereb h j is the average length (in bits) of tasks from EG j in sub-slot h. Considering that the data coming back from DCs is generally much less than that sent out by EGs, we ignore the transmission delay of back-stream from DCs to EGs.

2) SOJOURN TIME IN DCs
All tasks to a DC should be first put in a queue waiting for service with the principle of First Come First Served (FCFS) if no idle server is available immediately, where the queue backlog of a DC is proportional to the sojourn time in it. Let Q h i represent the queue backlog of DC i in sub-slot h. Suppose the service process of DC i in the considered time slot is stationary with average rateμ i , which is independent of the task arrivals and the queue backlog of DC i , then the queue evolution over time can be modeled as where λ h i is the total arrival rate of tasks (per sub-slot) to DC i in sub-slot h. (13) can 3 For clarity, the task dispatching to local edge servers is not incorporated, but our model can be easily extended. be rewritten by Our consideration of incorporating latency minimization with queue stability leads us to a design approach based on the Lyapunov optimization framework [31], which has been widely used in the stability control of dynamic systems. We define a Lyapunov function L(Q h ) as the sum of the squares of queue backlogs in DCs in sub-slot h where Q h is the vector of Q h i . L(Q h ) is a scalar measure of network congestion in the sense that it is small if all queues are small, and it is large if one or more queue is large. The difference in Lyapunov function from one sub-slot to the next and its upper bound can be given by . Then, we define the Lyapunov drift (Q h−1 ) as the expected difference in Lyapunov function over one sub-slot given the current queue backlog Q h−1 : With the hypothesis that the service rateμ i is independent of the queue backlog Q h−1 i , the upper bound of (Q h−1 ) can be obtained by Taking into account the mixed strategy φ h j , tasks dispatched from EG j to DC i can be expressed by λ h ji = l h j φ h ji . Let λ h −ji = m =j λ h mi represent tasks dispatched to DC i by EGs other than EG j , and then there is , which depends on actions of other EGs. To decouple it in a tractable manner, we estimate λ h −ji by historical data, e.g., parameter constrained by h−1 m=1 w m = 1. Besides, with the assumption that the task arrival rate l h j to EG j is independent of Q h−1 , Eq.(18) can be rewritten as

3) JOINT OPTIMIZATION OF COMMUNICATION LATENCY AND SOJOURN TIME
We build a drift-plus-penalty function for EG j as ji } by incorporating the communication latency and Lyapunov drift, where w j > 0 is a weighted parameter, j ∈ J . The upper bound of Following the design principle of Lyapunov optimization framework, we can set our objective to minimize the upper bound of the drift-plus-penalty term shown in (20) in each sub-slot. Ignoring the constant U h j , the task dispatching problem of EG j , j ∈ J based on the Lyapunov optimization framework can be formulated as problem Q j : Apparently, problem Q j is a quadratic programming problem, which can be easily solved by some standard algorithms, such as the interior-point method or the sequential quadratic programming (SQP) method.

III. TWO-TIMESCALE FRAMEWORK OF CAPACITY PLANNING AND TASK DISPATCHING
Problem P in (11) is non-convex, due to the non-linear equality constraint (11b) and the integer restriction for m act i [32], while problem Q j in (21) is a typical quadratic programming problem. Before giving the distributed algorithm for the DC capacity planning and task routing from EGs to DCs, we intend to study the solution of problem P.

A. SOLUTIONS OF PROBLEM P
In this part, we first apply continuity relaxation to problem P and propose an SCP algorithm to find its optimal non-integer solution. Without special statement, 'solution' mentioned in this subpart means non-integer solution. We find that the nonlinearity of constraint (11b) results from the non-linear function η i with respect to i , which is tightly associated with the part of power cost F PC i in . Therefore, we intend to solve F PC i (·) and plug the solution into (11) to eliminate variables q i,n , i ∈ I, n ∈ N i . Then, the continuity-relaxed problem P can be transformed into a sequence of convex problems. This process is summarized in the proposed SCP algorithm. It is shown that the mixed-integer solution derived from the optimal continuous solution is quasi-optimal.
First, ignoring the inequality constraint q i 0, the partial Lagrange function of F PC i with respect to q i constrained by (8a) can be established as where u L i is the Lagrange factor. According to the Lagrange Multiplier method [32], we set ∂L i ∂q i,n = 0 and ∂L i ∂u L i = 0 for ∀n ∈ N i . Then, we can obtain the partial Lagrange dual The corresponding function value can be given by where Then, the marginal power cost v * L i of DC i when q i,n = q * L i,n which is independent of n. Note that where q * L i is given. Then, we formulate a new problem which is denoted as problem P1. The feasible region of continuity-relaxed problem P1 is a convex set on account of the elimination of the nonlinear equality constraint (11b). In addition, we prove in Appendix B that the continuity-relaxed objective function of P1 is convex when the fitted function of η i meets the condition shown in Proposition 1. We observe from experiments that some cubic functions and exponential functions can both fit well the curve of η i shown in Fig. 2, under the convexity condition described in Proposition 1. Fig. 3 gives the optimal fitting curves of η i by cubic function and exponential function, respectively. In this sense, we can regard the continuity-relaxed problem P1 as a convex programming problem [32]. Third, we reconsider the inequality constraint q i,n ≥ 0. It is easily inferred that if all q * L i,n are non-negative, then they are the exactly optimal solutions to given i , i ∈ I, in which case problems P and P1 are equivalent. For ∀n ∈ N i , ∀i ∈ I, if q * L i,n < 0, then the optimal non-negative solution to the purchase amount of power q i,n is zero, which is proved in Appendix A. Under this principle, we propose the SCP algorithm to solve the global optimal solution of problem P. As is shown in Algorithm 1, lines 6 to 13 are used to revise the values of q i,n by updating the set N − i and making all q i,n = 0   (28) and obtain a new version of problem P1, which can be solved by standard convex programming methods, such as the SQP method shown in Line 4. In Appendix A, we give mathematical proofs that our proposed SCP algorithm can always find the globally optimal solution of problem P. Besides, the SCP algorithm converges fast, as it can find the optimal solution of problem P by solving problem P1 no more than 2 × I i=1 N i times. Finally, we consider the integer restrictions of m act i by rounding all m act i to the nearest integers. Simulation results show that the gaps between the function values of with respect to the optimal non-integer solution and this derived mixed-integer solution are on the order of 0.01%, which is negligible.

B. DISTRIBUTED SOLUTIONS OF PROBLEM P
Although problem Q j can be solved by each EG separately, problem P is associated with the whole DC population, so that we should first study the decomposition of problem P. We find that the coupling among DCs in problem P results from the linear equality constraint (11a), which can be easily decoupled by the method of partial Lagrange dual [32]. First, denoting π as the Lagrange multiplier associated with constraint (11a), the partial Lagrange function of problem P is given by Calculate X i , Y i and Z i on subset N − i according to (25) for ∀i ∈ I, and plug them into problem P1. 4: Solve problem P1 with a convex programming method, such as SQP.

7:
If there is any q * L i,n < 0, n ∈ N − i then 8: Move all indices n of q * L i,n < 0 from N − i to N − i 9:

10:
If there is any p i,n < ν * L i , n ∈ N − i then 11: Move all indices n of p i,n < ν * L i back to N − i 12:

End For
Then, the corresponding partial Lagrange dual function is where the subproblem D i (π) represents It is noted that D i (π) is only related to the individual decision variables of DC i , so that subproblem (31) can be separately solved by DC i , i ∈ I in a parallel manner. Furthermore, the partial dual problem of problem P can be formulated by From the perspective of the system operator, we can update the value of π in an iterative manner according to the method of gradient descent: where t = 1, 2, · · · is the count of iterations, l s > 0 is the step length and π(t) represents the updated value of π in the t−th iteration. Next, to make our proposed SCP algorithm valid for solving D i (π) shown in (31), the subproblem of the dual problem (32), we conduct the same operations described in subsection 3.1 on D i (π ) to eliminate q i , so that subject to the boundary restrictions (11d)-(11f). i − πλ i is convex given π, since the summation of a convex function and a linear function is still convex [32]. Then, for a given π, we can solve D i (π) by our proposed SCP algorithm as long as the set I in the SCP algorithm is initialized as {i}. In each iteration round of the SCP algorithm, an updated version of problem (34) will be solved rather than problem P1. We emphasize that although problem P1 is convex, the original problem P is non-convex, so that the optimal solution to the dual problem (32) can only provide a lower bound for the original problem P. The gap between the optimal function values of the dual problem (32) and the original problem P is regarded as the cost of distributed algorithm design. Fig. 4 shows the two-timescale framework of the distributed capacity planning and task dispatching algorithm. The procedure is as follows. (i) Large-slot layer. At the beginning of the considered large slot, the cloud layer should optimize the strategies of power purchasing and make decisions on the number of active servers by solving the partial dual problem (32) of problem P. The distributed solution procedure of problem (32) is iterative as follows.

C. TWO-TIMESCALE FRAMEWORK DESIGN
• Initialize the iteration count t = 1, the operating factor π(1) = 1, etc. Then, repeat the following three steps: -DC i , ∀i ∈ I receives the updated value of operating factor π(t), and then solves its individual VOLUME 8, 2020 problem (31) by the SCP algorithm. After that, each DC sends its updated solutionλ i to the coordinator. -Based on the feedbackλ i from each DC, the coordinator updates π(t + 1) according to (33). -The coordinator checks the termination condition.
If it holds, then this iteration procedure will be finished. Otherwise, the coordinator will broadcast the updated operating factor π(t + 1) to all DCs and make t = t + 1.
• Output the power purchasing strategy and planned number of active servers for each DC. It is noted that this iterative procedure is run by the coordinator. The termination condition referred to in the repeated round is that either the maximum iteration count has been reached or the descent distance of operation factor π is lower than a predetermined threshold.
For clarity, we give some necessary illustrations of the above procedure. Problem (32) is the partial dual problem of the original problem P, and can be solved in a distributed manner by iteratively solving its two subproblems (31) and (33). When we solve sub-problem (31) by our proposed SCP algorithm, an updated version of problem (34) will be involved in each iteration round of the SCP algorithm. Actually, the original problem P could be directly solved in a centralized manner, where a sequence of problems P1 will be involved in the SCP algorithm rather than problem (34). However, the centralized mechanism requires large amounts of individual information of the DCs to be published, such as the storage capacity, computation capacity, etc., which violates the principle of privacy protection, especially when the involved DCs are affiliated with different service providers.
(ii) Sub-slot layer. At first, EG j , j ∈ J obtains the number of active servers m act i of DC i , i ∈ I. Then, in each subslot, the EG population conducts the following two steps in a parallel manner.
• EG j , j ∈ J observes the processing speed µ h i and queue backlog Q h−1 i of DC i , and estimates the task arrivalsλ h −ji to DC i according to historical data, ∀i ∈ I.
• Accordingly, EG j , j ∈ J optimizes its dispatching strategy φ h j by solving problem Q j , and then dispatches tasks according to the optimized φ h j . We assume that DCs can dynamically control CPU frequencies of servers according to real-time queue backlogs of themselves, so that EGs should update the knowledge of the processing speed µ h i of DC i in each sub-slot, ∀i ∈ I. The cooperation of the dynamical CPU-frequency mechanism and the Lyapunov optimization technique can further improve the stability of queues in DCs, which brings significant benefits in terms of guarantees for the task latency. According to reference [18], the optimal strategy of a DC in the situation of dynamical CPU frequency is that all active servers have identical CPU frequencies. Due to space limitations, we will not elaborate on the dynamical CPU-frequency mechanism of DCs, which is an issue independent of our main concerns.
So far, we have finished the discussion on two-timescale framework design. It is observed that only little information needs to be exchanged among agents in the above distributed mechanism, including π(t),λ i (t), m act i , µ h i , λ h−1 i and Q h−1 i , which is beneficial for the protection of privacy and autonomy of agents. Besides, it can reduce communication costs in the control procedure and difficulties in the implementation. Although iterative updates of the operation factor π in the distributed scheme always consume extra communication costs compared to the centralized scheme, this is acceptable when the length of the large slot is large enough (e.g. one hour).
Finally, we briefly analyze the complexity of our proposed scheme. We roughly classify the complexity into three levels according to the adopted optimization technique and the problem scale, i.e., low, medium and high complexity. Table 4 gives a complexity analysis of some related state-of-the-art approaches. Specifically, the Lyapunov optimization technique is typically used in real-time system control, and the related schemes are of low complexity only if their real-time control problems are simple. The sub-slot scheme proposed in our work and those in [23] [24] only refer to small-scale linear or quadratic programming problems, so that all of them can be regarded as low-complexity approaches. Additionally, in the large-slot layer, the properties and solution approach (SCP algorithm) of problem P are similar to those in reference [17], while the problem scale is smaller. Thus, the largeslot scheme proposed in our work has lower complexity than that in [17]. Moreover, the problem scale of our proposed sub-slot scheme and large-slot scheme depends linearly on I × J and I , respectively, where I is the number of DCs and J is the number of EGs.

IV. NUMERICAL ANALYSIS
To verify our model, we use Matlab to solve problem P by our proposed SCP algorithm and give the performance analysis focusing on four aspects, including (i) the impact of PIF, (ii) the behavior of storage scheduling, (iii) the analysis of task dispatching, and (iv) the reduction of power cost and the improvement of clean power usage. Table 5 gives default settings of some parameters. Unless specified otherwise, we configure the parameters of all of our experiments according to Table 5. For example, the benchmarks of the electricity price used in the following experiments correspond to the default electricity price shown in Table 5. Specifically, we assume that each DC can buy power from N i = 3 ECs, whose power for sale mainly originates from thermal power (TP), solar power (SP) and wind power (WP), respectively. The corresponding electricity prices 4 always change at different hours and in different areas.

A. EFFECT OF THE POLLUTION INDICATOR FUNCTION
First, we aim to verify the effect of PIF on multi-source power buying. Fig. 5a shows the unit power costs for a 1-MW data center when buying power from TP, SP, WP alone and from multiple ECs in one time slot, respectively. It is observed that the unit cost increases with the growth of Q buy i . This means that the more the power we buy, the higher the unit cost we will pay, which can lead to an increasingly fast growth of total power costs and encourage power savings in order to reduce costs. In addition, we can see that the unit costs in the cases of buying power from TP, SP or WP alone are always higher than those when buying multi-source power, which is caused by the quadratic term in PIF. Therefore, our proposed PIF can encourage users to buy power from multiple suppliers rather than single suppliers. Fig. 5b shows the ratio of clean power, including SP and WP, to the total bought power Q buy i . It can be observed that 4 All electricity prices used in this paper are set according to those in China [33], [34], and converted from to $ at an exchange rate of 0.148. the fraction of clean power will increase with the growth of Q buy i when the electricity prices p i,n , ∀n ∈ N i are different, which means that the more the power we buy, the higher the fraction of clean power we will use, because the unit cost of clean power increases more slowly than that of brown power due to the smaller pollution factor γ i,n , although the price p i,n of clean power is higher in general. However, when p i,n , ∀n ∈ N i are identical, the ratio of clean power is a constant and is independent of Q buy i . In this case, the numerators of q * L i,n ∀n ∈ N i in Eq. (23) are equal, thus the ratio q * L i,1 : · · · : q * L i,N i can be calculated as 1 a i,1 : · · · : 1 a i,N i , which is constant.
Furthermore, we can see that the fraction of clean power is sensitive to the configuration of the electricity price and the pollution factor γ i,n . The fraction will be higher when the differences among γ i,n increase, or the differences among p i,n decrease, ∀n ∈ N i . We can conclude from Fig. 5 that our proposed PIF can encourage power savings, power buying from multiple sources and usage of clean power. The more the power we buy, the higher the usage ratio of clean power. If we set the unit pollution cost to a constant as in the previous works, then the total pollution cost will be a linear function rather than a non-linear function as PIF. In the linear case, users will always choose the EC with lowest unit cost, the sum of electricity price and unit pollution cost, and abandon the powers from other ECs that are a little more expensive. Different from the linear case, the PIF can motivate users to buy power from multiple ECs and the ratio of clean power will increase with VOLUME 8, 2020 the growth of the total bought power. In addition, the fraction of clean power is sensitive to several parameters, such as electricity price and pollution factor γ i,n .

B. ANALYSIS OF THE STORAGE SCHEDULING
We first fit the curve of η i shown in Fig. 2b with many different power and exponential functions. We find that the cubic function 0.873δ 3 + 1.830δ 2 + 1.495δ + 1.038 where δ ∈ [−1, 0.3] is the simplest one that meets the condition in Proposition 1 to guarantee the convexity of problem P1 with negligible fitting errors, so that we adopt the above cubic function in this paper.
Second, we aim to verify the impact of potential cost on storage scheduling and power cost. Fig. 6a shows the optimal decisions on charge and discharge without considering potential cost, in which case batteries discharge much in the first time slot and do not charge or discharge in the remaining slots in order to obtain the maximum present interests while ignoring the future interests. Fig. 6b shows the optimal decisions on charge and discharge considering potential cost, where batteries always tend to discharge in slots with higher unit cost and charge in slots with lower unit cost, where h is the index of a future slot. This reveals that the long-term storage scheduling according to the varying unit-cost of storage can only be optimized when the potential cost is taken into consideration. Additionally, we note that the potential cost can actually influence the storage scheduling in a duration longer than hours. It is observed from our many experiments that the storage scheduling controlled by potential cost generally FIGURE 6. Impact of potential cost on storage scheduling and power cost. exhibits a periodicity over days, following the corresponding electricity price. To clearly show how the fluctuations of the electricity price influence (dis)charging actions, we just show the storage scheduling during a day. Fig. 7a shows the daily monetary costs for a 1-MW DC with different accuracy ofε i and that withoutε i , where |ê| represents the prediction errors forε i measured as the differences between the actual values and the predicted values. The referred benchmark corresponds to the realization ofε i where |ê| = 0. It is observed that introducing potential cost can improve the power savings, but a small prediction errorê for ε i can cause a severe performance decline of storage scheduling. This means that the performance of storage scheduling is very sensitive to the parameterε i , and an accurate prediction method of the unit cost of storage in the future slots is needed, which will be further studied in our future work. Third, we aim to illustrate how the electricity price influences the power savings of storage scheduling, as shown in Fig. 7b. We adopt the default electricity price as a benchmark, and change the average values 'AVE' and the variance 'VAR' of the electricity price, and the square sum of differences between prices in adjacent time slots 'DIF 2 ', respectively, where '+' represents increase and '−' represents decrease. It is found that the ratios of monetary costs saved by storage scheduling to the total monetary costs will increase when 'VAR' or 'DIF 2 ' grows. This means that the higher the differences among electricity prices in different slots, 96480 VOLUME 8, 2020 the more the monetary costs saved by storage scheduling. However, the ratio of cost savings decreases when 'AVE' increases, because when the electricity prices in different slots are equally increased, the differences among them do not change and the monetary costs with and without storage scheduling increase equally. increases with the decrease ofv i , which means that more servers will be activated in areas with lower electricity cost in order to reduce power costs. It is observed thatλ i is proportional to m act i . Given the selfishness of players in distributed mechanisms, EGs always prefer low-latency generating DCs, so that the average request rate to DCs in a long duration can be balanced. However, if we observe in real time, periodical fluctuations and oscillations will be found, which is also caused by the selfishness of EGs along with the lack of cooperation among EGs. Benefiting from the form of squared costs in Lyapunov optimization techniques, the Price of Anarchy (PoA) of our proposed distributed task dispatching mechanism is maintained within a reasonable range, as shown in Fig. 9. PoA is a widely used scalar measure of the system performance degradation caused by the anarchy of players. In this paper, we define PoA as the ratio of task latencies derived from our proposed distributed mechanism and from the centralized mechanism. Besides, we can also see from Fig. 9 that the PoA increases with the growth of the number of EGs. Table 6 shows the monetary savings of the whole DC system in one time slot achieved by workload scheduling with different configurations of electricity price p i,n , number of DCs I and total request rate L σ , where the default parameters shown in Table 5 are used in the benchmark case and L max is the maximum request rate that the system can deal with to avoid overload. It is observed that the larger the geographical differences in the electricity price, measured by the variances, the more the monetary savings. However, when the prices in different areas increase equally, in which case only the means increase while the variances remain constant, monetary savings will not increase. That is to say the geographical difference in the electricity price is the main contributor to the power savings gained by workload scheduling. In addition, monetary savings will also increase with the growth of I and L σ , which means that the effect of workload scheduling on power cost savings will be improved when the number of geographical DCs or the total request rate increases.  Table 7 shows the average request delay of the whole DC system in the case of different numbers of DCs I and different weight parameters θ 1,i , where the upper-bound queuing delay D is set in terms of the web requests as 2 s. It can be observed that our formulation will never make the task latency approach or exceed the upper bound that users can tolerate, which is different from the approaches presented in [10]- [25]. That is to say our proposed method can realize the trade-off of power cost and delay cost by a flexible   weight parameter because the request delay is a subpart of the objective function. Fig. 10 and Fig. 11 show comparison experiments of our proposed scheme to another two-timescale scheme investigated by Yao [24], [25] and to a typical single-timescale scheme investigated by Wang [18]. Specifically, we make observations of the average sojourn time experienced by tasks in DCs and average power of DCs in different situations of workload density, measured by the fraction of the actual request rate and by the system capacity. We can see that Yao's scheme has better latency performance than ours when the workload density is low, but its latency performance exhibits a rapid degradation with the increase of the workload density. On the other hand, our scheme exhibits a slower growth of task latency and is much better than Yao's scheme in situations of medium and high workload densities. This phenomenon mainly benefits from our diverse-timescale decisions on the number of active servers and on task routing, which can better adapt to the uncertainty of task arrivals than Yao's scheme. Additionally, it is observed that the latency performance of Wang's scheme exhibits obvious fluctuations. This is because the task routing strategy and the CPU-cycle frequency can not be adjusted in time to adapt to the time-varying workload density. We note that dynamical CPU-frequency strategies are conducted in small timescales in both Yao's scheme and our scheme. Moreover, the power consumption performance of our scheme is between those of Yao's and Wang's.  Table 8 shows the monetary cost savings obtained by the joint optimization of workload scheduling and power storage scheduling in problem P. It is observed that the monetary cost savings obtained by the joint optimization are larger than those obtained by workload scheduling in Table 6 or power storage scheduling in Fig. 7b alone, and will increase with the growth of I , L σ and the variance of the electricity price. This means that (i) the joint scheduling of workload and power storage can further improve the power cost savings, and (ii) the amount of cost savings is directly proportional to the number of DCs, the request rate, and the spatial and temporal change rate of electricity price. Specifically, we find that cost savings of up to 10%-30% of the total monetary costs can be obtained by using workload and storage scheduling. However, this kind of power cost reduction depends on the consideration of spatial and temporal differences of unit power cost, which has nearly no social benefits if the unit cost is decoupled from clean power generation [13]. Therefore, we will next check the efficiency of our proposed scheme from the point of view of clean power usage.   Table 9 shows the corresponding fractions of bought clean power. It is observed that when introducing the pollution factors γ i,n or increasing the differences among the γ i,n of different ECs, the monetary costs grow, while the fraction of clean power is improved. This means that our proposed method can realize a trade-off between monetary costs and pollution costs, by which the fraction of clean power can be improved to 50%-60% of the total bought power.

V. CONCLUSION AND FUTURE WORK
In this paper, we focused on the eco-friendly and delay-aware power cost minimization of geo-distributed DCs in a deregulated electricity market. We formulated this problem as problem P and proposed the SCP algorithm to obtain the globally optimal non-integer solution. Based on this, we obtained the quasi-optimal mixed-integer solution of problem P. Simulation results revealed that our method could effectively cut down the total power cost and encourage an eco-friendly use of power, as well as reduce the delay of requests, achieving monetary cost savings of up to 10%-30%. More importantly, our proposed PIF-based power cost model can greatly improve the use of cleaner power and encourage power saving. For example, for a 1-MW DC, the amount of clean power it uses can reach 60% of the total power it buys.
Some interesting improvements toward the adjustment strategies of the integer-constrained variables, not reported here due to space constraints, can be found in [35]. As part of our future work, we will explore better methods to predict the future unit cost of power storage and distributed algorithms to realize workload scheduling under transmission delay considerations.

APPENDIX A PROOFS ABOUT THE OPTIMALITY OF THE SCP ALGORITHM
In Proposition 1, we prove that the optimal solution for those q i,n with q * L i,n < 0 in problem (7) is zero. Based on this, we prove that our proposed SCP algorithm can obtain the globally optimal solution of problem P via Proposition 2 and Proposition 3.
Proposition 1: Ignoring i , ∀q i,n in problem (7), if the optimal Lagrange solution q * L i,n obtained by Eq. (23) is less than zero, then its optimal solution q * i,n of problem (7) with the constraint q i,n ≥ 0 is zero, where i ∈ I, n ∈ N i .
Proof: First, we assume that there is only one element of the optimal Lagrange solutions obtained by Eq. (23) that is less than zero, whose index is denoted by (i, m), m ∈ N i , i ∈ I, i.e., q i,m is the only variable with q * L i,m < 0. Then it is true that is the optimal Lagrange solution, obtained by Eq. (23) after removing m from N i , with the constraint n∈N i &n =m q i,n = Q cons i + Q bat i − q i,m . Secondly, we denote f n (q i,n ) = a i,n q 2 i,n + p i,n q i,n , and rewrite (7) as F PC i (q i ) = n a i,n q 2 i,n +p i,n q i,n = n f n (q i,n ). Then (A-1) can be obtained as Here, q * L i,n (0) is the q * L i,n obtained by Eq. (23) when q i,m = 0, and α i,n is equal to a i,n X i,m according to (23), where X i,m is X i after eliminating m from N i .
Thirdly, because all f n (·) are strictly increasing and convex, and q * L i,m < 0 < q + i,m where q + i,m represents q i,m > 0, we can obtain (A-2) as follows.
In addition, it can be proved by contradiction that Thus we can obtain (A-3) as below.
Finally, based on(A-1), (A-2) and (A-3), we can obtain (A-4) as This means that the optimal solution for q i,m in (7) whose q * L i,m < 0 with the constraint q i,m ≥ 0 is zero. Based on the above, the case that there is only one q L * i,n < 0 has been proved. Then we can similarly prove the case with two, three or more q L * i,n < 0. Proposition 2: The inner iteration procedure of the SCP algorithm shown in lines 7-9 can obtain the optimal solution of problem (7) for the ith data center.
Proof: First, in the t-th, t ∈ {1, · · · , T } iteration of Algorithm 1, where T is the total number of iterations, the optimal solution of any q i,n which has q L * i,n < 0 is equal to zero according to Proposition 1.
Secondly, ∀q i,n , n ∈ N i whose q L * i,n ≥ 0, the constraint q i,n ≥ 0 is a slack constraint [32], so that its optimal solution in (7) is the q L * i,n obtained by Eq. (23). Thirdly, we need to prove that zero is the final optimal solution of the q i,n set to zero in the t-th iteration, although it is optimal in the t-th iteration. When t = T − 1, the optimal solutions of the remaining q i,n can be obtained in the T th iteration, all of which are no less than zero, and there are no changes of the q i,n set to zero in the (T −1)-th iteration. When 1 ≤ t < T − 1, there is at least one new n added to N − i in the (t +1)-th iteration. We denote q i,n 1 as the one set to zero in the t-th iteration, and then denote q * L,t i,n 1 and ν * L,t i as the optimal Lagrange solution of q i,n 1 and the corresponding marginal cost obtained by (25) before setting q i,n 1 to zero in the t-th iteration. Considering q * L i,n (q + i,m ) < q * L i,n (0) < q * L i,n as stated in the proof of Proposition 1, when we set q i,n 1 = 0 in the t-th iteration, the optimal Lagrange solutions of the remaining q i,n obtained in the (t + 1)-th iteration are all less than those obtained in the t-th iteration. Thus ν * L,t i > ν * L,t+1 i , which can be further extended to ν * L,t 1 i > ν * L,t 2 i , t 1 < t 2 . The marginal cost of q * L,t i,n 1 is ν * L,t i . If q i,n 1 = 0 is not the finally optimal solution, then it would be recalculated as in (23) in the (t + m)-th iteration, m = {2, · · · , T − t − 1}, which must be larger than zero, and the corresponding marginal cost is ν * L,t+m i . Then ν * L,t+m i > ν * L,t i , because q * L,t+m i,n 1 , which is a contradiction. Therefore, we can conclude that zero is the final optimal solution of the q i,n set to zero in the t-th iteration. Based on the above, the complete N − i can be obtained in the (T − 1)-th iteration, in which the optimal solutions of all q i,n are zero, and the optimal solutions of all q i,n ∈ N − i can be solved in the T -th iteration. In addition, it can be easily proved that T ≤ max{N 1 , · · · , N I } + 1, so that it is true that the inner iteration procedure of the SCP algorithm shown in lines 7-9 can obtain the optimal solution of problem (7) for the ith data center with finite repeated rounds.
Proposition 3: Algorithm 1 is guaranteed to obtain the optimal non-integer solution of problem P shown in (11) Proof: Denote the continuity-relaxed problem P and problem P1 as problem P and problem P1 , respectively. Assume that problem P1 is convex, which will be proved in Appendix B. If the optimal values of all q i,n substituted into problem P1 satisfy the non-negativity constraints, then problem P1 is equal to problem P and can be solved by the standard SQP algorithm. Then we will prove that all optimal q i,n ≥ 0 can be obtained by Algorithm 1 as follows.
Different from the individual power cost problem for a single DC in problem (7), for problem P , the final optimal solution of the q i,n set to zero in the t-th iteration may be larger than zero. Thus we add lines 10-12 into Algorithm 1 to move the n whose q L * ,t i,n < 0 but q * i,n > 0 back to N − i . We denote q i,m as the one whose q L * ,t i,n < 0 but q * i,n > 0, denote + > 0 as an infinitesimal positive number, and denote In addition, it can be easily inferred that T ≤ 2 * ( There are two steps leading to the iterations in the SCP algorithm, including (i) setting those q i,n = 0 whose optimal Lagrange solution q * L i,n < 0 temporarily, (ii) recalculating their optimal values if zeros are not the optimal solutions. In the worst case, all q * L i,n are less than zero, which causes In conclusion, Algorithm 1 is guaranteed to obtain the optimal solutions of all q i,n limited by q i,n ≥ 0, so that it is sure to obtain the optimal non-integer solution of problem P .

APPENDIX B CONVEXITY CONDITION FOR PROBLEM P1
In this section, we will give proofs about the condition that η i (δ i ), or η i for short, needs to meet to guarantee the convexity of the continuity-relaxed problem P1 defined on the convex set C. Here, C is a convex set composed of all constraints in (28) except the integer restrictions for m act i . Proposition 4: For any function of η i (δ i ) whose first and second partial derivatives are all existent and continuous, as long as the condition 2 · dη i d i + d 2 η i d( i ) 2 · i ≥ 0 is met in the convex set C, where δ i = i THC i and ∀i ∈ I, then problem P1 is a convex programming problem.
Here, T and H are the length and the number of sub-slots, respectively. Proof: [Proof] For any function of η i (δ i ) whose first and second partial derivatives are all existent and continuous, ∀i ∈ I, we can obtain the first and second partial derivatives of i . Based on this, the Hessian Matrix of i in problem P1 can be easily obtained, and the Leading Principal Minors i,1 , i,2 and i,3 are given as follows.
As we have stated that m act iû i −λ i > 0, Q cons i + η i i ≥ 0, X i > 0, Y i > 0, etc. above, it can be concluded that both i,1 > 0 and i,2 > 0 are true in the convex set C. i,3 ≥ 0 is also true in the case of 2 · i ≥ 0, the Hessian Matrix of i is positive semidefinite in the convex set C. Furthermore, because a sum of convex functions is also a convex function according to [32], is convex in the convex set C and the continuity-relaxed problem P1 is a convex programming problem when 2 · dη i d i + d 2 η i d( i ) 2 · i ≥ 0 , ∀i ∈ I for any η i (δ i ) whose first and second partial derivatives are all existent and continuous.