Minimizing the Maximum Expected Waiting Time in a Periodic Single-Server Queue with a Service-Rate Control

Full terms and conditions of use: https://pubsonline.informs.org/page/terms-and-conditions This article may be used only for the purposes of research, teaching, and/or private study. Commercial use or systematic downloading (by robots or other automatic processes) is prohibited without explicit Publisher approval, unless otherwise noted. For more information, contact permissions@informs.org. The Publisher does not warrant or guarantee the article’s accuracy, completeness, merchantability, fitness for a particular purpose, or non-infringement. Descriptions of, or references to, products or publications, or inclusion of an advertisement in this article, neither constitutes nor implies a guarantee, endorsement, or support of claims made of that product, publication, or service. Copyright © 2019, The Author(s)


Introduction 1.A Nonstationary Stochastic Design Problem
In this paper, we address an open problem from Whitt (2015), which considered the problem of stabilizing performance over time, that is, making a specified time-dependent performance measure a target constant function, in a single-server queue with unlimited waiting space, the first-come, first-served (FCFS) discipline, and a time-varying arrival-rate function.The stabilization is to be achieved with a deterministic service-rate function, under the assumption that the customer service requirements are specified independently of the service-rate control.This is a stochastic design problem instead of a real-time stochastic control problem; that is, the service-rate control is to be determined in advance, assuming full knowledge of the model, but without knowledge of the system state (e.g., the value of the stochastic queue-length process) that will actually prevail at any time.
As explained in section 1 of Whitt (2015), variants of this service-rate control are performed in response to time-varying demand, in many service operations, such as hospital surgery rooms and airport inspection lines, but little is known about the ideal timing and extent of service-rate changes.Service-rate controls for singleserver queues are also of current interest within more complex systems, such as in energy-efficient data centers in cloud computing (Kwon and Gautam 2016) and in business process management (Suriadi et al. 2017).
In Whitt (2015) it was shown that a rate-matching control, whereby the service rate is made proportional to the arrival rate, stabilizes the queue-length process but not the (virtual) waiting-time process.In this paper we develop an algorithm to approximately stabilize the expected waiting time at a target level.It uses a modification of the service-rate control involving two parameters: a time lag and a damping factor.

Related Literature
There is a large literature on similar stochastic design problems involving setting staffing levels (the number of servers) in a multiserver queue to stabilize performance in the face of time-varying demand (e.g., Jennings et al. 1996, Feldman et al. 2008, Stolletz 2008, Liu and Whitt 2012b, Defraeye and van Nieuwenhuyse 2013, He et al. 2016, Pender and Massey 2017, Liu 2018, Whitt 2018).For a single-server queue, the direct analog would be turning on and off the server, which is a restrictive extreme version of the service-rate control we consider.
The dynamic control problem of turning on and off the server in specified system states has received considerable attention in the stationary setting, starting with Yadin and Naor (1963) and Heyman (1968).Similar dynamic control problems for single-server queues including service-rate controls have been analyzed as Markov decision processes in George and Harrison (2001) and Adusumilli and Hasenbein (2010) and references therein.We emphasize that our design problem is different; our service-rate control is for a nonstationary model and must be set in advance, without knowledge of the system state.In many cases, our new problem is more realistic because arrival rates are often strongly time varying and can be reasonably well estimated in advance, whereas changes to the service rate may be difficult to implement without advance planning.Of course, in general, both problems are important.
Given the extensive research on the staffing design problem for many-server queues, it is natural to consider variants of the successful staffing algorithms, but it is now well known that the behavior of many-server queues tends to be dramatically different from that of single-server queues.That difference can be seen by comparing the many-server fluid models in Liu and Whitt (2012a) with the single-server fluid models in Chen and Mandelbaum (1991), as discussed in Liu and Whitt (2011, p. 836).A simple fluid model supporting the rate-matching control in Whitt (2015) is supported by our heavy-traffic weak law of large numbers in Theorem 4 (see Corollary 3), but we are working to go beyond that.
Hence, it should not be surprising that service-rate controls using variants of the established many-server staffing algorithms are no longer effective for single-server queues.For example, a natural analogue of the square-root staffing function from Jennings et al. (1996) was considered as a candidate service-rate control in equation 2.3 of Whitt (2015) but was found to be ineffective, as illustrated by figure 2 of Whitt (2015).Additionally, variants of the iterated staffing algorithm in Feldman et al. (2008) and Defraeye and van Nieuwenhuyse (2013) were found to be ineffective, evidently because the controls have impact over greater time intervals (are less "local") with single-server systems.
As indicated in Whitt (2015), controlling the service rate to meet time-varying demand is analogous to Kleinrock's classic service-capacity-allocation problem in a stationary Markovian Jackson network (Kleinrock 1964); we allocate service capacity over a time instead of over space (different queues within the network).

The Rate-Matching Service-Rate Control
Given that the service requirements are specified independently, the actual service times resulting from a timevarying control are relatively complicated, but a construction is given in section 3.1 of Whitt (2015).In Whitt (2015), several controls were considered, but most attention was given to the rate-matching control, which chooses the service rate to be proportional to the arrival rate; that is, for a given target traffic intensity ρ, the service-rate function is with ≡ denoting equality by definition.In Whitt (2015), theorem 4.2 shows that the rate-matching control stabilizes the queue-length process, theorem 5.1 gives an expression for the waiting time with the ratematching control, whereas theorems 5.2 and 5.3 establish heavy-traffic limits showing that the queue length is asymptotically stable, but the waiting time is not, being asymptotically inversely proportional to the arrivalrate function.

The Open Problem: Stabilizing the Expected Waiting Time
The open problem from Whitt (2015) is developing a service-rate control that can stabilize the expected waiting time.(We only discuss the continuous-time virtual waiting-time process in this paper, which is the waiting time of a potential or hypothetical customer if it were to arrive at that time, and so omit "virtual.") Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control Toward that end, we now study a modification of the rate-matching control.Without loss of generality, we write the periodic arrival-rate function as where 0 < ρ < 1, and s is a periodic function with period c satisfying s ≡ 1 c c 0 s(u) du ≡ 0. (3) As a regularity condition, we require that Most of our numerical examples will be for a sinusoidal function, where s(t) β sin(γt) for s(t) in ( 2), so that and β is the relative amplitude, with 0 ≤ β ≤ 1, and the period is c 2π/γ.In the periodic setting of Equations ( 2)-( 4), we consider the rate-matching control in (1) modified by a time lag η and damping factor ξ; in particular, for 0 < ξ ≤ 1 and η > 0. Thus, the average arrival rate and service rate are λ ρ and μ 1, so the long-run traffic intensity is ρ ≡ λ/ μ ρ.However, the instantaneous traffic intensity ρ(t) ≡ λ(t)/μ(t) can satisfy ρ(t) > 1 for some t in each periodic cycle, for example, if in the setting of Equations ( 5) and ( 6).

Formulation of Optimal Control Problems
Because it is directly of interest, and because we want to allow for imperfect stabilization, we formulate our control problem as minimizing the maximum expected waiting time over a periodic cycle [0, c].We formulate the main optimization problem as a min-max problem; that is, where E[W y ] is the expected (periodic) steady-state (virtual) waiting time starting at time yc within a cycle of length c, 0 ≤ y < c, and }(m) is the set of all periodic service-rate functions with average rate m, which we take to be m ≡ 1.
Given that the average arrival rate is ρ < 1, the obvious reference case is the mean waiting time E[W] in the associated stationary model, which for the M/GI/1 model is and thus 1.For the general periodic problem, what is the solution (value of w * and set of optimal service-rate functions μ * (t) as a function of the model)?
2. For the sinusoidal special case in Equation ( 5), what is the solution?3. To what extent do the optimal solutions stabilize the expected waiting time E[W y ] over time?In particular, is it possible to stabilize E[W y ] perfectly?
Remark 1 (Stabilizing the Full Waiting-Time Distribution).Theorems 4.1 and 4.2 of Whitt (2015) show that ratematching control stabilizes the delay probability P(W y > 0), whereas corollary 5.1 of Whitt (2015) shows that ratematching control cannot stabilize the mean waiting time.Theorem 5.2 of Whitt (2015) establishes a heavy-traffic limit (with periodicity in the stronger fluid scale in Section 6 here) that shows that it is not possible to stabilize the queue-length and waiting-time processes simultaneously.Thus, we conclude that it is not possible to stabilize the full waiting-time distribution.Hence, the open problems above are only for the mean.In this paper we primarily focus on the mean, but we also show that it stabilizes the entire distribution to some extent in Section 9.1.
In this paper we only consider the restricted set of controls in Equation ( 6).Now our goal is For practical purposes, this two-parameter control is appealing for its simplicity.We also find that it is quite effective, even though it cannot stabilize E[W y ] perfectly.
We also consider the associated stabilization control, whereby Equation ( 9) is replaced by In our sinusoidal examples, in which there is strong symmetry, we find that the solutions to Equations ( 9) and ( 10) are the same (but we have no proof), but neither stabilizes perfectly.For more general periodic arrival-rate functions, we detect differences.

Organization of This Paper
This paper involves some challenging technical methods.Hence, we present the more accessible results first.
We start in Section 2 by presenting two simulation examples to illustrate the effectiveness of our new algorithm.Then, in Section 3, we introduce the two technical tools we will apply: (1) an extension of the rareevent simulation algorithm for the GI t /GI/1 model from Ma and Whitt (2018) to the GI t /GI t /1 model with a general service-rate control and (2) heavy-traffic limits involving scaling of the underlying deterministic periodic arrival-rate function.
We start in earnest in Section 4. We elaborate on the model and key processes representing the workload and the waiting time in Section 4. Theorem 1 shows that the rate-matching control stabilizes the workload process as well as the queue-length process.We discuss the extension of the rare-event simulation algorithm from Ma and Whitt (2018) to our setting and its application to perform simulation search in Section 5.In both Sections 4 and 5, we will be brief because we can draw on Whitt (2015) and Ma and Whitt (2018).
We establish our main heavy-traffic limits with periodicity in the stronger fluid scaling [see Equation ( 13)] in Section 6.We present the proof of the main heavy-traffic limit (Theorem 3) in Section 7. We establish heavytraffic limits with periodicity in the weaker diffusion scaling [see Equation ( 14)] in Section 8.
We give simulation examples in Section 9.In Section 9.1, we present simulation results using the fluid scaling in Section 6; in Section 9.2, we present simulation results using the diffusion scaling in Section 8. We draw conclusions in Section 10.

Simulation Examples for the M t /M t /1 Model
To illustrate the effectiveness of our new algorithm, we show results for two simulation examples.We consider the Markovian M t /M t /1 model with the sinusoidal arrival-rate function in Equations ( 2)-(5).The first example has model parameters (ρ, β, γ) (0.8, 0.2, 0.1), so the average arrival rate is ρ 0.8, the average service time is 1, and the cycle length is c 2π/γ 62.8. Figure 1 (left) shows the expected steady-state waiting time E[W y ] together with the corresponding expected workload E[L y ] and the product λ(y)E[W y ], all for 0 ≤ y < 1.The second example on the right differs only by increasing ρ from 0.8 to 0.95. Figure 1 also shows the upper and lower 95% confidence interval bounds for E[L y ] and E[W y ] with black dashed lines, but these can only be seen by zooming in.
Figure 1 shows that the expected waiting time E[W y ] is well stabilized at a value somewhat higher than the expected steady-state waiting time for the stationary M/M/1 model, which is ρ/(1 − ρ) (4 on the left and 19 on the right).The maximum deviation (maximum − minimum) over a cycle is 0.0335 is for ρ 0.8 and 0.4653 for ρ 0.95.Thus, the maximum relative errors are approximately 0.8% for ρ 0.8 and 2.2% for ρ 0.95, clearly adequate for practical applications.Nevertheless, careful simulations and statistical analysis allow us to conclude that it is impossible to stabilize the expected waiting time perfectly with this control.To see the contrasting view with rate-matching control for this same model, see figure 6 of Whitt (2015).(See Ma and Whitt 2015 for more examples.) It is natural to wonder whether there is any order in the optimal controls found for ρ 0.8 and ρ 0.95 in Figure 1.The dependence on ρ is revealed by the main heavy-traffic limit theorem (Theorem 3).Remark 2 (The Cost of Periodicity).The difference between the stable average waiting time in Figure 1 and the value ρ/(1 − ρ) for the stationary model (4 on the left and 19 on the right) might be called the "average cost of periodicity," but we point out that the overall average waiting time with a service-rate control could be much less than in the stationary model.The classic results for the periodic M t /GI/1 queue in Rolski (1981Rolski ( , 1989) ) do not apply because, in general, the service times are neither independent of the arrival process nor independent and identically distributed (i.i.d.); see Example 1.
Example 1 (Small Expected Waiting Times with Periodicity).To illustrate a nonstationary model with a low average expected waiting time, consider the M t /M t /1 model with the two-level arrival-rate function with period c: where δ < c/2, and 1 A is the indicator function of the set A; that is, 1 A (t) 1 if t ∈ A and 0 otherwise.Let the service-rate function be as in Equation ( 6) with η 2δ and ξ 1.Then the number of arrivals in the interval (c/2) − δ, (c/2) + δ [ ) has a Poisson distribution with mean ρc, whereas the number of potential departures in the interval (c/2) + δ, (c/2) + 3δ [ ) has a Poisson distribution with mean c.Thus, for ρ < 1 and c bδ suitably large, the net input over the interval (c/2) − δ, (c/2) + δ [ )is approximately Gaussian with mean −(1 − ρ)bδ and variance (1 + ρ)bδ, which is unlikely to be positive.By choosing δ suitably small and bδ suitably large, subject to specified ρ, we can make the maximum steady-state expected waiting time, and thus the average, approach 0. One way to explain this phenomenon is to observe that the interarrival times and service times will be highly correlated.

The Key Technical Tools
In this section, we discuss the two technical tools that we use.

The Primary Tool: A Simulation Search Algorithm
Our primary tool for finding good (η, ξ) controls is a simulation search algorithm.For that purpose, we extend the rare-event simulation algorithm for the time-varying workload process in the periodic GI t /GI/1 model in Ma and Whitt (2018) to the GI t /GI t /1 model, where the service rate is time varying as well.(The notation GI t means that the process is a deterministic time transformation of a renewal process; see Section 4.) The workload L(t) represents the amount of work in service time in the system at time t, whereas the waiting time can be represented as the first-passage time The waiting time W(t) coincides with the workload L(t) when μ(t) 1 for all t, but not otherwise.
As in Ma and Whitt (2018), the rare-event simulation algorithm calculates the periodic steady-state workload L y and waiting time W y , starting at time yc within a cycle of length c, 0 ≤ y < 1.We use a search over the parameters (η, ξ), as discussed in Section 5, to solve the optimization problems [Equations ( 9) and ( 10)].The search part is relatively elementary because we have only two control parameters.For background on simulation optimization, see Fu (2015) and Jian and Henderson (2015) and the references therein.
The computational complexity for one control vector (η, ξ) is essentially the same as in Ma and Whitt (2018).In particular, the program running time tends to be proportional to the number of replications and number of y values, which for the case ρ 0.8 in Figure 1 were taken to be 40,000 and 40, respectively.This required approximately 100 minutes on a desktop computer.As indicated in section 4.7 of Ma and Whitt (2018), the run time tends to be of order (1 − ρ) −1 , so the cases with high traffic intensity are more challenging.The simulation search is performed in stages, with fewer y values and replications in the early stages but the full long run at the end to confirm performance.

Gaining Additional Insight: Heavy-Traffic Limits
To better understand how the control parameters and performance depend on the model parameters, we establish heavy-traffic (HT) limits, which involve considering a family of models indexed by ρ and letting ρ ↑ 1, drawing on our previous work in Whitt (2014Whitt ( , 2015) ) and Ma and Witt (2018).That previous work shows that the scaling is very important because there are several possibilities.We use the conventional HT scaling of time by (1 − ρ) −2 (usually denoted by n) and space by 1 − ρ (usually denoted by 1/ n √ ), as in chapters 5 and 9 of Whitt (2002), but if we do so without also scaling the arrival-rate function, then the HT limit is easily seen to be the same as if the periodicity were replaced by the constant long-run average, as shown by Falin (1989).
To obtain insight into the periodic dynamics, it is thus important to also scale the arrival-rate function, which is initially specified in Equation (2) with Equations (3) and (4).However, the work of Whitt (2014Whitt ( , 2015) ) actually uses two different HT scalings of the arrival-rate function.Our main HT scaling in Section 6 follows Whitt (2015) and has periodicity in the fluid scale; that is, but in Section 8 we also consider the scaling from Whitt (2014) and Ma and Whitt (2018), which has the periodicity in diffusion scale; that is, The extent of the periodicity is stronger in Equation ( 13) than in Equation ( 14) because of the extra factor (1 − ρ) before s in Equation ( 14).The workload and the waiting time have the same HT limit with the diffusionscale scaling in Equation ( 14) but different limits with the fluid-scale scaling in Equation ( 13).To capture the clear differences shown in Figure 1, we obviously want the stronger fluid scaling in Equation ( 13).The HT functional central limit theorem (FCLT) in Theorem 3 for the scaling in Equation ( 13) in Section 6 helps interpret Figure 1.
It is important to note that if we have a constant service rate with this scaling, then the waiting times explode as ρ ↑ 1 because the instantaneous traffic intensity ρ(t) ≡ λ(t) > 1 over intervals growing as ρ ↑ 1; this case is analyzed in Choudhury et al. (1997).
We also establish an HT functional weak law of large numbers (FWLLN) in Theorem 2, which yields a deterministic fluid approximation.However, it is not very useful because it shows that our proposed control with ξ 1 stabilizes the waiting time perfectly for all η as ρ ↑ 1 (but it helps to see that nothing bad happens).

The Model
In this section, we specify the general model, defining the arrival process in Section 4.1 and the basic queuing stochastic processes in Section 4.2.We specialize to the periodic G t /G t /1 model in Section 4.3.We show that the workload is stabilized by the rate-matching control in Equation ( 1), extending the results for the queuelength process in Whitt (2015).

The Arrival Process
We represent the periodic arrival counting process A as a deterministic time transformation of an underlying rate-1 counting process N with associated sequence of interarrival times {U k : k ≥ 1} by where λ is the arrival-rate function.This is a common representation when N is a rate-1 Poisson process; then A is a "nonhomogeneous Poisson process" (NHPP).For the G t /G t /1 model, N is understood to be a rate-1 stationary point process.Hence, for the GI t /GI t /1 model, N is an equilibrium renewal process with time between renewals having mean 1, for which the first interrenewal time U 1 has the equilibrium distribution.
The representation in Equation ( 15) has been used frequently for processes N more general than NHPPs, an early source being by Massey and Whitt (1994).
For the sinusoidal arrival-rate function in Equation ( 5), the associated cumulative arrival-rate function is We only consider the case ρ < 1, under which a proper steady state exists under regularity conditions (which we do not discuss here).Behavior differs for short cycles and long cycles.For the case of a constant service rate, there are two important cases for the relative amplitude: (1) 0 In the first case, we have ρ(t) < 1 for all t, where ρ(t) ≡ λ(t) is the instantaneous traffic intensity, but in the second case, we have intervals with ρ(t) ≥ 1, where significant congestion can build up.If there is a long cycle as well, the system may be better understood from fluid and diffusion limits, as in Choudhury et al. (1997).However, this difficulty can be avoided by a service-rate control.

The General
We consider a modification of the standard single-server queue with unlimited waiting space where customers are served in order of arrival.Let {V k } be the sequence of service requirements.As in Whitt (2015), we separately define the rate at which service is performed from the service requirement.Given the arrival counting process A(t) defined in Section 4.1, let the total input of work over the interval [0, t] be the random sum Let service be performed at time t at rate μ(t) whenever there is work to perform.Paralleling the cumulative arrival rate Λ(t) defined in Equation ( 15), let the cumulative available service rate be Let the net-input process of work be X(t) ≡ Y(t) − M(t), t ≥ 0. Then we can apply the reflection map to the net input process X(t) to represent the workload (the remaining work in service time) at time t, starting empty at time 0, as In this setting it is elementary that the continuous-time (virtual) waiting time (before starting service) at time t, which we denote by W(t), can be related to L(t).
Lemma 1 (Waiting Time Representation).The waiting time at time t can be represented as Proof.By definition, for M t (u) above, as claimed in Equation ( 19).□  2).For that purpose, we exploit the arrival-process construction in Equation ( 15) in terms of the stationary processes N ≡ {N(t) : t ≥ 0} and V ≡ {V k : k ≥ 1} in Equation ( 15).Let the associated service-rate function μ(t) also be periodic with cycle length c, average service rate μ 1, and bounds 0 As in Ma and Whitt (2018) and earlier in Loynes (1962) and chapter 6 in Sigman (1995), we now convert the standard representation of the workload process in Section 4 to a simple supremum by using a reverse-time construction.To do so, we extend the stationary processes {N(t)} and {V k } to the entire real line.We regard the periodic arrival rate and service rate as defined on the entire real line as well, with the functions fixed by their position within the periodic cycle at time 0. With those conditions, the reverse-time construction is achieved by letting the interarrival times and service times be ordered in reverse time going backward from time 0. Then Ã(t) counts the number of arrivals in To exploit the reverse-time representation, let be the reverse-time cumulative arrival-rate function starting at time yc within the periodic cycle [0, c], 0 ≤ y < 1, and Λ−1 y is its inverse function, which is well defined because Λy (t) is continuous and strictly increasing.As an analogue of Equation ( 21) for the cumulative service rate, let We let the service requirements V k come from a general stationary sequence with With this reverse-time representation, the workload at time yc in the system starting empty at time yc − t can be represented as where Xy is the reverse-time net input of work starting at time yc within the cycle of length c.The other quantities in Equation ( 23) are the reverse-time cumulative arrival-rate function Λy (t) in Equation ( 21) with inverse Λ−1 y (t) and the reverse-time cumulative service-rate function My in Equation ( 22) with inverse M−1 y .The equality in distribution in Equation ( 23) holds because N is a stationary point process, which is a point process with stationary increments and a constant rate.
As t → ∞, L y (t) ↑ L y (∞) ≡ L y with probability one (w.p.1) as t → ∞, for Even though Equation ( 23) is valid for all t, we think of the system starting empty at times −kc, for k ≥ 1, so we let yc − t −kc or, equivalently, we stipulate that t c(k + y), 0 ≤ y < c, and consider successive values of k and let k → ∞ to get Equation ( 24).That makes Equation ( 23) valid to describe the distribution of L(c(k + y)) for all k ≥ 1.We now observe that the time transformation in Equation ( 23) shows that the periodic G t /G t /1 model is actually equivalent to a G/G t /1 model with a stationary arrival process and a new cumulative service-rate function My ( Λ−1 y (t)).Corollary 1 (Conversion of G t /G t /1 to an Equivalent G t /G/1).In addition to representing the periodic steady-state workload L y in a periodic G t /G t /1 model as a periodic steady-state workload in a periodic G/G t /1 model, which has a stationary stochastic input and a deterministic service rate, as shown in Equation (24), we can represent it as a periodic steady-state workload in a periodic G t /G/1 model, which has a periodic stochastic input and a constant service rate, via Corollary 2 (The Associated Periodic Steady-State Waiting Time).The periodic steady-state waiting time associated with the periodic steady-state workload in Equation ( 24) is Proof.Apply the reasoning of Lemma 1. □ In Whitt (2015), we showed that the rate-matching service-rate control in Equation ( 1) stabilizes the queuelength process.Now we establish the corresponding result for the workload.
Theorem 1 (Stabilizing the Periodic Workload).If the rate-matching control in Equation (1) is used, then L y d L for L y in Equation ( 24), where L is the steady-state workload in the associated (stable) stationary G/G/1 model; that is, which is independent of y.
Proof.With the rate-matching control, we have M(t) cΛ(t) and My (t) c Λy (t), t ≥ 0. As a consequence,

The Simulation Search Algorithm
The rare-event simulation algorithm from Ma and Whitt (2018) exploits the classic rare-event simulation algorithm for the GI/GI/1 queue, exploiting importance sampling using an exponential change of measure, as in chapter XIII of Asmussen (2003) and chapter VI of Asmussen and Glynn (2007).Hence, our simulation algorithm applies to the GI t /GI t /1 queue.It was shown in Ma and Whitt (2018) that the algorithm is effective for estimating the mean as well as small tail probabilities (also see Ma and Whitt 2016).
5.1.The GI t /GI t /1 Model In the GI t /GI t /1 setting, the underlying rate-1 process N is an equilibrium-renewal process, which means that U 1 has the stationary-excess or equilibrium distribution U e , which may be different from the i.i.d.distributions of U k , k ≥ 2. Also in the GI t /GI t /1 setting, the service times V k are i.i.d. with distribution V and are independent of the arrival process.The simulation algorithm exploits the discrete-time representation of the workload L y in Equation ( 24) and the waiting time W y ; that is, where M y is the same as M t , which is the forward integral of the service rate starting from position y within a cycle.We exploit the rare-event simulation algorithm in Ma and Whitt (2018), which is based on an exponential change of measure; we refer to Ma and Whitt (2018) for background.In that setting, we use the underlying measure P θ * determined for GI/GI/1 queue.We again use the same notations X k (ρ) V k − ρ −1 U k and partial sum process S n ≡ n k 1 X k for GI/GI/1 and define the new associated process which is the process inside the supremum function.To avoid duplication of notation, we let the likelihood function here be denoted by Ψ instead of L. Then the estimator of the rare-event probability for W y can be derived as follows: where (θ * ) is the exponentially tilted likelihood ratio for process Q n at n M y (b), and m X (θ * ) is the moment-generating function of X at θ * .As in Ma and Whitt (2018), the first 5.2.The Extended Algorithm for the GI t /GI t /1 Model Here is a summary of the extended algorithm to estimate the tail probabilities in the GI t /GI t /1 queue with average service rate 1 and average arrival rate ρ: 1. Construct a table of the inverse cumulative arrival-rate function ρ Λ−1 y (same as for GI t /GI/1).2. Determine the required length of partial sums (n s ) needed in each application (same as for GI t /GI/1).3.For each replication, generate the required vectors of exponentially tilted interarrival times ρ −1 Ũ and service times Ṽ from F −θ * ρ −1 U and F θ * V , respectively (same as for GI t /GI/1).4. Calculate the associated vectors of S n and Q n , and find out the stopping time τ Q M y (b) , which is the hitting time of Q n at level M y (b).This step is different from that for GI t /GI/1 in that first we need to calculate M y (b) as the hitting level instead of b, and second we calculate vector Q n different from R n in an additional function My in the second term.
5. Use the preceding estimator to calculate the tail probability P(W y > b) for each replication (same as for GI t /GI/1).
6. Run N i.i.d.replications, and calculate the mean of the estimated values of P(W y > b) (same as for GI t /GI/1).

Explicit Representations for the Sinusoidal Case
Here we summarize the expressions for all the basic deterministic rate functions in our sinusoidal examples, extending Equations ( 5), (6), and (16): (31)

The Search Algorithm
We use an elementary iterative search algorithm, fixing an initial value of η at the mean for the steady-state model, ρ/(1 − ρ), and searching first over ξ and then over each variable until we get negligible improvement.This simple approach is substantiated by estimating the structure of the objective function.Figure 2 shows that the function is not convex as a function of η but suggests that it is unimodal with a unique global minimum, supporting our simple procedure.The plots for the maximum deviation max 0≤y≤c {E[W y ]} − max 0≤y≤c {E[W y ]} are similar.We perform the search with fewer points y and replications in the initial stages and then confirm with more points, 40 values of y and 40, 000 replications, which yield excellent statistical precision, as can be seen from the narrow confidence interval bands in Figure 1.

Supporting Heavy-Traffic Limits for Periodic Queues
In this section, we obtain an HT FWLLN and an HT FCLT for the periodic G t /G t /1 model with a general servicerate control of the form in Equation ( 6).The HT FCLT produces a limit depending on an asymptotic time lag η and damping factor ξ, which arise from HT limits; see condition (58) in Theorem 3 and the conclusion in Equation ( 50).Thus, we reduce the optimization problems over the parameter pairs (η ρ , ξ ρ ) in Equations ( 9) and ( 10), asymptotically as ρ ↑ 1, to diffusion-control problems with the parameter pairs ( η, ξ).

The Underlying Rate-1 Processes
As in much of the HT literature, we start by introducing basic rate-1 stochastic processes, but here we consider service requirements instead of service times.We assume that the rate-1 arrival and service-requirements processes N and V specified in Section 4 are independent and each satisfies an FCLT.To state the result, let Na n and Ŝv n be the scaled processes defined by with ≡ denoting equality in distribution and x denoting the greatest integer less than or equal to x.We assume that where $ is the usual function space of right-continuous real-valued functions on [0, ∞) with left limits, and ⇒ denotes convergence in distribution, as in Whitt (2002), whereas B a and B s are independent standard (mean 0, variance 1) Brownian motion processes (BMs).The assumed independence implies joint convergence in Equation ( 33) by theorem 11.4.4 of Whitt (2002).We emphasize that GI assumptions are not needed, but this is an important special case.If the service times V k are i.i.d.mean-1 random variables with variance, also the squared coefficient of variation (scv), c 2 s , then the limit in Equation (33) holds with service variability parameter c s .Similarly, if the base arrival process is a renewal process or an equilibrium renewal process with times between renewals having mean 1 and variance (and scv) c 2 a , then the limit in Equation ( 33) holds with arrival variability parameter c a .(See Nieuwenhuis 1989 for theoretical support in the case of an equilibrium renewal process.) For the queuing HT FCLT, we will apply theorem 9.3.4 of Whitt (2002), which refers to the conditions of theorem 9.3.3.Those conditions require a joint FCLT for the partial sums of the arrival and service processes, notably (3.9) on p. 295.This convergence follows from the FCLTs we assumed for Na n and Ŝv n in Equation ( 33).In particular, the assumed FCLT for N a n implies the associated FCLT for the partial sums of the interarrival times by theorem 7.3.2and corollary 7.3.1 of Whitt (2002).

A Family of Models
As a basis for the HT FCLT, we create a model for each ρ, 0 < ρ < 1.We do so by defining the arrival-rate and service-rate functions.
6.2.1.The Arrival-Rate and Service-Rate Functions.Let the arrival-rate function in model ρ be as in Equation ( 13) in the setting of Equations ( 2)-( 4).As a further regularity condition, we also require that the function s be an element of the function space $, as in Whitt (2002).Then the associated cumulative arrival-rate function in model ρ is where for s again being the periodic function in Equations ( 2)-( 4).From Equations ( 34) and ( 35), we see that the associated arrival-rate function obtained by differentiation in Equation ( 34) is indeed λ ρ (t) in Equation ( 33).
The time scaling in Equations ( 33) and ( 34) implies that the period in model ρ with arrival-rate function λ ρ (t) in Equation ( 33) is c ρ c(1 − ρ) −2 , where c is the period of s(t) in Equations ( 2)-( 4).Thus, the period c ρ in model ρ is growing with ρ.This scaling follows lemma 5.1 and theorem 5.2 of Whitt (2015), with n there replaced by (1 − ρ) −2 .In particular, the scaling here is in fluid or FWLLN scale and thus is different from the diffusion or FCLT scaling in theorem 3.2 of Whitt (2014) and theorem 2 of Ma and Whitt (2018).
Let A ρ (t) ≡ N a (Λ ρ (t)) be the arrival process in model ρ, which is obtained by using the cumulative arrivalrate function Λ ρ in Equation ( 34) in place of Λ in Equation ( 15).Given that definition, we see that the cumulative arrival rate is indeed We now define associated scaled time-varying service-rate functions.These are the rate-matching service-rate functions in Whitt (2015) modified by a time lag and a damping factor.In particular, where s is the periodic function with period c in Equation (3), whereas η ρ is the ρ-dependent time lag and ξ ρ is the ρ-dependent damping factor.From Equations ( 37) and (3), we see that the average service rate is μρ 1 for all ρ.As a consequence, the average traffic intensity is λρ / μρ ρ for all ρ, whereas the instantaneous traffic intensity at time t is λ ρ (t)/μ ρ (t), t ≥ 0, which is a more complicated periodic function, again with period c.
6.2.2.The Associated Queuing Processes.Having defined the family of arrival processes A ρ (t) and deterministic service-rate functions M ρ (t) above, we define the other queuing processes Y ρ (t), X ρ (t), L ρ (t), and W ρ (t) as in Section 4.2.Let the completed-work process be defined by We now can apply Lemma 1 in Section 4 to express the waiting-time process as The (virtual) waiting time W ρ (t) represents the time that a hypothetical arrival at time t would have to wait before starting service.As in equations (3.7) and (3.8) of Whitt (2015), we can define the queue-length process (number in system) and the departure process in model ρ jointly.We can also express the departure process in terms of the workload process instead of the queue-length process by but we do not focus on the departure and queue-length processes here.
Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control

The Scaled Queuing Processes
We start with the FWLLN-scaled processes.First, let the scaled deterministic rate functions be for Λ ρ (t) and M ρ (t) in Equation ( 37).We immediately see that where for S(t) in Equation ( 35).
Let the FWLLN-scaled arrival stochastic process be defined by Let the input, net-input, workload, completed-work, and waiting-time components of the FWLLN-scaled the vector ( Āρ , Ȳρ , Xρ , Lρ , Cρ , Wρ ) be defined in the same way.
Then let the associated FCLT-scaled deterministic rate functions be defined by for Λ f in Equation ( 43).Let the associated FCLT-scaled stochastic processes be defined by

The HT FWLLN
We start with the HT FWLLN.The limit provides a deterministic fluid approximation.However, simple fluid approximations evidently are too crude to provide much help.Corollary 3 shows that rate-matching control stabilizes both the workload and the waiting time for the fluid approximation.Let $ k be the k-fold product space of $ with itself, let ⇒ denote convergence in distribution, and let x • y be the composition function defined by (x • y)(t) ≡ x(y(t)).Let a ∧ b ≡ min {a, b}, and let ψ : D → D be the standard one-dimensional reflection map as in section 13.5 of Whitt (2002); that is, Theorem 2 (HT FWLLN).Under the definitions and assumptions in Section 6 above, if ξ ρ → ξ and η ρ → η as ρ ↑ 1, and the system starts empty at time 0, then and for ( Āρ , Ȳρ , Xρ , Lρ , Cρ , Wρ ) defined in Equation (44), where for Λ f (t) in Equation (43) with S(t) in Equation (35), M f (t) in Equation (48), and ψ being the reflection map in Equation (47).
Proof.We successively apply the continuous mapping theorem (CMT) using the functions in sections 12.7 and 13.2-13.6 of Whitt (2002).First, observe that Equation ( 48) is a minor modification of Equation ( 41).Let Na ρ and Sρ denote Na n and Sv n , respectively, where, paralleling Equation (32), we let Na n (t) ≡ n −1 N a (nt) and Sv n ≡ n −1 S v nt , t ≥ 0, and then let n (1 − ρ) −2 .Then observe that Āρ Na ρ • Λρ and Ȳρ Sρ • Āρ , so we can apply the CMT with the composition map.The limit for Xρ follows from the CMT with addition, and then the limit for Lρ follows from the CMT with the reflection map in Equation (47).To establish the limit for the scaled waiting time Wρ (t) in $, we apply the CMT with the inverse function.Finally, the limit for Cρ again follows from the CMT with addition.□ We obtain stronger results in special cases.
Remark 4 (Stabilization Achieved by Many Fluid Models).It is evident that the conclusion of Corollary 3 holds for any single-server fluid model with arrival rate λ(t) and service rate μ(t), provided that μ(t) ≥ λ(t) for all t.The (η, ξ) controls are intended to address the time-varying arrival rate in the more general stochastic setting.
As a modification of Corollary 3, we can have all customers wait exactly η if we provide no service until time η.
Corollary 4 (Stabilizing the Waiting Time at Any Positive Value).In addition to the conditions of Theorem 2, if ξ 1 and M f (t) 0, 0 ≤ t < η, then M f (t) Λ f (t − η), t ≥ η, for a fixed time lag η > 0, so that and Corollary 4 (Sinusoidal with Damped Time Lag).In addition to the conditions of Theorem 2, suppose that for positive constants β and γ with β < 1, so that s(t) is periodic with period c ≡ c γ 2π/γ.Then For the special case ξ 1, W(t) η.If in addition, and η < π/γ, the supremum in ( 55) is attained at s * (π/2γ) − η/2 , so that Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control for t ≥ c + η.As η ↓ 0, Remark 5 (The Impact of High or Low Frequency).Corollary 5 shows the impact of high or low frequency.First, it is well known that high frequency has negligible impact, because performance tends to be determined by the behavior of the cumulative arrival rate function Λ(t) in ( 15) rather than the rate function λ(t).From ( 54) and ( 55), we see that S(t) → 0 and L(t) → 0 as γ → ∞.On the other hand, for any fixed t, s(t) → 0 as γ → 0.

The HT FCLT
We now state our main HT result: the HT FCLT with periodicity in fluid scale, as in (13).We present the proof in Section 7 after discussing consequences here.
We now draw attention to some important consequences.First, Theorem 3 establishes an HT time-varying (TV) Little's law (LL), paralleling the many-server heavy-traffic (MSHT) TV LL in Sun and Whitt (2018) and exposed for the rate-matching control in theorem 5.2 of Whitt (2015).This is a time-varying version of the familiar state-space collapse, which goes back to the early HT papers (e.g., Whitt 1971).We remark that the relation is different from the time-varying LL discussed in Bertsimas and Mourtzinou (1997), Fralix and Riano (2010), Sigman and Whitt (2019), and Whitt and Zhang (2019).
Corollary 6 (HT Time-Varying Little's Law).Under the conditions of Theorem 3, the limit processes are related by We now consider an alternative deterministic limit to the HT FWLLN in Theorem 2. Now we assume that the FCLT holds with the variability parameter set equal to 0. For this purpose, we assume that s(t) is differentiable and let ṡ(t) be its derivative.
Corollary 7 (The Case of No Variability).If c x 0 and s(t) is differentiable in addition to the conditions of Theorem 3, then In the sinusoidal case with s(t) ≡ β sin γt in (5), For β 1 and γ → 0, which is strictly positive over subintervals if ξ > 0.
For the nondegenerate sinusoidal arrival rate function, the derivative in (64) of Corollary 7 implies it is not always possible to stabilize the limiting time-varying diffusion process Ŵ with ξ > 0 in Theorem 3. We conjecture that it is never possible to stabilize it perfectly.
We now establish conditions for the optimality of an ( η * , ξ * ) control for the limiting diffusion control problem for either formulation ( 9) or (10).Our proof will exploit uniform integrability (UI); see p. 31 of Billingsley (1999).
Corollary 8 (Optimality for the Limiting Diffusion Process).Consider the special case of the GI t /GI t /1 model with is the optimal control for problem (9) or (10), then the limiting control ( η * , ξ * ) is optimal for the corresponding diffusion control problem.
Proof.We let ( η, ξ) be any alternative control for the limiting diffusion process.Then let ( ηρ , ξρ ) be an associated control for model ρ, 0 < ρ < 1, where ηρ ≡ η/(1 − ρ) and ηρ ≡ 1 + (1 − ρ) ξ.Then, by this construction, condition (58) holds for the family ( ηρ , ξρ ).We next want to show that the convergence in distribution can be extended to convergence of the means for all t, which requires uniform integrability uniformly in t; see p. 31 of Billingsley (1999).We use the bounds on the second moments to show that it holds.
Toward that end, we exploit the upper bounds for the workload process in the G t /G t /1 model in terms of the associated workload process in the stationary G/G/1 model from section 3 of Ma and Whitt (2018).These bounds extend directly to the G t /G t /1 model by virtue of Corollary 1.These bounds show that the mean workload is bounded above uniformly in y over the interval [0, c].These bounds also apply to the waiting time process because W(t) ≤ L(t)/μ L , where μ L > 0 is a lower bound on the service rate, which follows from ( 4) and ( 6).For the stationary GI/GI/1 model, finite second moments imply the existence of the first moments of the waiting time and uniform integrability needed for convergence; see p. 31 of Billingsley (1999) and sections X.2 and X.7 of Asmussen (2003).
Finally, we observe that our optimal policy (η * ρ , ξ * ρ ) has expected value greater than or equal to the alternative policy ( ηρ , ξρ ) for all ρ, whereas both converge as ρ → 1.Hence, the limit of the optimal policies, ( η * , ξ * ) must be at least as good as ( η, ξ).□ We apply Corollary 8 to support our numerical calculations by observing that (η * ρ , ξ * ρ ) when scaled as in (58) converges to a limit.We thus deduce that the limit must be the optimal policy for the diffusion.However, this numerical evidence is not a mathematical proof.Moreover, although the numerical evidence is good, it is not exceptionally good, especially for ξ * ρ , as can be seen from Table 1 in Section 9.1.

Proof of Theorem 3
To establish (59), apply ( 37) and ( 45) to obtain Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control where on the third line we have subtracted and added the term ξ ρ S(t) and on the last line we have differentiated using by assumption (58).We used the assumed continuity of s to have S be continuously differentiable, so the derivative of S(t) holds uniformly in t over bounded intervals.
We next establish Equation (60).First, the limit for Âρ is given in lemma 5.1 of Whitt ( 2015), but we need to make an adjustment because the arrival rate in model ρ is chosen to be ρ here as opposed to 1 before.From Equations ( 43), ( 42), and (45), we see that for all ρ, where Λ f (t) is defined in Equation ( 43).Then the limit for Âρ follows from the standard argument for random sums.The key is to observe that where Nρ is defined to be Nn in Equation (32) for n (1 − ρ) −2 .So we can start with the joint convergence We then apply convergence preservation with the map g(x, y, z) x • y + z (composition plus addition) as in section 13.3 of Whitt (2002) A variant of the random-sum argument holds for Ŷρ too.In particular, we start with the joint convergence The joint convergence holds by virtue of theorems 11.4.4 and 11.4.5 of Whitt (2002).We then apply convergence preservation with the map g(x, y, z) x • y + z (composition plus addition) as in section 13.3 of Whitt (2002) to get Then the limits for Xρ and Lρ follow from the continuous mapping theorem with standard reflection map reasoning, for example, as in chapter 9 of Whitt (2002), even though the service-rate function is now more general.However, the waiting time requires a new treatment.The limit follows from the definition of the scaled servicerate control in Equation ( 37) and the first-passage-time representation of the waiting time in Equation ( 39).The structure and result are similar to the Puhalskii (1994) theorem and related results in section 13.7 of Whitt (2002), but they evidently do not apply directly.
We apply Taylor's theorem to a perturbation of S in Equation ( 35).The essential idea is that uniformly in t and u over bounded intervals.Just as in Equation ( 67), we use the assumed continuity of s to have S be continuously differentiable, so the derivative of S(t) holds uniformly in t and u over bounded intervals.
For the specific application, let and By combining Equation ( 74) and the two limits in condition (58), we see that ζ ρ (t, u) is asymptotically negligible as ρ → 1 uniformly in t and u over bounded intervals.We will use this at the critical final step in the following representation.
For the final step, to simplify, we make the entire argument deterministic by using the Skorohod representation theorem, as in theorem 3.2.2 of Whitt (2002), to replace the stochastic convergence Lρ ⇒ L in $ by associated convergence w.p.1.Then we see from line 6 of Equation ( 77) that in the infinimum it suffices to consider u only just beyond L(t)/λ f (t), which for t in a bounded interval is bounded for each sample path because λ f (t) has been assumed to be bounded below, whereas L(t) is bounded above, for t in a bounded interval.Thus, we can write for t and u over specified bounded intervals, K an appropriate positive constant, and for an appropriate ū.Given that Lρ ⇒ L and ζ ↑ ρ → 0 in $, we can use the standard sandwiching argument (uniformly over bounded time intervals) to obtain convergence Ŵρ (t) ⇒ L(t)/λ f (t) ≡ Ŵ(t) in $, which completes the proof.□

An HT FCLT with Periodicity in the Weaker Diffusion Scaling
In this section, we establish an HT FCLT with periodicity holding in the weaker diffusion scale instead of in the fluid scale, as was done in Section 6.The scaling here follows Whitt (2014) and Ma and Whitt (2018) instead of Whitt (2015).In this scaling, the HT limit of the waiting time coincides with the HT limit for the workload process and so does not capture the differences we see in the simulations in preceding sections.

An Alternative Family of Models
We start with the same basic rate-1 processes in Section 6.1.We then create a model for each ρ, 0 < ρ < 1, now using Equation ( 14) instead of Equation ( 13).This yields the family of cumulative arrival-rate functions for S in Equation ( 35).Differentiating in Equation ( 79) yields the arrival-rate function in Equation ( 14).Just as before, the time scaling in Equations ( 14) and ( 79) implies that the period in model ρ with arrival-rate function λ ρ (t) in Equation ( 14) is c ρ c(1 − ρ) −2 , where c is the period of s in Equations ( 2)-( 4).Thus, the period c ρ in model ρ is growing with ρ.
Proof.We will be brief because most of the argument is essentially the same as in Whitt (2014) and Ma and Whitt (2018).First, the limit for Âρ is given in theorem 3.2 of Whitt (2014).Then the limit for Ŷρ follows from theorem 9.3.4 of Whitt (2002), as noted in the proof of theorem 2 in Ma and Whitt (2018).(See C(t) in equation 9.2.4 and C n in equation 9.3.4 and theorem 9.3.4 of Whitt 2002.)Then the limits for Xρ and Lρ follow from the standard reflection mapping argument even though the service-rate function is now more general.Again, the waiting time requires a new treatment.The limit follows from the first-passage-time representation in Equation (39).In particular, paralleling Equation ( 77), letting Mρ (t, u) ≡ M ρ ((1 − ρ) −2 t + u), we have for t ≥ 0, where which is asympotically negligible as ρ → 1 uniformly in compact intervals, given the conditions on η ρ and ξ ρ .
As technical support for the last step, note that for s U in Equation ( 4).Also add and subtract ξ ρ S(t) and treat the two terms separately; that is, Hence, we can apply the continuous mapping theorem for the inverse in section 13.6 of Whitt (2002) to get Ŵρ ⇒ L in $ as ρ → 1, jointly with the other limits.□

Simulation Examples in the Setting of Section 6
In this section, we report results of simulation experiments to evaluate the new optimal (η * ρ , ξ * ρ ) controls as a function of ρ for models scaled according to Theorem 3, specifically by Equations ( 13), (34), and (37), so that we can see the systematic behavior.9.1.1.Sinusoidal Examples.We start with sinusoidal examples and then consider nonsinusoidal examples.Table 1 shows results for four values of the traffic intensity ρ with ρ ↑ 1 for the sinusoidal model in Equations ( 2)-( 6) with HT scaling in Equation ( 13) with parameters (ρ, β ρ , γ ρ ) (ρ, 0.2, 2.5(1 − ρ) 2 ).For this case, we found that the solutions to optimization problems ( 9) and ( 10) are identical, to within our statistical precision.Hence, our solutions are for both problems.
Table 1 shows the estimated optimal controls η * ρ and ξ * ρ in each case plus scaled versions consistent with condition (58).Table 1 shows that the relative error is roughly independent of ρ, being less than 1% in each case.Table 1 also shows that the limit η * ≈ 1.45 is rapidly approached by (1 − ρ)η * ρ /ρ, whereas the limit ξ * ≈ 1.8 is roughly approached by (ξ * ρ − 1)/(1 − ρ), both of which are consistent with condition (58).The results support Theorem 3, but unfortunately, the rate of convergence in the control parameters is not fast.Evidently the optimal damping control ξ * ρ is more problematic.For the model in Table 1, Figure 3 shows the expected periodic steady-state virtual waiting time (solid blue line), the expected steady-state workload (the dashed red line), and arrival rate multiplied by the mean waiting time (the dotted green line) for ρ 0.8 (left) and ρ 0.95 (right).As in Figure 1, the 95% confidence interval bands are included, but they can only be seen by zooming in.
We also considered alternative values of the relative amplitude β.Table 2 shows the solutions to the minimum-deviation optimization problem in Equation ( 10) for the sinusoidal model in Table 1 except that β has been increased to β 0.8 from 0.2.Table 2 shows that the relative error is roughly independent of ρ, but the relative error has increased to approximately 10% from approximately 1% in Table 1.Unlike in Figure 3, it Notes.The first row shows the expected periodic steady-state virtual waiting time (solid blue line), the expected steady-state workload (the dashed red line), and arrival rate multiplied by the mean waiting time (the dotted green line) for ρ 0.8 (left) and ρ 0.95 (right) in the base case (β, γ) (0.2, 2.5).The second row shows arrival rate divided by ρ and service rate for ρ 0.8 (left) and ρ 0.95 in the same base case.The optimal control parameters are (η * ρ , ξ * ρ ) (5. is evident that the (η * ρ , ξ * ρ ) control does not stabilize the expected waiting time perfectly, either for fixed ρ or asymptotically as ρ → 1.
Unlike the rate-matching control in Whitt (2015), which stabilizes the entire queue-length distribution, the optimal modified (η, ξ) control neither stabilizes the mean perfectly nor does it stabilize the entire waiting-time distribution.However, it seems to do a reasonable job of both.Figures 4, 5, and 6 illustrate by showing plots of the time-varying (i) standard deviation SD[W y ], (ii) delay probability P(W y > 0), and (iii) full complementary cumulative distribution function (CCDF) {P(W y > x) : x ≥ 0} for the two cases in Figure 3, that is, for ρ 0.8 (left) and ρ 0.95 (right).
Very roughly, Figures 3, 4, 5, and 6 are consistent with the time-varying waiting-time distribution being exponential, as in the M/M/1 stationary model.We should not be surprised that the results look similar for ρ 0.8 and ρ 0.95 because they are scaled to be part of the family of systems satisfying the HT limit.To consider  7 shows plots of the estimated waiting-time CCDF (left) and probability density function (PDF, right) for ρ 0.95 in Figure 1, which is not scaled.Without the scaling, the cycles are relatively short.Figure 7, especially the PDF, shows much more varied behavior in this difficult short-cycle setting.
We now show two candidate modifications of the control used in Figure 1.First, Figure 8 shows the analogue of Figure 1 where we fix ξ 1 and use only the single control parameter η.As stated in Remark 2 in Section 2, if we let ξ 1 and optimize over η, then for ρ 0.8 we get η * ρ 5.93 and a maximum deviation of 0.4109, which yields approximately 10% relative error instead of less than 1%.For ρ 0.95, η * ρ 28.3, the maximum deviation is 3.034, and the relative error is approximately 14%.
On the right appears the diffusion limit for the net input X(t) −t − M(t) when c x 0, showing that condition (65) holds.into the difficult cases with zero or near-zero λ f (t) can be gained from the time-varying Little's law in Corollary 6 and steps ( 77) and ( 78) in the proof of Theorem 3.

Simulation Examples for the Alternative Scaling in Section 8
We now consider four simulation examples in the alternative HT scaling in Section 8.This is the same HT scaling as in Ma and Whitt (2018).We consider the base case of β 1, γ 2.5, and use Specifically, we consider cases with ρ 0.84, 0.92, 0.96, 0.98.Here we use the lags η ρ 5.25, 11.5, 24, 49 calculated by ρ/(1 − ρ), the scaler ξ ρ ρ. (These are consistent with Theorem 4.) Figures 13 and 14 show the expected periodic steady-state waiting time (the solid blue line) and the expected steady-state workload (the dashed red line).Figures 13 and 14 show that stabilization is not achieved

Conclusions
In this paper, we extended the rare-event simulation algorithm for the periodic GI t /GI/1 model in Ma and Whitt (2018) to the periodic GI t /GI t /1 model and applied the new algorithm to study methods to stabilize the expected (virtual) waiting time over time.We studied the modification in Equation ( 6) of the rate-matching service-rate control in Equation (1a) to include a time lag η and a damping factor ξ. We developed and applied a simulation search algorithm to find optimal pairs of control parameters (η, ξ) for the control problems in Equations ( 9) and ( 10).Thus, we obtained a practical solution to the open problem in Whitt (2015) of developing an effective way to stabilize the expected waiting time in the periodic single-server model.We also established supporting HT limits for the general periodic G t /G t /1 model and showed that the control problems in Equations ( 9) and (10) converge to associated diffusion control parameters with  appropriate scaling.The scaling involves the conventional HT scaling associated with ρ ↑ 1, so time is scaled by (1 − ρ) −2 , whereas space is scaled by 1 − ρ, but in addition, to gain insight into the time-varying behavior, we identify and study three different scalings of the arrival-rate function.
As observed by Falin (1989), if the arrival-rate function is left unscaled, then the HT limit is the same as if the periodicity were not present at all.A major conclusion is that important insight into the time-varying performance can be gained by scaling the arrival-rate function as well.Moreover, as illustrated by Sections 6 and 8, there are two different natural scalings: first, there is the stronger scaling in the fluid scale in Section 6 as in Whitt (2015), and second, there is the weaker scaling in the diffusion scale in Section 8 as in Whitt (2014) and Ma and Whitt (2018).In the weaker scaling, the rate-matching control from Whitt (2015) stabilizes both the queue length and the waiting time, but in the stronger fluid scaling we see significant differences, consistent with the simulation results in Figure 1.This insightful scaling the fluid scale also yields a limiting diffusion control problem.
We conducted extensive simulation algorithms showing that the new (η, ξ) control is effective in stabilizing the expected waiting time.However, unlike the rate-matching control for the queue length process in  Whitt (2015), the new modified rate-matching control does not stabilize the expected waiting time perfectly, either for fixed ρ or in the HT limit.However, Figures 1 and 3 show that it stabilizes it remarkably well, whereas Figures 4, 5, and 6 show that it stabilizes the full waiting-time distribution quite well too.
It is interesting to consider the performance impact of time-varying arrivals.In Section 1, we observed that the difference between the stable average waiting time in Figure 1 and the value ρ/(1 − ρ) for the stationary model (4 on the left and 19 on the right) might be called the "average cost of periodicity," but Example 1 showed that the overall average expected waiting time with a service-rate control could be much less than in the stationary model.It remains to investigate more carefully.
Indeed, there remain many opportunities for future research, including the open problems mentioned in Section 1.5.It also remains to directly solve the diffusion-control problems with objectives ( 9) and ( 10) resulting from Theorem 4, and there are other methods worth carefully studying, such as modifications of the iterated staffing algorithm from Feldman et al. (2008) for single-server models.
Finally, we mention that the methods in this paper generalize and can be applied to other problems.First, the rare-event simulation algorithm in Section 5 applies to any GI t /GI t /1 model with other service-rate controls.Second, the heavy-traffic limits in Section 6 and 8 evidently extend to general G t /G t /1 models with other service-rate controls.More generally, simulation of converging stochastic processes is a promising way to numerically solve complex diffusion-control problems.

Figure 7 .
Figure 7. Plots of the Waiting Time CCDF (Left) and PDF (Right) for Four Places in the Cycle for the Example in Figure 1 with ρ 0.95 but Without Heavy-Traffic Scaling

Figure 14 .
Figure 14.The Expected Periodic Steady-State Virtual Waiting Time (Blue Line) and the Expected Steady-State Workload (Red Line) for the Alternative Scaling in Section 9.2 for ρ 0.96, β 0.04, γ 0.004, η ρ 24, ξ ρ 0.96, Yielding a Maximum Deviation 0.0228 (Left) and ρ 0.98, β 0.02, γ 0.001, η ρ 49, ξ ρ 0.98, Yielding a Maximum Deviation 0.0070 (Right) We have not yet solved this general optimization problem in Equation(7).Here are open problems, applying to the Markovian M t /M t /1 model and generalizations: 1 model.However, in general, E[W] is not a lower bound for the average of the periodic steady-state mean E[W y ] over a cycle (see Remark 2 and Example 1).
Figure 1.Estimates of the Periodic Steady-State Values of E[W y ] (Blue Solid Line), E[L y ] (Red Dashed Line), and λ(y)E[W y ] (Green Dotted Line) for the Optimal Control (η * , ξ * ) for the Sinusoidal Example with Parameter Triples (ρ, β, γ) (0.8, 0.2, 0.1) (Left) and (0.95, 0.2, 0.1) (Right), So That the Cycle Length Is c 2π/γ 62.8 Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control Stochastic Systems, Articles in Advance, pp.1-30, © 2019 The Author(s) 4.3.The Periodic G t /G t /1 Model As in Ma and Whitt (2018), we consider the periodic steady state of the periodic G t /G t /1 model with arrivalrate function in Equation ( Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control Stochastic Systems, Articles in Advance, pp.1-30, © 2019 The Author(s) Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control Stochastic Systems, Articles in Advance, pp.1-30, © 2019 The Author(s) Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control Stochastic Systems, Articles in Advance, pp.1-30, © 2019 The Author(s)

Table 1 .
The (Identical) Solutions to the Minimax and Minimum-Deviation Optimization Problems in (9) and (10) for the Sinusoidal Model in (2)-(6) with HT Scaling in (13) with Parameters (ρ, β ρ , γ ρ ) (ρ, 0.2, 2.5(1 − ρ) 2 ).Notes.The four cases are determined by ρ as specified in the first row.The mean waiting times are reported with and without space scaling in the bottom rows.Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control 80, 0.84) for ρ 0.8 and (27.7, 0.93) for ρ 0.95.The maximum minus minimum of EW y over a cycle equals 0.0321 for ρ 0.8 and 0.1425 for ρ 0.95.Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control Stochastic Systems, Articles in Advance, pp.1-30, © 2019 The Author(s)

Table 2 .
Solutions to the Minimum-Deviation Optimization Problem in (10) for the Sinusoidal Model in Table 1, Except β Has Been Increased to β 0.8 from 0.2 Notes.The four cases are determined by ρ as specified in the first row.The mean waiting times are reported with and without space scaling in the bottom rows.
Figure 4. Plots over One Cycle in Table 1: the Standard Deviation for ρ = 0.8 (Left) and 0.95 (Right) Notes.Estimates of the periodic steady-state standard deviation SD[W y ] for ρ 0.8 (left) and ρ 0.95 (right), shown in red solid line.Also displayed are the fluid arrival rate λ f 1 + s(t) (blue dashed line) and fluid service rate μ f 1 + ξ * ρ s(t − η * ρ ) (blue dotted line).Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control a very different case, Figure

Table 1 :
the Delay Probability for ρ = 0.8 (Left) and 0.95(Right)Notes.Estimates of the periodic steady-state probability of delay P(W y > 0) for ρ 0.8 (left) and ρ 0.95 (right), shown in red solid line.Also displayed are the fluid arrival rate λ f 1 + s(t) (blue dashed line) and fluid service rateμ f 1 + ξ * ρ s(t − η * ρ ) (blue dotted line).Plots of the Periodic Waiting Time CCDF for 4 places in the Cycle for the Example in Table 1: ρ = 0.8 (Left) and 0.95 (Right) Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control Stochastic Systems, Articles in Advance, pp.1-30, © 2019 The Author(s) Note.Estimates of the periodic steady-state ccdf P(W y > x) for four values of y: 0, 1/4, 1/2, 3/4 for ρ 0.8 (left) and ρ 0.95 (right).
Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control Stochastic Systems, Articles in Advance, pp.1-30, © 2019 The Author(s)

Table 3 .
Optimal η and ξ for ρ 0.8, β 0.5, Cycle Length c 60 Note.The four cases are determined by the value of p in (88) and the objective for the time-varying waiting time in (9) or (10) in the first two rows of the table.

Table 4 .
Optimal η and ξ for ρ 0.8, β 1, Cycle Length c 60 Ma and Whitt: Periodic Single-Server Queue with a Service-Rate Control Stochastic Systems, Articles in Advance, pp.1-30, © 2019 The Author(s) Note.The four cases are determined by the value of p in (88) and the objective for the time-varying waiting time in (9) or (10) in the first two rows of the table.