Transient error approximation in a Lévy queue

Abstract

Motivated by a capacity allocation problem within a finite planning period, we conduct a transient analysis of a single-server queue with Lévy input. From a cost minimization perspective, we investigate the error induced by using stationary congestion measures as opposed to time-dependent measures. Invoking recent results from fluctuation theory of Lévy processes, we derive a refined cost function, that accounts for transient effects. This leads to a corrected capacity allocation rule for the transient single-server queue. Extensive numerical experiments indicate that the cost reductions achieved by this correction can be significant.

1 Introduction

The issue of matching a service system’s capacity to stochastic demand induced by its clients arises in many practical settings. Typically, the resources available to satisfy demand are scarce and hence expensive. This forces the manager to consider a trade-off between the system efficiency and the quality of service perceived by its clients. In this paper, we focus on this trade-off in the context of the M / G / 1 queue, in which the variable amenable for optimization is the server speed $\mu$.

In general, optimizing the server speed $\mu$ in a single-server queue in a time-homogeneous environment, while trading off congestion levels against capacity allocation costs, does not pose any technical challenges. Typically, the objective function to be minimized, the total cost function, has the shape

$$\begin{aligned} \varPi _\infty (\mu ) = \mathbb {E}[Q_\mu (\infty )] + \alpha \mu = \frac{\lambda \mathbb {E}[B^2] }{2(\mu -\lambda \mathbb {E}[B])} + \alpha \mu , \end{aligned}$$

(1.1)

where $\mathbb {E}[Q_\mu (\infty )]$ denotes the expected steady-state amount of work given server speed $\mu$, and B describes the service requirement per arrival. The parameter $\alpha >0$ represents the relative capacity allocation costs incurred by deploying service rate $\mu$. This one-dimensional optimization problem yields the optimizer

$$\begin{aligned} \mu _\infty ^\star = \lambda \mathbb {E}[B] + \sqrt{\frac{\lambda \mathbb {E}[B^2]}{2\alpha }}. \end{aligned}$$

Despite the simplicity and tractability of the problem described above, the presence of the steady-state measure in the cost function in (1.1) should be handled carefully. By employing this particular cost structure, one automatically agrees with the underlying assumption of the system being sufficiently close to its steady state. However, referring to the practical applications of the single-server model, system parameters rarely remain constant over time. Moreover, planning periods for the optimization problem are naturally finite. Hence, the true expected costs incurred, which we denote by $\varPi _T(\mu )$, in addition depend on the length of the planning period T. Consequently, the usage of steady-state models for decision making needs to be justified by a more elaborate time-dependent or transient analysis for these type of settings.

The time-dependent behavior of the single-server queue received much attention in queueing theory. First efforts to analyze the time-dependent properties of the M / G / 1 queue date back to the 1950s and 1960s; for example, [7, 10, 17, 28, 29]. The analyses in these papers mostly yield implicit expressions for performance characteristics through Laplace transforms, integro-differential equations and infinite convolutions. More specifically, there is vast literature on the transient analysis of the M / M / 1 queue, with the goal to derive explicit expressions for queue length characteristics; see, for example, [3, 9, 23, 24]. These works provide a variety of explicit expressions for the transient dynamics, although the complexity of the resulting expressions, typically involving Bessel functions, exposes the intricate intractability of the matter. Consequently, approximation methods for insightful quantification of the dynamics based on numerical [20] or asymptotic methods have become prevalent in more recent literature. The asymptotic methods either exploit knowledge of the evolution of the queueing process as time t grows large [3, 21, 22], or as the arrival rate $\lambda$ is increased to infinity [1, 2, 11]. It is noteworthy that a substantial contribution to the transient literature is made by Abate and Whitt [1,2,3,4], who exploit the existence of a decomposition of the mean transient queue length and obtain expressions for the moments of the queue length and virtual waiting through probabilistic arguments in several queueing models. More recently, asymptotic methods have been used to justify the application of stationary performance measures in Markovian environments or to refine them; see, for example, [12, 30]. Other approximative methods known as uniform acceleration expansions [19] have been developed to reveal the asymptotic behavior of the single-server queue as a function of t, which are moreover able to capture time-varying arrival rates.

The majority of the works mentioned above do reflect on the error imposed by usage of steady-state performance metrics instead of the correct time-dependent counterpart. However, no light has been shed on the accumulation of this error over a finite period of time. To the best of our knowledge, the only work that addresses this issue is the paper by Steckley and Henderson [27], who compute an approximation for the error accumulated between the steady-state and transient delay probability. Our analysis on the other hand is centered around the mean workload, which requires a different approach. In addition, the focus in [27] is on performance measures only, while the main goal of our paper is to investigate the quality of staffing rules.

Although the M / G / 1 queue serves as the leading example in our analysis, we choose to use a more general framework for the arrival process of the queue. Namely, we let the server face a Lévy process. This gives the advantage that once we have obtained the results, we can apply them to broader queue input classes, such as Brownian motion and the Gamma process.

To shed light on the influence of the transience of the queueing process on traditional staffing questions, we will study the capacity allocation problem in the context of cost minimization in which the objective function is $\varPi _T(\mu )$, i.e., a function of both $\mu$ and T. We investigate how the invalidity of the stationary assumption is echoed through the operational cost accounting for congestion-related penalties.

Furthermore, we establish a result on the strict convexity of the function $\varPi _T(\mu )$, for almost all values of T (with a few minor exceptions for certain deterministic initial states), which is an essential property for convergence of both cost function and corresponding minimizer to their stationary counterparts.

As it will appear that an exact analysis of this disparity is intractable, we will present an explicit approximate correction to the conventional stationary objective function given by $\varPsi (\mu )/T$ and prove that

$$\begin{aligned} \varPi _T(\mu ) = \varPi _\infty (\mu ) + \frac{\varPsi (\mu )}{T} + O(1/T^2), \end{aligned}$$

with the help of recent results from the fluctuation theory of Lévy processes. Based on this refinement, we ultimately examine how incorporating transient effects changes the optimal capacity level and propose a refinement to the steady-state capacity allocation rule,

$$\begin{aligned} \mu _T^\star = \mu _\infty ^\star + \frac{\mu _\bullet }{T} + o(1/T). \end{aligned}$$

We moreover deduce an explicit expression for $\mu _\bullet$ in terms of the initial state and the first three moments of the service requirement per arrival. It is noteworthy that similar refined square-root staffing rules have been proposed for multi-server queues in the Halfin–Whitt regime; see, for example, [14,15,16, 25, 31]. In those cases, the relevant decision value is the number of servers and refinements are derived for $\lambda \rightarrow \infty$, whereas we consider the regime $T\rightarrow \infty$.

Building upon the insights gained through the analysis of this optimality gap, we reflect on the parameter settings of the underlying queueing process in which our refined capacity sizing rule yields significant improvement and in which cases it has little effect. Special emphasis is put on the relationship between the accuracy of the standard procedure and the length of the planning period.

The remainder of the paper is structured as follows. Section 2 is devoted to the model description and presents some preliminary results. The main result will be given in Sect. 3, and results regarding the optimization problem will be discussed in Sect. 4, followed by the validation of our novel techniques through numerical experiments in Sect. 5. We will give some concluding remarks and topics for further research in Sect. 6. We have deferred all proofs to the appendices.

2 Model description

2.1 A queueing model with Lévy input

The model that inspired our study is the standard M / G / 1 queue starting out of equilibrium. Customers arrive to the queue according to a Poisson process with rate $\lambda$, and each arrival has iid service requirement $B_i$, stemming from a common random variable B. Without loss of generality, we will assume $\mathbb {E}[B] = 1$ throughout. The server is able to remove $\mu$ amounts of work from the system per time unit; a variable we will refer to as the server speed. For example, if $\mu = 3$ and two customers are in the system with remaining service times 4 and 2, then the queue will be empty 2 time units later, provided that no new arrivals occur in the meantime. Let $N_\lambda (t)$ denote the number of arrivals until time t. Accordingly, the total work generated by the customers is given by

$$\begin{aligned} Z_\lambda (t) = \sum _{i=1}^{N_\lambda (t)} B_i. \end{aligned}$$

Furthermore, define $X_{\lambda ,\mu }(t) := Z_\lambda (t) - \mu t$. We call $X_{\lambda ,\mu }$ the net-input process. More generally, we assume throughout the paper that $X_{\lambda ,\mu }$ is a Lévy process. Specifically, we let $Z_\lambda$ be of the form $Z_\lambda (t) = U(\lambda t)$, where U is a spectrally positive Lévy process generated by the triplet $(a,\sigma ,\nu )$ and $\mathbb {E}[U(1)] = 1$. This restriction to spectrally positive processes is equivalent to stating $\nu (-\infty ,0)=0$ and is a vital assumption in our analysis. Subsequently, we assume the net-input process $X_{\lambda ,\mu }$ to be

$$\begin{aligned} X_{\lambda ,\mu }(t) = U(\lambda t) - \mu t, \qquad t \ge 0. \end{aligned}$$

(2.1)

Note that by setting $a=\sigma =0$ and $\nu = \lambda \, F_B$, where $F_B$ is the cumulative distribution function of B, we retrieve the original M / G / 1 queue. The stochastic process central to our analysis is the workload process $Q_{\lambda ,\mu }(t)$, $t\ge 0$, which describes the amount of work the server is facing at time t. The net-input process $X_{\lambda ,\mu }$ completely determines the trajectory of $Q_{\lambda ,\mu }$, namely

$$\begin{aligned} Q_{\lambda ,\mu }(t) = \max \left\{ Q(0) + X_{\lambda ,\mu }(t), \sup _{s\in [0,t]} [X_{\lambda ,\mu }(t)-X_{\lambda ,\mu }(s)]\right\} , \qquad t\ge 0, \end{aligned}$$

(2.2)

where Q(0) is the initial workload in the system. In fact, $Q_{\lambda ,\mu }$ is the reflected version of $X_{\lambda ,\mu }$ with reflection barrier at zero. Careful inspection of the structure also reveals that $X_{\lambda ,\mu }(t) \equiv X_{\lambda /\mu ,1}(\mu t) \equiv X_{1,\mu /\lambda }(\lambda t)$, so that

$$\begin{aligned} Q_{\lambda ,\mu }(t) \,{\buildrel d \over =}\,Q_{\lambda /\mu ,1}(\mu t) \,{\buildrel d \over =}\,Q_{1,\mu /\lambda }(\lambda t) \end{aligned}$$

(2.3)

for all $\lambda ,\mu ,t>0$. This identity will prove to be convenient for the numerical analysis in Sect. 5. For reasons of clarity, we omit the subscript $\lambda$ in our expressions if no ambiguity is possible.

The process $Q_{\mu }$ is a natural indicator of the level of congestion in the system and therefore a good choice for quantifying the Quality of Service (QoS) received by a client. We remark that alternative processes characterizing congestion in the system can be deduced directly from $Q_{\mu }(t)$. For example, consider the virtual waiting time process $V_{\mu }(t)$, which is the waiting time a customer would experience if he arrives at time t. This satisfies the relation $V_{\mu }(t) \equiv Q_{\mu }(t)/\mu$ for all $t\ge 0$. Likewise, the expected number of customers in the system $L_{\mu }(t)$ at time $t\ge 0$ is given by Little’s law:

$$\begin{aligned} \mathbb {E}[L_{\mu }(t)] = \lambda \, \mathbb {E}[V_{\mu }(t)] = \frac{\lambda }{\mu }\, \mathbb {E}[Q_{\mu }(t)]. \end{aligned}$$

To facilitate our investigation of the queueing model, we end this subsection by introducing some notation regarding the net-input and workload process and by stating a useful preliminary result concerning the stationary process $Q_{\mu }(\infty )$. Throughout the paper, we assume $\mu >\lambda$ to ensure ergodicity of the queue and convergence in distribution to the limit

$$\begin{aligned} Q_{\mu }(\infty ) := \lim _{t\rightarrow \infty } Q_{\mu }(t), \end{aligned}$$

for any initial state $Q(0)<\infty$. This random variable necessarily coincides with the stationary distribution of $Q_{\mu }(t)$. By $\kappa _U(\cdot )$ and $\kappa _{\mu }(\cdot )$, we denote the Lévy exponents of the processes U and $X_{\mu }$, respectively:

$$\begin{aligned} \kappa _{\mu }(\theta ) = \log \mathbb {E}[e^{\theta X_{\mu }(1)}] = \log \mathbb {E}[e^{\theta (U(\lambda ) - \mu )}] = \lambda \kappa _U(\theta ) - \mu \theta . \end{aligned}$$

Furthermore, define $u_k = \mathbb {E}[\{U(1) - \mathbb {E}U(1)\}^k]$ for $k=2,3,...$. Using this representation, we obtain the following preliminary result.

Lemma 1

Let $\mathbb {E}|U(1)|<\infty$, $u_2, u_3 < \infty$ and $\mu > \lambda$. If $Q_{\mu }(\infty )$ represents the steady-state distribution of the workload process, then

$$\begin{aligned} \mathbb {E}\left[ Q_{\mu }(\infty )\right] = \frac{\lambda u_2}{2(\mu -\lambda )},\qquad \mathbb {E}\left[ Q_{\mu }^2(\infty )\right] =\frac{\lambda ^2u_2^2}{2(\mu -\lambda )^2} + \frac{\lambda u_3}{3(\mu -\lambda )}. \end{aligned}$$

The proof of Lemma 1 follows directly by differentiation of the Laplace transform of $Q_\mu (\infty )$ and is therefore straightforward.

2.2 Finite horizon

For the purpose of this paper, we are interested in the dynamics of the workload process within a fixed time frame of length $T>0$. For all $0\le t \le T$, we assume that the parameters of the queue, $\lambda ,\mu ,u_2,u_3$, remain unchanged. If at $t=0$ the queue is not in steady-state corresponding to the specified parameters of the starting period, the process $\{ Q_{\mu }(t)\,:t\in [0,T] \}$ differs from its stationary counterpart $Q_{\mu }(\infty )$. To illustrate this, Fig. 1 depicts the expected value $Q_{\mu }$ in a M / M / 1 queue as a function of time for several initial workloads Q(0) for a particular setting of $\lambda$ and $\mu$. Clearly, transient behavior of $\mathbb {E}[Q_{\mu }(t)]$, for $Q(0) \ne Q_{\mu }(\infty )$, differs significantly from the steady-state mean with the same system parameters. Note that even if $Q(0) \equiv \mathbb {E}[Q_{\mu }(\infty )]$, the time-dependent mean does not coincide with the steady-state mean. Moreover, $\mathbb {E}[Q_{\mu }(t)]$ is not even a strictly increasing nor decreasing function of time. This phenomenon is a consequence of the decomposition of the transient mean into one strictly increasing, and a strictly decreasing term for $Q(0)>0$, as discussed in [3]. Nonetheless, $Q_{\mu }(t)$ converges in distribution to $Q_{\mu }(\infty )$ as $t\rightarrow \infty$, if $\mu >\lambda$.

Fig. 1

Time-dependent mean workload in a M / M / 1 queue with $\lambda = 10$ and $\mu =11$ for different initial states Q(0). The dashed line depicts $\mathbb {E}Q_{\mu }(\infty )$

Since the time horizon of our analysis is limited to $t\le T$, the process may not approach the steady-state distribution sufficiently close to appropriately use its steady-state properties for capacity allocation. To overcome this disparity, we propose a way to include the influence of this transient phase in the capacity allocation problem.

2.3 Cost structure

As mentioned before, we are interested in balancing the QoS and efficiency of the queue by choosing the optimal server speed $\mu$. The adjective optimal indicates that we intend to choose the speed according to some objective function. In our case, we conduct our analysis based on a cost function, which consists of a part accounting for the penalty for congestion in the system and a part for staffing cost. The cost value of both parts is governed by the variable $\mu$. The instantaneous cost incurred at time t equals

$$\begin{aligned} \mathbb {E}[Q_{\mu }(t)] + \alpha \mu , \end{aligned}$$

where $\alpha$ is a positive constant defining the relative staffing cost. Hence, the cost structure we apply is a combination of the transient mean of the workload process and a linear staffing cost. Accumulated and normalized over the period [0, T], the cost function on which the rest of this paper will be based equals

$$\begin{aligned} \varPi _{T}(\mu ) := \frac{1}{T}\int _0^T\left( \mathbb {E}[Q_{\mu }(t)] + \alpha \mu \, \right) \mathrm{d}t = \frac{1}{T} \int _0^T \mathbb {E}[Q_{\mu }(t)]\, \mathrm{d}t + \alpha \mu . \end{aligned}$$

(2.4)

We use shorthand notation for the normalized congestion costs:

$$\begin{aligned} C_{T}(\mu ) := \frac{1}{T}\int _0^T \mathbb {E}[Q_{\mu }(t)] \mathrm{d}t, \end{aligned}$$

(2.5)

and $C_{\infty }(\mu ) := \mathbb {E}[Q_{\mu }(\infty )]$. In order to compare the actual costs incurred over the interval [0, T] to the cost function of the queue in stationary conditions, we define

$$\begin{aligned} \varPi _{\infty }(\mu ) := C_{\infty }(\mu ) + \alpha \mu = \mathbb {E}[Q_\mu (\infty )] + \alpha \mu , \end{aligned}$$

(2.6)

which allows an explicit expression by Lemma 1. Under mild conditions on the net-input process and the distribution of the initial state, the cost functions coincide for $T\rightarrow \infty$.

Proposition 1

Let $\mu >\lambda$ and assume $\mathbb {E}[U(1)],\, \mathbb {E}[Q_\mu (0)] < \infty$. Then

$$\begin{aligned} \lim _{T\rightarrow \infty } \varPi _{T}(\mu ) = \varPi _{\infty }(\mu ). \end{aligned}$$

Rewriting (2.4) gives the relation

$$\begin{aligned} \varPi _{T}(\mu )&= \frac{1}{T}\int _0^{T} \left( \mathbb {E}[Q_{\mu }(t)] - \mathbb {E}[Q_{\mu }(\infty )] \right) \, \mathrm{d}t + \mathbb {E}[Q_{\mu }(\infty )] + \alpha \mu \nonumber \\&= \varOmega _{T}(\mu ) + \varPi _{\infty }(\mu ). \end{aligned}$$

(2.7)

Section 3 is concerned with the analysis of the correction factor $\varOmega _{T}(\mu )$.

Ultimately, we are concerned with the additional costs incurred by choosing the server speed through minimization of $\varPi _{\infty }(\mu )$ instead of $\varPi _{T}(\mu )$. Therefore, we formulate the exact and approximate optimization problems as follows

$$\begin{aligned} \mu _T^\star := \arg \min _{\mu \ge 0} \varPi _{T}(\mu ), \qquad \qquad \mu _\infty ^\star := \arg \min _{\mu \ge 0} \varPi _{\infty }(\mu ), \end{aligned}$$

(2.8)

$$\begin{aligned} \varPi _{T}^\star := \varPi _{T}(\mu _T^\star ), \qquad \qquad \varPi _{\infty }^\star := \varPi _{T}(\mu _\infty ^\star ). \end{aligned}$$

(2.9)

In Sect. 4, we turn to the comparison of $\mu _T^{\star }$ and $\mu _\infty ^\star$ as well as the optimality gap $\varPi _{\infty }^\star - \varPi _{T}^\star$.

3 Analysis of the objective function

From (2.7) it is evident that, for finding an explicit characterization of $\varPi _{T}(\mu )$, it suffices to study the term $\varOmega _T(\mu )$ in more detail. We start by stating the main result of this section, which describes the leading order behavior of $\varOmega _T(\mu )$ as T increases.

Theorem 1

Let $X_\mu (t)$ be of the form (2.1). If $\mathbb {E}[\max (Q(0),Q_\mu (\infty ))^3] < \infty$ and $u_2,u_3 < \infty$, then

$$\begin{aligned} \varOmega _T(\mu )&= \frac{\mathbb {E}[Q(0)^2] - \mathbb {E}[Q_\mu (\infty )^2]}{2T(\mu -\lambda )} + O\left( \frac{1}{T^2}\right) \\&= \frac{1}{2T(\mu -\lambda )}\left( \mathbb {E}[Q(0)^2] - \frac{\lambda ^2 u_2^2}{2(\mu -\lambda )^2} - \frac{\lambda u_3}{3(\mu -\lambda )}\right) + O\left( \frac{1}{T^2}\right) , \end{aligned}$$

for $\mu >\lambda$.

Note that this expression provides an approximation of the actual cost function $\varPi _T(\mu )$. We elaborate on the implications of this additional information on the optimization problem in Sect. 4.

In the remainder of this section, we provide a detailed description of the steps taken to obtain this outcome. We assume a fixed service rate $\mu$ throughout the analysis in this section and therefore omit the subscript $\mu$. Proofs of the intermediate results can be found in Appendix 2.

3.1 Constructing a coupling

Before starting our analysis of the correction term $\varOmega _{T}(\mu )$, we introduce some auxiliary notation. By $Q^A(t)$ we denote the workload process as described in Sect.  2.1 with initial state A and $\mathbb {E}_A$ the expectation with respect to any nonnegative random variable A, which is independent of the net-input process X. To be able to compare $\mathbb {E}[Q(t)]$ and $\mathbb {E}[Q(\infty )]$ as in $\varOmega _T(\mu )$, we will use a coupling technique. Observe that by the definition of the stationary distribution $Q(\infty ) \,{\buildrel d \over =}\,Q^{Q(\infty )}(t)$ for all $t \ge 0$ and therefore $\mathbb {E}[Q(\infty )] = \mathbb {E}_{Q(\infty )}[Q^{Q(\infty )}(t)]$. Furthermore, $\mathbb {E}[Q(t)] = \mathbb {E}_{Q(0)}[Q^{Q(0)}(t)]$. Hence, quantifying the difference between the transient and stationary mean is equivalent to comparing the workload processes of two queues starting in two different (random) states at $t=0$.

We starting our analysis for two queues starting in two deterministic states $x,y\ge 0$, respectively. At the end of our analysis, we will obtain the original form by replacing x with Q(0) and y with $Q(\infty )$.

Equation (2.2) shows that all randomness in the workload process originates from the process X(t). With this in mind, we couple the processes $Q^x(t)$ and $Q^{y}(t)$ on a sample path level by feeding both queues the same net-input process X(t) for $t\ge 0$. This allows us to compare the processes in the same probability space, so that $\mathbb {E}[Q^x(t)] - \mathbb {E}[Q^y(t)] = \mathbb {E}[Q^x(t) - Q^y(t)]$ for all $t\ge 0$. Define

$$\begin{aligned} Y^{x,y}(t) := Q^x(t) - Q^y(t) \end{aligned}$$

and

$$\begin{aligned} \varOmega _{T}^{x,y} := \frac{1}{T}\,\int _0^T \mathbb {E}\left[ Y^{x,y}(t)\right] \, \mathrm{d}t. \end{aligned}$$

A possible sample path triple for $Q^x(t)$, $Q^0(t)$ and $Y^{x,0}(t)$ is depicted in Fig. 2. As we see from this figure, $Y^{x,0}(t)$ has nice structural properties which we will exploit in the next subsection.

Fig. 2

Sample path visualization of the processes $Q^x(t)$ (solid), $Q^0(t)$ (gray) and $Y^{x,0}(t)$ (dashed)

3.2 Difference process and leading order behavior of the correction term

We further examine the difference process $Y^{x,y}(t)$ with $x>y$. Recall from (2.2),

$$\begin{aligned} Q^z(t) = \max \{ z + X(t),\, \sup _{0<s\le t} [X(t)-X(s)]\} = X(t) + \max \{ z, -\inf _{0\le s\le t} X(s)\}, \end{aligned}$$

(3.1)

for any initial state $z\ge 0$, where X(t) is a Lévy process with no negative jumps. Let $\tau ^z(w)$, $0\le w<z$, denote the first passage time of level w by the process starting in z, i.e.

$$\begin{aligned} \tau ^z(w) := \inf \left\{ t \ge 0\, |\, Q^z(t) \le w \,\right\} . \end{aligned}$$

Then it is easily seen that for all $z\ge 0$,

$$\begin{aligned} Q^z(t) = \left\{ \begin{array}{ll} z + X(t), &{} \mathrm{if }\ t<\tau ^z(0), \\ \sup _{0<s\le t} [X(t)-X(s)], &{} \mathrm{if }\ t \ge \tau ^z(0). \end{array}\right. \end{aligned}$$

Consequently,

$$\begin{aligned} Y^{x,y}(t) = \left\{ \begin{array}{ll} x - y, &{} \text {if }t< \tau ^y(0),\\ \inf _{0<s\le t} \{ x+X(s)\}, &{} \text {if }\tau ^y(0) \le t < \tau ^x(0),\\ 0, &{} \text {if }t \ge \tau ^x(0). \end{array}\right. \end{aligned}$$

(3.2)

Using this representation, we can identify

$$\begin{aligned} \varOmega ^{x,y}_T = \frac{1}{T}\,\mathbb {E}\left[ \int _0^{\tau ^x(0)\wedge T} Y^{x,y}(t) \mathrm{d}t\right] , \end{aligned}$$

where $\wedge$ denotes the minimum operator, due to the fact that $Y^{x,y}(t) = 0$ for $t\ge \tau ^x(0)$. Subsequently, we decompose $\varOmega _T^{x,y}$ into two terms:

$$\begin{aligned} \varPsi ^{x,y}_T := \frac{1}{T} \int _0^\infty \mathbb {E}[Y^{x,y}(t)]\, \mathrm{d}t \qquad \text {and} \qquad \varDelta _T^{x,y} := \varOmega _T^{x,y} - \varPsi _T^{x,y}. \end{aligned}$$

(3.3)

Note that $\varPsi _T^{x,y}$ is obtained by replacing T by $\infty$ only in the integration bound. It is customary in the literature, particularly in the area of stochastic simulation, to compare the truncated integral to its natural expansion of the integration range to a semi-infinite interval; see, for example, [6, Prop. 2.1].The truncated integral connects to the long-run average estimator of a certain performance metric, whereas the infinite integral reflects its exact expectation. The decomposition in (3.3) is insightful, because $\varPsi _T^{x,y}$ prescribes the leading order behavior of $\varOmega _T^{x,y}$, while $\varDelta _T^{x,y}$ captures the smaller order error term. In this section, we only consider $\varPsi _T^{x,y}$. Sect. 3.3 investigates the magnitude of $\varDelta _T^{x,y}$. The next preliminary result presents a useful property of $\varPsi _T^{x,y}$.

Lemma 2

Let $x>y$. If $\mathbb {E}[\tau ^x(0)]<\infty$, then

$$\begin{aligned} \varPsi ^{x,y}_T = \frac{1}{T}\,\mathbb {E}[\tau ^{y}(0)](x-y) + \varPsi ^{x-y,0}_T. \end{aligned}$$

(3.4)

This leaves us with two unknowns: $\mathbb {E}[\tau ^y(0)]$ and $\varPsi _T^{x-y,0}$. The next lemma gives an equivalent form for the latter.

Lemma 3

If $\mathbb {E}[\tau ^z(0)] < \infty$, then for all $z\ge 0$,

$$\begin{aligned} \varPsi ^{z,0}_T = \int _0^z \mathbb {E}[\tau ^w(0)]\, \mathrm{d}w. \end{aligned}$$

(3.5)

Since the term $\mathbb {E}[\tau ^z(0)]$, for several values of z, appears in many of the preliminary results, we devote our attention to this in the next subsection.

First passage time

When studying the first passage time of level $0\le w < z$, $\tau ^z(w)$, of the workload process starting in z, we first observe that $\{\tau ^z(z-w)\}_{w=0}^z$ is a spectrally positive Lévy process itself, also visible through Fig. 2. More precisely, it is a subordinator, i.e., a Lévy process whose paths are almost surely non-decreasing [18]. In order to calculate $\mathbb {E}[\tau ^z(z-w)]$, we use theory presented in [26, Section 46], although results presented there are valid for spectrally negative Lévy processes, as opposed to the absence of negative jumps in our case. Nonetheless, our setting is easily transformed into this framework by observing that $\hat{X} \equiv -X$, that is $\hat{X}(t) = -X(t)$ for all $t\ge 0$, is spectrally negative. Furthermore, let

$$\begin{aligned} \hat{\tau }^0(w) := \inf \{ t \ge 0\,:\, \hat{X}(t) \ge w\} = \inf \{ t \ge 0\,:\, z+X(t) \le z-w\} = \tau ^z(z-w). \end{aligned}$$

(3.6)

For completeness, we cite [26, Thm. 46.3].

Theorem 2

Let $\hat{X}(t)$ be a spectrally negative Lévy process with generating triplet $(-a,\sigma ,\hat{\nu })$ and $\hat{\tau }^0(y)$ its corresponding hitting time process. Define $\varUpsilon (\theta )$ for $\theta \ge 0$ as

$$\begin{aligned} \varUpsilon (\theta ) = -a\theta + \tfrac{1}{2}\sigma ^2\theta ^2 + \int _{-\infty }^0 \left( e^{\theta x}-1-\theta x\mathbf{1}_{[-1,0)}(x)\right) \, \hat{\nu }(\mathrm{d}x). \end{aligned}$$

(3.7)

Then $\varUpsilon (\theta )$ is strictly increasing and continuous, $\varUpsilon (0)=0$, and $\varUpsilon (\theta )\rightarrow \infty$ as $\theta \rightarrow \infty$. For $w\ge 0$ and $0\le u < \infty$, we have

$$\begin{aligned} \mathbb {E}\left[ \exp (-u\hat{\tau }^0(w))\right] = \exp (-w\,\varUpsilon ^{-1}(u)), \end{aligned}$$

(3.8)

where $\theta =\varUpsilon ^{-1}(u)$ is the inverse function of $u=\varUpsilon (\theta )$.

This immediately induces an expression for $\mathbb {E}[\tau ^w(0)]$ and henceforth $\varPsi ^{z,0}$.

Corollary 1

Let X(t) be a spectrally positive Lévy process defined as in (2.1) with $\mu > \lambda$. Let $\varPsi ^{z,0}_T$ as in (3.5). Then

$$\begin{aligned} \varPsi ^{z,0}_T = \frac{z^2}{2T(\mu -\lambda )}. \end{aligned}$$

Furthermore, if $x,y\ge 0$, then

$$\begin{aligned} \varPsi ^{x,y}_T = \frac{x^2-y^2}{2T(\mu -\lambda )}. \end{aligned}$$

(3.9)

Randomization

As we stated before, we easily obtain the original $\varOmega _T$ from $\varOmega _T^{x,y}$ through substitution of x and y by Q(0) and $Q(\infty )$, respectively, and taking the expectation. In the previous paragraph, we deduced an explicit expression for $\varPsi _T^{x,y}$, the leading order term for $\varOmega _T^{x,y}$. Therefore, we equivalently get an approximation for $\varOmega _T$, given by

$$\begin{aligned} \varPsi _T := \frac{1}{T} \int _0^\infty \left( \mathbb {E}[Q(t)]-\mathbb {E}[Q(\infty )] \right) \, \mathrm{d}t, \end{aligned}$$

through randomization of x and y in $\varPsi _T^{x,y}$. By combining the results in Corollary 1, Lemma 1 and Proposition 2, which is given at the end of this section, we directly prove the result in Theorem 1.

3.3 Truncation error

In order to get a better comprehension of the properties of $\varPsi _T$, we depict the value in terms of the (infinite) region between the curves $\mathbb {E}[Q(t)]$, $\mathbb {E}[Q(\infty )]$ and the vertical axis for the case $Q(0)\equiv 0$ in Fig. 3. In this figure, $\varOmega _T$ is given by the area enclosed by the two curves, the vertical axis and the line $t=T$. One can see that the main contribution to the correction term $\varOmega _T$ is given for small t. As t increases, the difference between transient and stationary mean decreases. Hence, for moderate values of T, the contribution to the integral in (3.3) is only minor compared to the contribution over the interval [0, T].

Fig. 3

Visualization of $\varOmega _T$ and $\varPsi _T$ as the area between the curves E[Q(t)], $\mathbb {E}[Q(\infty )]$ for $Q(0) = 0$

Recall the definition of $\varDelta ^{x,y}_T$ as in (3.3). As we alluded to in Sect. 3.2, we claim the contribution of $\varDelta ^{x,y}_T$ to $\varOmega _T^{x,y}$ is negligible compared to $\varPsi ^{x,y}_T$. Also note that

$$\begin{aligned} \varDelta _T := \varOmega _T - \varPsi _T = {-}\frac{1}{T} \int _T^\infty \mathbb {E}[Q(t)] - \mathbb {E}[Q(\infty )]\,\mathrm{d}t \end{aligned}$$

(3.10)

can be derived through $\varDelta ^{x,y}_T$ in a similar manner as we did for $\varPsi ^{x,y}_T$ to obtain $\varPsi _T$. To substantiate our claim, we compute an upper bound for $\varDelta ^{x,y}_T$ of order $1/T^2$. The existence of such an upper bound poses a limit on the error this tail integral contributed to the cost structure as a whole.

Proposition 2

Let $x,y\ge 0$ and $\mathbb {E}[\max (Q(0),Q_\mu (\infty ))^3] < \infty$. Then

$$\begin{aligned} |\varDelta ^{x,y}_T| \le \frac{1}{T^2}\left( \frac{\max (y,x)^3}{3(\mu -\lambda )^2}+\frac{u_2 \max (y,x)^2}{2(\mu -\lambda )^3}\right) \end{aligned}$$

and

$$\begin{aligned} |\varDelta _T| \le \frac{1}{T^2}\left( \frac{\mathbb {E}\left[ \max (Q(0),Q_\mu (\infty ))^3\right] }{3(\mu -\lambda )^2}+\frac{u_2 \mathbb {E}\left[ \max (Q(0),Q_\mu (\infty ))^2\right] }{2(\mu -\lambda )^3}\right) . \end{aligned}$$

Remark

In the case where the net-input process X is light-tailed, that is, there exists $u>0$ such that $\mathbb {E}[\mathrm{e}^{u X(1)}] < \infty$, it can be shown that the truncation error is of order $\mathrm{e}^{-\beta T}/T$ for some $\beta >0$. See Appendix 2 for details.

4 Optimization

The result in Theorem 1, characterizing the leading order behavior of $\varOmega _T(\mu )$, also reveals the behavior of $\varPi _T(\mu )$ in leading order. Namely,

$$\begin{aligned} \varPi _T(\mu ) = \varPi _\infty (\mu ) + \varPsi _T(\mu ) + O(1/T^2). \end{aligned}$$

In fact, this representation naturally gives rise to an approximation of the actual cost function:

$$\begin{aligned} \hat{\varPi }_{T}(\mu ) := \varPi _{\infty }(\mu ) + \varPsi _T(\mu ). \end{aligned}$$

(4.1)

Denote the corresponding minimizer of $\hat{\varPi }_{T}$ by

$$\begin{aligned} \hat{\mu }_T^\star := \arg \min _{\mu \ge 0} \hat{\varPi }_{T}(\mu ), \qquad \hat{\varPi }_{T}^\star := \hat{\varPi }_{T}(\hat{\mu }_T^\star ) \end{aligned}$$

(4.2)

in addition to the definitions in (2.8) and (2.9). This section is devoted to the analysis of the minimizers $\mu _T^\star$, $\hat{\mu }_T^\star$ and $\mu _\infty ^\star$, and the optimality gap for the two approximations.

Throughout this section, we assume that $u_2, u_3 <\infty$ and $\mathbb {E}[Q(0)^2] <\infty$.

By its definition in (2.6) and Lemma 1, we have an exact expression for the steady-state cost function:

$$\begin{aligned} \varPi _{\infty }(\mu ) = \frac{\lambda u_2}{2(\mu -\lambda )} + \alpha \mu . \end{aligned}$$

It is easily verified that $\varPi _{\infty }$ is strictly convex in $\mu$, for instance by observing that $\varPi _{\infty }''(\mu ) > 0$ for all $\mu > \lambda$. Therefore, $\varPi _{\infty }$ has a unique global minimizer and

$$\begin{aligned} \mu _\infty ^\star = \lambda + \sqrt{\frac{\lambda u_2}{2\alpha }}, \qquad \varPi _{\infty }^\star = \alpha \lambda + \sqrt{2\alpha \lambda u_2}. \end{aligned}$$

(4.3)

We are interested in the relation between $\mu _\infty ^\star$ and $\mu _T^\star$, and between $\hat{\mu }_T^\star$ and $\mu _T^\star$. Since $\varPi _{T}(\mu ) = \varPi _{\infty }(\mu ) + O(1/T)$ for all $\mu > \lambda$, we have pointwise convergence of the sequence $\varPi _{T}$, as well as $\hat{\varPi }_{T}$, to $\varPi _{\infty }$ for $T\rightarrow \infty$; we also expect $\mu _T^\star \rightarrow \mu _\infty ^\star$ and $\hat{\mu }_T^\star \rightarrow \mu _\infty ^\star$ for $T\rightarrow \infty$. Before proving that this convergence indeed holds, we present a result on the strict convexity of the function $\varPi _{T}$.

Lemma 4

Let $\mu \ge 0$. The function $\varPi _{T}(\mu )$ is

• convex in $\mu$, if $Q(0)\equiv x$, $T<x/\mu$ and $\sigma =0$,

• strictly convex in $\mu$, otherwise.

Building upon strict convexity of both $\varPi _T(\mu )$ and $\varPi _\infty (\mu )$ for $\mu >\lambda$, we derive the following convergence result.

Proposition 3

Let $\mu _T^\star$, $\hat{\mu }_T^\star$ and $\mu _\infty ^\star$ be as defined in (2.8) and (4.2). Then

$$\begin{aligned} \mu _T^\star \rightarrow \mu _\infty ^\star \, \qquad \text { and } \qquad \hat{\mu }_T^\star \rightarrow \mu _\infty ^\star , \end{aligned}$$

for $T\rightarrow \infty$.

The next result describes a refinement of $\mu _T^\star$ in terms of $\mu _\infty ^\star$.

Proposition 4

For T sufficiently large,

$$\begin{aligned} \mu _T^\star = \mu _\infty ^\star + \frac{\mu _\bullet }{T} + o(1/T), \end{aligned}$$

where

$$\begin{aligned} \mu _\bullet = \frac{\mathbb {E}[Q(0)^2]}{\sqrt{8\lambda u_2\alpha }} - \frac{u_3}{3 u_2} - 3\sqrt{\frac{\alpha \lambda u_2}{8}}. \end{aligned}$$

(4.4)

Based on Proposition 4, we propose a corrected staffing rule, accounting for the finite horizon:

$$\begin{aligned} \tilde{\mu }_T^\star = \left[ \mu _\infty ^\star + \frac{\mu _\bullet }{T}\right] ^+, \end{aligned}$$

(4.5)

with $\mu _\bullet$ as in (4.4). Here $[x]^+ := \max \{x,0\}$, which ensures the value of $\tilde{\mu }_T^\star$ is nonnegative and thus is a feasible solution of the optimization problem. This refined capacity allocation rule is expected to reduce the costs incurred in transient settings. However, the value of particular interest to us is the cost penalty for using either one of the approximations rather than the actual minimum $\mu _T^\star$, that is, the optimality gap. As it happens, we deduce the order of the optimality gap for $\mu _\infty ^\star$ with the help of the explicit form of $\mu _\bullet$ given in (4.4), which is stated in the next proposition. The proof is given in Appendix 3.

Proposition 5

Let $\mu _\infty ^\star$ be as in (4.3). Then,

$$\begin{aligned} \varPi _\infty ^\star - \varPi _T^\star = O(1/T^2). \end{aligned}$$

5 Numerical experiments

5.1 Influence of $\varOmega _{T}(\mu )$

We first assess the contribution of the correction to the cost function provided by Theorem 1. In other words, we investigate whether $\hat{\varPi }_{T}(\mu )$ as in (2.4) yields a significantly better fit to $\varPi _{T}(\mu )$ than $\varPi _{\infty }(\mu )$ does. Note that these three functions only differ in the costs describing the congestion. Therefore, we limit our study in this subsection to the evaluation of $C_T(\mu )$ as in (2.5) with stationary equivalent $C_{\infty }(\mu ) = \mathbb {E}[Q_{\mu }(\infty )]$. Our novel approximation hence reads

$$\begin{aligned} \hat{C}_{T}(\mu ) := C_{\infty }(\mu ) + \varOmega _{T}(\mu ), \end{aligned}$$

with $\varOmega _{T}(\mu )$ given in Theorem 1. We conduct our numerical experiments based on three models, namely:

• 1. M / M / 1 queue: U(t) is a unit rate compound Poisson process with exponentially distributed increments. We have $u_2 = 2$, $u_3=3$, so that

  $$\begin{aligned} \hat{C}_{T}(\mu ) = \frac{\lambda }{\mu -\lambda } + \frac{1}{T(\mu -\lambda )} \left( \frac{x^2}{2} - \frac{\lambda ^2}{(\mu -\lambda )^2} - \frac{\lambda }{\mu -\lambda } \right) . \end{aligned}$$

  (5.1)

• 2. $M/\mathrm{Pareto}/1$ queue: U(t) is a unit rate compound Poisson process with Pareto increments. The Pareto distribution deserves special attention due to its heavy tailed nature, having tail probability $\bar{F}(x) = (x/k)^{-\gamma }$, if $x\ge k$, and 1 otherwise. It is well-known that heavy-tailed service times lead to long relaxation time. For our purposes, we fix shape parameter $\gamma = 16/5$ and scale parameter $k=11/16$, so that $\beta = 1$, $u_2 = 121/96$, $u_3 = 1331/256$ and $u_k=\infty$ for all $k>3$. Hence,

  $$\begin{aligned} \hat{C}_{T}(\mu ) = \frac{121\lambda }{192(\mu -\lambda )} + \frac{1}{2T(\mu -\lambda )} \left( x^2 - \frac{(121\lambda /96)^2}{2(\mu -\lambda )^2} - \frac{ 1331\lambda /256 }{2(\mu -\lambda )}\right) . \end{aligned}$$

  (5.2)

• 3. Reflected Brownian motion: U(t) is Brownian motion with drift 1 and infinitesimal variance $\sigma ^2$. We have $u_2 = \sigma ^2$, $u_3=0$, so that

  $$\begin{aligned} \hat{C}_{T}(\mu ) = \frac{\lambda \sigma ^2}{2(\mu -\lambda )} + \frac{1}{2T(\mu -\lambda )} \left( x^2 - \frac{\lambda ^2\sigma ^4}{2(\mu -\lambda )^2}\right) . \end{aligned}$$

  (5.3)

In light of the equivalence relations in (2.3), we only consider the case $\lambda =1$. The cost values for general values of $\lambda$ follow by appropriate rescaling of $\mu$ and T.

For the M / M / 1 and $M/\mathrm{Pareto}/1$ queues, we obtained the function $C_{T}(\mu )$ through simulation and the results are accurate up until a 95% confidence interval of width $10^{-3}$. For reflected Brownian motion, we used the explicit distribution function given in [13] for double numerical integration. The results for several values of T and two different starting states are depicted in Figs. 4, 5 and 6. These plots also include the approximated functions $\hat{C}_{T}(\mu )$. We make a few observations based on these figures.

First, we indeed note the pointwise convergence of $\hat{C}_{T}(\mu )$ to $\hat{C}_{\infty }(\mu )$ as T grows, for all $\mu$ in all three cases. However, the difference between the stationary costs and those for small values of T can be significant. This is most clear in the plots with $x=2.5$ and when $\mu$ is close to $\lambda$, i.e., it is in heavy traffic. In these scenarios, it is evident that refinements to the stationary cost function are needed. $\hat{C}_{T}(\mu )$ does a fairly good job at providing such correction, especially for moderate values of $\mu$.

Furthermore, we note that $C_{T}(\mu )$ approaches $C_{\infty }(\mu )$ from below for $x=0$ for any value of $\mu$, while this is not strictly the case for $x>0$. $\hat{C}_{T}(\mu )$ correctly captures the sign of this correction.

Fig. 4

Comparison of exact waiting cost function $C_T(\mu )$ against corrected cost function $\hat{C}_T(\mu )$ and PSA cost function $C_\infty (\mu )$ for $T=2,5$ and 10 for the M / M / 1 queue with $\lambda =1$. a $x=0$. b $x=2.5$

Fig. 5

Comparison of exact waiting cost function $C_T(\mu )$ against corrected cost function $\hat{C}_T(\mu )$ and PSA cost function $C_\infty (\mu )$ for $T=2,5$ and 10 for the M / Pareto / 1 queue with $\lambda =1$. a $x=0$. b $x=2.5$

Fig. 6

Comparison of exact waiting cost function $C_T(\mu )$ against corrected cost function $\hat{C}_T(\mu )$ and PSA cost function $C_\infty (\mu )$ for $T=2,5$ and 10 for reflected Brownian motion with $\sigma =1$. a $x=0$. b $x=2.5$

Finally, observe that $\hat{C}_{T}(\mu )\rightarrow -\infty$ as $\mu$ approaches $\lambda$ from above. This divergence is clear from the expressions in (5.1)-(5.3). Our correction term relies on the premise that under the coupling scheme, the sample paths of the two queues starting from different states have hit with high probability. This is equivalent to stating that the ‘largest’ of the two queues is emptied at least once before time T. However, as $\mu$ approaches $\lambda$, the system enters heavy traffic, and hence the hitting time of the zero barrier is set to run off to infinity. Consequently, this causes our approximation to be inaccurate for small values of $\mu$.

5.2 Validation of corrected staffing rule

In this section we examine whether the corrected staffing rule $\tilde{\mu }_T^\star$ as in (4.5) indeed yields a significant cost reduction over the choice of $\mu _\infty ^\star$ by comparing their true costs $\varPi _{T}(\tilde{\mu }_T^\star )$ and $\varPi _{T}(\mu _\infty ^\star )$. We conduct this comparison for different values of the parameters, $\alpha$, T and starting state x through numerical experiments. The three models on which we do our calculations are the M / M / 1 queue, the M / Pareto / 1 queue and the reflected Brownian motion, as introduced in the previous subsection. We again focus on $\lambda =1$ only.

For each of the three models, we adhere to the following setup. The quality of both staffing rules is assessed for $\alpha = 0.1, 1$ and 2, resembling three modes of valuation of the QoS in the system. As a benchmark, observe that the expected workload in steady-state conditions with staffing level $\mu _\infty ^\star$ equals

$$\begin{aligned} C_\infty (\mu _\infty ^\star ) = \sqrt{\frac{\alpha \lambda u_2}{2}}. \end{aligned}$$

For each value of $\alpha$, we consider two scenarios: One in which the system starts empty, i.e., $x=0$, and one in which the initial state is double this benchmark value, thus $x=\sqrt{2\alpha \lambda u_2}$. The numerics are presented for each model separately. We discuss general conclusions drawn from these results afterwards.

M / M / 1 queue

As we discussed before, if U is a unit rate compound Poisson process with exponentially distributed increments, then $Q_{\mu }$ describes the workload process in an M / M / 1 queue. For this setting, we get

$$\begin{aligned} \mu _\infty ^\star = \lambda + \sqrt{\frac{\lambda }{\alpha }},\qquad \tilde{\mu }_T^\star = \left[ \lambda + \sqrt{\frac{\lambda }{\alpha }} + \frac{1}{T}\left( \frac{x^2}{4\sqrt{\lambda \alpha }} - 1 - \frac{3}{2} \sqrt{\lambda \alpha }\right) \right] ^+. \end{aligned}$$

Table 1 presents the actual costs corresponding to these two staffing levels for different value of x and $\alpha$.

Table 1 Comparison of costs for the M / M / 1 queue for steady-state and corrected staffing rules and relative cost improvement (r.c.i.)

M/Pareto/1 queue

In the case where the service requirements follow a Pareto distribution with shape parameter $\gamma = 16/5$, the staffing rule becomes

$$\begin{aligned} \mu _\infty ^\star = \lambda + \frac{11}{8}\sqrt{\frac{ \lambda }{3 \alpha }}, \ \tilde{\mu }_T^\star = \left[ \lambda + \frac{11}{8}\sqrt{\frac{ \lambda }{3 \alpha }} + \frac{1}{T}\left( \frac{2 x^2}{11\sqrt{\lambda \alpha /3}} - \frac{11}{8} - \frac{11\sqrt{3\lambda \alpha }}{16}\right) \right] ^+. \end{aligned}$$

The numerical results are given in Table 2.

Table 2 Comparison of costs for the $M/\mathrm{Pareto}/1$ queue for steady-state and corrected staffing rules and relative cost improvement (r.c.i.)

Just as in the results for the M / M / 1 queue, we observe a higher reduction for larger value of $\alpha$ and T. Also, again $\tilde{\mu }_T < \mu _\infty ^\star$. Hence, the conclusions for the $M/\mathrm{Pareto}/1$ queue are similar to those of the M / M / 1 queue.

Reflected Brownian motion

In the case where the input process U is Brownian motion with drift 1 and infinitesimal variance $\sigma ^2$, the steady-state staffing rule and its corrected version reduce to

$$\begin{aligned} \mu _\infty ^\star = \lambda + \sqrt{\frac{\lambda \sigma ^2}{2\alpha }}, \qquad \tilde{\mu }_T^\star = \left[ \lambda + \sqrt{\frac{\lambda \sigma ^2}{2\alpha }} + \frac{1}{2\sqrt{2}\,T}\left( \frac{x^2}{\sqrt{\lambda \alpha }\sigma } - 3\sigma \sqrt{\alpha \lambda } \right) \right] ^+. \end{aligned}$$

In Tables 3 and 4, the costs obtained through numerical evaluation are presented for several values of x, T. We also vary $\sigma$ to examine the influence of the volatility of the arrival process on the quality of the staffing rules.

The observations on the influence of $\alpha , x$ and T are similar to those for the M / M / 1 queue and the $M/\mathrm{Pareto}/1$ queue. However, here we see little improvement by the corrected staffing rule for small values of $\alpha$ for both values of x. The results in Tables 3 and 4 also suggest that the reduction is smaller for larger values of $\sigma$.

Table 3 Comparison of costs for RBM with $\sigma = 1$ for steady-state and corrected staffing rules and relative cost improvement (r.c.i.)

Table 4 Comparison of costs for RBM with $\sigma = 2$ for steady-state and corrected staffing rules and relative cost improvement (r.c.i.)

5.3 Discussion

Based upon these numerical results in Tables 1, 2, 3 and 4, we make a few remarks. The three models roughly exhibit similar behavior as T, x and $\alpha$ are varied.

Non-surprisingly, we note that $\tilde{\mu }_T$ approaches $\mu _\infty ^\star$ with increasing T, which also implies that the cost reduction achieved by the corrected staffing rule vanishes as $T\rightarrow \infty$. Also, we observe that in all scenarios examined, the cost reduction increases with $\alpha$. This can be explained through investigation of the objective function $\varPi _T$ as a function of $\mu$. Namely, for $\alpha$ small, the curve is relatively flat around the true optimum $\mu _T^\star$. Hence, in this case a moderate deviation from $\mu _T^\star$ will likely not lead to a significant cost increase. However, as $\alpha$ becomes larger, i.e., server efficiency is valued more than minimization of congestion, the curve becomes more sharp around $\mu _T^\star$, and hence more accurate approximations of $\mu _T^\star$ are required to achieve an acceptable cost level. Hence, the corrected staffing rule (4.5) proves particularly useful in these cases.

Another point we highlight is that the relative improvement is higher for $x=0$ than for $x=\sqrt{2\alpha \lambda u_2}$. Moreover, even though the initial state of the system is above the optimal equilibrium, $\tilde{\mu }_T$ is smaller than $\mu _\infty ^\star$. This is somewhat counter-intuitive. In fact, from (4.4) it follows that $\mu _\bullet$ positively contributes to the corrected staffing function if

$$\begin{aligned} \mathbb {E}[Q^2(0)] > 3\alpha \lambda u_2 + \frac{2 u_2}{3 u_3}\,\sqrt{2\alpha \lambda u_2}. \end{aligned}$$

6 Conclusion and further research

Motivated by the time-varying nature of queues in practical applications, we studied the impact that the transient phase has on traditional capacity allocation questions. By defining a cost minimization problem in which the objective function contains a correction accounting for the transient period, we identified the leading and second-order behavior of the cost function as a function of the interval length T. As a by-product, this result yields an approximation for the actual cost function, which is a refinement to its stationary counterpart. Our numerical experiments in Sect. 5.1 demonstrate the improved accuracy achieved by this approximation in a number of settings. By perturbation analysis of the optimization problem, this furthermore gives rise to a correction to the steady-state optimal capacity allocation of order 1 / T. The necessity of the refined capacity allocation level is substantiated by the numerics in Sect. 5.2, which show the cost reduction that can be achieved in a number of settings, compared to settings in which stationary metrics are used. Particularly for small values of T and large values of $\alpha$, this reduction is significant. Additionally, these results also indicate that it is relatively safe to use the stationary cost when T is moderate, or $\alpha$ is small. The latter reflects the scenario in which QoS is much more valued than service efficiency. This observation links to the flat nature of the cost function around its optimal value for $\alpha$ small, a statement on the optimality gap that we formally proved in Proposition 4.

Besides the validation of our theoretical results of Sects. 3 and 4, the numerical results also reveal some phenomena that require more investigation.

As noted, our corrected capacity allocation level $\tilde{\mu }_T^\star$ is in most cases studied less than the steady-state optimal value $\mu _{\infty }^\star$. This implies that congestion levels tend to be higher under our staffing scheme then under stationary staffing. A possible explanation for this may be the fact that the planning period under consideration is finite. Clearly, in the setting we analyzed, anything that happens after time T is neglected. Therefore, it might be beneficial from the cost perspective to end the period with a higher expected congestion level, as it does not need to be canceled out in the future. Related to this observation, it would be interesting to look at the setting in which staffing decisions need to be made in consecutive periods of equal length, in which the arrival rate changes at the start of each period. This case requires careful consideration of the correlation among the staffing decisions within the separate periods.

Another question that arises concerns the translation of our (qualitative) findings to more general queues, in particular the M / G / s queue. Whereas in our analysis the central decision variable is the server speed $\mu$, the variable of interest in multi-server queues is typically the number of servers. It may well be that similar explicit corrections to staffing levels can be deduced to account for transience. Since our analysis heavily relies on the comparability of the sample paths of two single-server queues, which is due to the equal negative drift for the two processes, another approach must be taken to tackle this extension.

The analysis and findings for the single-server queue with Lévy input presented in this paper may serve as a stepping stone for investigation of these more elaborate problems.

Appendices

Appendix 1: Proofs of Section 2

1.1 Proof of Proposition 1

Proof

We prove the limit by showing that the difference

$$\begin{aligned} \varPi _T(\mu ) - \varPi _\infty (\mu ) = \frac{1}{T} \int _0^T \mathbb {E}[Q_\mu (t)] - \mathbb {E}[Q_\mu (\infty )] \, \mathrm{d}t \end{aligned}$$

converges to zero as $T\rightarrow \infty$ for $\mu >\lambda$ fixed. The assumption $\mathbb {E}[U(1)], \mathbb {E}[Q(0)] < \infty$ implies by [4, Prop. 1] that $\mathbb {E}[Q_\mu (t)]<\infty$ for all $t\ge 0$. Following [4], we use the decomposition

$$\begin{aligned} \mathbb {E}[Q_\mu (t)] = \mathbb {E}[Q^0_\mu (t)] + \left\{ \mathbb {E}[Q_\mu (t)] - \mathbb {E}[Q_\mu ^0(t)]\right\} , \end{aligned}$$

where $Q_\mu ^0(t)$ represents the workload process if the system starts empty. From this decomposition, it is revealed that $\mathbb {E}[Q^0_\mu (t)]$ and $\left\{ \mathbb {E}[Q_\mu (t)] - \mathbb {E}[Q_\mu ^0(t)]\right\}$ are nonnegative monotonically increasing and decreasing functions of t, respectively; see [4, Prop. 2,Thm. 11]. Recall that $\mathbb {E}[Q_\mu (t)] \rightarrow \mathbb {E}[Q_\mu (\infty )]$ for $t\rightarrow \infty$ by ergodicity of the workload process for any initial state $\mathbb {E}[Q(0)]< \infty$, if $\mu >\lambda$. Hence,

$$\begin{aligned} \mathbb {E}[Q_\mu (t)]&\le \sup _t \mathbb {E}[Q_\mu ^0(t)] + \sup _t \left\{ \mathbb {E}[Q_\mu (t)] - \mathbb {E}[Q_\mu ^0(t)]\right\} \\&= \mathbb {E}[Q_\mu (\infty )] + \left\{ \mathbb {E}[Q_\mu (0)] - \mathbb {E}[Q_\mu ^0(0)]\right\} = \mathbb {E}[Q_\mu (\infty )] + \mathbb {E}[Q(0)], \end{aligned}$$

for all $t\ge 0$, which proves that the expected workload is bounded. Fix $\varepsilon >0$. By convergence of $\mathbb {E}[Q_\mu (t)]$ for $t\rightarrow \infty$, there exists a value $t^* := t^*(\varepsilon )$ such that for all $t\ge t^*$

$$\begin{aligned} \left| \mathbb {E}[Q_\mu (t)] - \mathbb {E}[Q_\mu (\infty )] \right| < \varepsilon /2. \end{aligned}$$

(7.1)

Next, set

$$\begin{aligned} T^* := T^*(\varepsilon ) = \frac{2\,t^*(\varepsilon )}{\varepsilon }\, \left( 2 \mathbb {E}[Q_\mu (\infty )] + \mathbb {E}[Q(0)]\right) . \end{aligned}$$

Then, for $T\ge \hat{T}:= \max \{ t^*,T^* \}$, we have

$$\begin{aligned} \left| \frac{1}{T} \int _0^T \mathbb {E}[Q_\mu (t)] - \mathbb {E}[Q_\mu (\infty )] \, \mathrm{d}t \right|&\le \frac{1}{T} \int _0^{t^*} \left| \mathbb {E}[Q_\mu (t)] - \mathbb {E}[Q_\mu (\infty )] \right| \, \mathrm{d}t\\&\quad + \frac{1}{T} \int _{t^*}^T \left| \mathbb {E}[Q_\mu (t)] - \mathbb {E}[Q_\mu (\infty )] \right| \, \mathrm{d}t \\&\le \frac{1}{T} \int _0^{t^*} \mathbb {E}[Q_\mu (t)] + \mathbb {E}[Q_\mu (\infty )] \, \mathrm{d}t + \frac{1}{T} \int _{t^*}^T \frac{\varepsilon }{2}\, \mathrm{d}t \\&< \frac{t^*}{T}\,\left( 2 \mathbb {E}[Q_\mu (\infty )] + \mathbb {E}[Q(0)]\right) + \frac{T-t^*}{T} \,\frac{\varepsilon }{2}\\&< \frac{t^*}{T^*} \,\left( 2 \mathbb {E}[Q_\mu (\infty )] + \mathbb {E}[Q(0)]\right) + \frac{\varepsilon }{2} = \varepsilon . \end{aligned}$$

Hence, for any choice of $\varepsilon >0$ we can find a value $\hat{T}$ such that $\varPi _{\hat{T}}(\mu )$ approaches $\varPi _\infty (\mu )$ within distance $\epsilon$, which proves the limit. $\square$

Appendix 2: Proofs of Section 3

1.1 Proof of Lemma 2

Proof

Using the representation in (3.2), we write

$$\begin{aligned} \varPsi ^{x,y}_T&= \frac{1}{T}\int _0^{\infty } \mathbb {E}[Y^{x,y}(t)]\mathrm{d}t \\&= \frac{1}{T}\,\mathbb {E}\left[ \int _0^{\tau ^y(0)}Y^{x,y}(t)\right] \mathrm{d}t + \frac{1}{T}\,\mathbb {E}\left[ \int _{\tau ^y(0)}^{\tau ^x(0)} Y^{x,y}(t) \mathrm{d}t\right] \\&\quad + \frac{1}{T}\,\mathbb {E}\left[ \int _{\tau ^y(0)}^\infty Y^{x,y}(t) \mathrm{d}t\right] ,\\&= \frac{1}{T}\,\mathbb {E}\left[ \int _0^{\tau ^y(0)}(x-y) \mathrm{d}t\right] + \frac{1}{T}\,\mathbb {E}\left[ \int _{\tau ^y(0)}^{\tau ^x(0)} Y^{x,y}(t) \mathrm{d}t\right] \\&= \frac{1}{T}\,\mathbb {E}[\tau ^y(0)](x-y) + \frac{1}{T}\,\mathbb {E}\left[ \int _{\tau ^y(0)}^{\tau ^x(0)} Y^{x,y}(t) \mathrm{d}t\right] . \end{aligned}$$

By (3.2) and the strong Markov property holding for Lévy processes [5], observe that $Y^{x-y,0}(t) \,{\buildrel d \over =}\,Y^{x,y}(\tau ^y(0)+t)$, whereby

$$\begin{aligned} \frac{1}{T}\,\mathbb {E}\left[ \int _{\tau ^y(0)}^{\tau ^x(0)} Y^{x,y}(t)\,\mathrm{d}t\right] = \frac{1}{T}\,\mathbb {E}\left[ \int _{0}^{\tau ^{x-y}(0)} Y^{x-y,0}(t) \mathrm{d}t\right] = \varPsi ^{x-y,0}_T, \end{aligned}$$

which completes the proof. $\square$

1.2 Proof of Lemma 3

Proof

Note that $Y^{z,0}(t)$ and $\tau ^z(w)$ are intimately related. Namely, due to the fact that X has no negative jumps,

$$\begin{aligned} \{ \tau ^z(w) \le t\} = \{Y^{z,0}(t) \le w \}. \end{aligned}$$

In fact, $Y^{z,0}(\tau ^z(w)) = w$, which implies that $\tau ^z$ is a right inverse for $Y^{z,0}(t)$. Therefore, the following equality holds:

$$\begin{aligned} \int _0^{\tau ^z(0)} Y^{z,0}(t)\, \mathrm{d}t = \int _0^z \tau ^z(w)\, \mathrm{d}w, \end{aligned}$$

which implies with the help of Fubini’s theorem that

$$\begin{aligned} \varPsi ^{z,0}_T = \frac{1}{T}\,\int _0^z \mathbb {E}[\tau ^z(w)]\, \mathrm{d}w = \frac{1}{T}\,\int _0^z \mathbb {E}[\tau ^{z-w}(0)]\,\mathrm{d}w =\frac{1}{T}\, \int _0^z \mathbb {E}[\tau ^{w}(0)] \,\mathrm{d}w. \end{aligned}$$

$\square$

1.3 Proof of Corollary 1

Proof

From (3.8),

$$\begin{aligned} \mathbb {E}\left[ \hat{\tau }^0(w)\right] = -\tfrac{\mathrm{d}}{\mathrm{d}u} \left. \mathbb {E}\left[ \exp (-u\,\hat{\tau }^0(w))\right] \right| _{u=0} = w\left. \frac{\mathrm{d}}{\mathrm{d}u} \varUpsilon ^{-1}(u)\right| _{u=0}. \end{aligned}$$

(8.1)

Since $\varUpsilon (\theta )$ is strictly increasing and $\varUpsilon (0)=0$, we get $\varUpsilon ^{-1}(0)=0$ and

$$\begin{aligned} \left. \tfrac{\mathrm{d}}{\mathrm{d}u}\varUpsilon ^{-1}(u)\right| _{u=0} = \frac{1}{\varUpsilon '(\varUpsilon ^{-1}(0))} = \{ \varUpsilon '(0) \}^{-1}. \end{aligned}$$

Furthermore,

$$\begin{aligned} \varUpsilon '(\theta )= & {} -a+ \sigma ^2\theta + \int _{-\infty }^0 (x\, e^{\theta x} - x\mathbf{1}_{[-1,0)}(x)) \hat{\nu }(\mathrm{d}x)\\= & {} -a + \sigma ^2\theta - \int _0^\infty (y\, e^{-\theta y} - y\mathrm{1}_{(0,1]}(y)) \nu (\mathrm{d}y). \end{aligned}$$

Thus, $\varUpsilon '(0) = -\mathbb {E}[X(1)] = \mu -\lambda$ and $\mathbb {E}[\hat{\tau }^0(w)] = w/(\mu -\lambda )$. By (3.5) and (3.6), we deduce that

$$\begin{aligned} \varPsi ^{z,0}_T = \frac{1}{T}\, \int _0^z \mathbb {E}[\tau ^w(0)] \,\mathrm{d}w = \frac{1}{T}\, \int _0^z \mathbb {E}[\hat{\tau }^0(w)] \mathrm{d}w = \frac{z^2}{2T(\mu -\lambda )}. \end{aligned}$$

For $x>y$, we use Lemma 2 to conclude

$$\begin{aligned} \varPsi ^{x,y}_T = \frac{y(x-y)}{T(\mu -\lambda )} + \frac{(x-y)^2}{2 T(\mu -\lambda )} = \frac{x^2-y^2}{2T(\mu -\lambda )}. \end{aligned}$$

The result for $x<y$ follows directly by the observation $\varPsi ^{x,y}_T = -\varPsi _T^{x,y}$. $\square$

1.4 Proof of Proposition 2

Proof

To derive the upper bound for $\varDelta ^{x,y}_T$, we apply the same coupling argument as described in Sect. 3. Let us assume without loss of generality that $x>y$. In this case,

$$\begin{aligned} |\varDelta ^{x,y}_T| = \frac{1}{T} \int _T^\infty \mathbb {E}[Q^x(t)-Q^y(t)]\, \mathrm{d}t \le \frac{1}{T}\int _T^\infty \mathbb {E}[Q^x(t)-Q^0(t)]\, \mathrm{d}t. \end{aligned}$$

By the decomposition in (3.2),

$$\begin{aligned} \int _T^\infty \mathbb {E}\left[ Q^x(t) - Q^0(t)\right] \, \mathrm{d}t&= \int _T^\infty \mathbb {E}\left[ \left( x+\inf _{s\le t} X(s)\right) \mathbf 1 _{\{\tau ^x(0)>t\}}\right] \, \mathrm{d}t \\&= \int _T^\infty \int _0^x P( x-u + \inf _{s\le t}X(s)> 0) \, \mathrm{d}u\, \mathrm{d}t \\&= \int _T^\infty \int _0^x P( \tau ^{x-u}(0) > t ) \,\mathrm{d}u\,\mathrm{d}t \\&\le \int _T^\infty \int _0^x \frac{\mathbb {E}[\tau ^{x-u}(0)^2]}{t^2}\,\mathrm{d}u \,\mathrm{d}t\\&= \int _0^x \int _T^\infty \frac{\mathbb {E}[\tau ^{x-u}(0)^2]}{t^2}\,\mathrm{d}t \,\mathrm{d}u = \int _0^x \frac{\mathbb {E}[\tau ^{w}(0)^2]}{T}\,\mathrm{d}w. \end{aligned}$$

We obtain $\mathbb {E}[\tau ^w(0)^2]$ with the help of its Laplace transform in (3.8). Namely,

$$\begin{aligned} \mathbb {E}[\tau ^w(0)^2]= & {} \left. \tfrac{\mathrm{d}^2}{\mathrm{d}u^2}\mathbb {E}[\exp (-u \tau ^w(0))]\right| _{u=0} \\= & {} w^2\,\left( \left. \tfrac{\mathrm{d}}{\mathrm{d}u}\varUpsilon ^{-1}(u)\right| _{u=0}\right) ^2 - w\left. \tfrac{\mathrm{d}^2}{\mathrm{d}u^2}\varUpsilon ^{-1}(u)\right| _{u=0}. \end{aligned}$$

As in the previous subsection, we have $\left. \tfrac{\mathrm{d}}{\mathrm{d}u}\varUpsilon ^{-1}(u)\right| _{u=0} = (\mu -\lambda )^{-1}$, and

$$\begin{aligned} \left. \tfrac{\mathrm{d}^2}{\mathrm{d}u^2}\varUpsilon ^{-1}(u)\right| _{u=0} = {-}\frac{\varUpsilon ''(\varUpsilon ^{-1}(0))}{\varUpsilon '(\varUpsilon ^{-1}(0))^3} = {-}\frac{\varUpsilon ''(0)}{\varUpsilon '(0)^3}. \end{aligned}$$

Since $\varUpsilon '(0) = \mu -\lambda$ and

$$\begin{aligned} \varUpsilon ''(0) = \sigma ^2 + \int _0^\infty x^2\,\nu (\mathrm{d}x) = u_2, \end{aligned}$$

we conclude

$$\begin{aligned} \mathbb {E}[\tau ^w(0)^2] = \frac{w^2}{(\mu -\lambda )^2} + \frac{u_2w}{(\mu -\lambda )^3}, \end{aligned}$$

so that

$$\begin{aligned} |\varDelta ^{x,y}_T| \le \frac{1}{T^2}\int _0^x \frac{w^2}{(\mu -\lambda )^2} + \frac{u_2w}{(\mu -\lambda )^3} \mathrm{d}w = \frac{1}{T^2}\left( \frac{x^3}{3(\mu -\lambda )^2}+\frac{u_2 x^2}{2(\mu -\lambda )^3}\right) . \end{aligned}$$

(8.2)

For general $x,y\ge 0$,

$$\begin{aligned} |\varDelta ^{x,y}_T| \le \frac{1}{T^2}\left( \frac{\max (y,x)^3}{3(\mu -\lambda )^2}+\frac{u_2 \max (y,x)^2}{2(\mu -\lambda )^3}\right) . \end{aligned}$$

As a direct consequence,

$$\begin{aligned} |\varDelta _T| \le \frac{1}{T^2}\left( \frac{\mathbb {E}\left[ \max (Q(0),Q_\mu (\infty ))^3\right] }{3(\mu -\lambda )^2}+\frac{u_2 \mathbb {E}\left[ \max (Q(0),Q_\mu (\infty ))^2\right] }{2(\mu -\lambda )^3}\right) . \end{aligned}$$

$\square$

Remark

Observe that if X is light-tailed, that is $\mathbb {E}[\exp \{ -\theta X(1) \}] = \mathbb {E}[\exp \{\kappa (\theta )\}] < \infty$ for some $\theta <0$, then $\varUpsilon (\theta )$ as in (3.8) has an analytic continuation in the negative half-plane, and in this region $\varUpsilon (\theta )<0$. Consequently, we can replace the upper bound on the tail probability of $\tau ^{x-u}(0)$ by

$$\begin{aligned} \mathbb {P}\left( \tau ^{x-u}(0)> t\right) = \mathbb {P}\left( e^{\beta \tau ^{x-u}(0)} > e^{\beta t} \right) \le e^{-\beta t} \, e^{ (x-u)\varUpsilon ^{-1}(-\beta )}, \end{aligned}$$

for some $\beta > 0$, so that

$$\begin{aligned} \int _T^\infty \mathbb {E}[Q^x(t) - Q^0(t)]\, \mathrm{d}t \le e^{-\beta T}\, \frac{e^{x\varUpsilon ^{-1}(-\beta )}-1}{\beta \,\varUpsilon ^{-1}(-\beta )}. \end{aligned}$$

Along similar lines, we deduce

$$\begin{aligned} |\varDelta ^{x,y}_T| \le \frac{ e^{-\beta T}}{T}\, \frac{e^{x\varUpsilon ^{-1}(-\beta )} + e^{y\varUpsilon ^{-1}(-\beta )} -2}{\beta \,\varUpsilon ^{-1}(-\beta )} \end{aligned}$$

and

$$\begin{aligned} |\varDelta _T| \le \frac{e^{-\beta T}}{T}\, \frac{\mathbb {E}\left[ e^{Q(0)\varUpsilon ^{-1}(-\beta )}\right] + \mathbb {E}\left[ e^{Q_\mu (\infty )\varUpsilon ^{-1}(-\beta )}\right] -2}{\beta \,\varUpsilon ^{-1}(-\beta )}, \end{aligned}$$

assuming that $\mathbb {E}[\mathrm{e}^{-y Q(0)}] < \infty$ for all $y>0$. The condition $\mathbb {E}[\mathrm{e}^{Q_\mu (\infty )\varUpsilon ^{-1}(-\beta )}]<\infty$ follows from Lemma 1. Hence, the error decays exponentially fast for light-tailed input processes.

Appendix 3: Proofs of Section 4

1.1 Proof of Lemma 4

Proof

Since the term $\alpha \mu$ is convex, the strictness should come from the term $\frac{1}{T} \int _0^T \mathbb {E}[Q_\mu (t)]\,\mathrm{d}t$. Furthermore, observe that if a function $f_\mu (t)$ is convex for all $t\ge 0$, and strictly convex for all $t\ge \varepsilon$ for some $\varepsilon \in [0,T)$, i.e., for any $\mu _1,\mu _2>0$ and $a\in (0,1)$

$$\begin{aligned} a\, f_{\mu _1}(t) + (1-a) f_{\mu _2}(t) > f_{a\mu _1+(1-a)\mu _2}(t), \end{aligned}$$

then

$$\begin{aligned}&a\int _0^T \, f_{\mu _1}(t)\, \mathrm{d}t + (1-a)\int _0^T f_{\mu _2}(t) \mathrm{d}t = \int _0^T \, a f_{\mu _1}(t) + (1-a)f_{\mu _2}(t) \mathrm{d}t\\&\quad = \int _0^\varepsilon a f_{\mu _1}(t) + (1-a)f_{\mu _2}(t) \mathrm{d}t + \int _\varepsilon ^T \, a f_{\mu _1}(t) + (1-a)f_{\mu _2}(t) \mathrm{d}t\\&\quad > \int _0^\varepsilon f_{a\mu _1+(1-a)\mu _2}(t) \, \mathrm{d}t + \int _\varepsilon ^T f_{a\mu _1+(1-a)\mu _2}(t) \, \mathrm{d}t = \int _0^T f_{a\mu _1+(1-a)\mu _2}(t) \, \mathrm{d}t. \end{aligned}$$

Hence, it suffices to prove the convexity of $\mathbb {E}[Q_\mu (t)]$ as a function of $\mu$ for all $t\ge 0$, and strict convexity for $t\ge \varepsilon$ for some $\varepsilon \in [0,T)$.

Let $\tau ^x_{\mu }(0)$ denote the first passage time of level 0 in the process $Q_\mu$ with $Q(0)=x$. Then,

$$\begin{aligned} Q_\mu (t)&= U(t)-\mu t + \max \left\{ x , -\inf _{s\le t} [ U(s)-\mu s] \right\} \end{aligned}$$

(9.1)

$$\begin{aligned}&= \left\{ \begin{array}{ll} x+U(t)-\mu t, &{} \text {if } t<\tau ^x_{\mu }(0),\\ U(t)-\mu t -\inf _{s\le t} [ U(s)-\mu s], &{} \text {if } t\ge \tau ^x_{\mu }(0) , \end{array}\right. \end{aligned}$$

(9.2)

where

$$\begin{aligned} \tau ^x_{\mu }(0) := \inf \{ t \ge 0\,:\, x+U(t)-\mu t \le 0\} \end{aligned}$$

and U(t) is a spectrally positive Lévy process. Fix $\mu _1, \mu _2>0$ and $a\in (0,1)$. Define $\mu _3 := a\mu _1+(1-a)\mu _2$, and

$$\begin{aligned} D(t) := a Q_{\mu _1}(t) + (1-a) Q_{\mu _2}(t) - Q_{\mu _3}(t). \end{aligned}$$

In order to prove strict convexity, we have to show that $D(t) \ge 0$ for all $t\ge 0$, thereby implying $\mathbb {E}[D(t)] \ge 0$, i.e., convexity, for all $t\ge 0$, and $D(t)>0$ with positive probability for $t\in [\varepsilon ,T]$, for some $\varepsilon \in [0,T)$.

We distinguish two cases: $x>0$ and $x=0$.

Case $x>0$. We start by noticing that if $Q_{\mu _1}$, $Q_{\mu _2}$ and $Q_{\mu _3}$ experience the same input process U(t), then by absence of negative jumps in U(t), it holds that

$$\begin{aligned} \tau ^x_{\mu _2}(0)< \tau ^x_{\mu _3}(0) < \tau ^x_{\mu _1}(0). \end{aligned}$$

(9.3)

We use shorthand notation

$$\begin{aligned} I_k(t) := \inf _{0\le s\le t}[U(s)-\mu _k s], \end{aligned}$$

for $k=1,2,3$. Using representation (9.2) of the workload process, we obtain

$$\begin{aligned} D(t) = \left\{ \begin{array}{ll} 0, &{} \text {if } t< \tau ^x_{\mu _2}(0),\\ -(1-a)\left( x+I_2(t) \right) , &{} \text {if } \tau ^x_{\mu _2}(0) \le t< \tau ^x_{\mu _3}(0),\\ a x - (1-a)I_2(t) + I_3(t), &{} \text {if } \tau ^x_{\mu _3}(0) \le t < \tau ^x_{\mu _1}(0),\\ - a I_1(t) - (1-a) I_2(t) + I_3(t), &{} \text {if } t \ge \tau ^x_{\mu _1}(0). \end{array} \right. \end{aligned}$$

This partition allows us to spot when strict convexity can occur. Note that by definition $t \ge \tau ^x_{\mu _2}(0)$, $I_2(t) = \inf _{0\le s\le t} [U(s)-\mu _2 s] \le -x$, so that $D(t)\ge 0$ if $\tau ^x_{\mu _2}(0) \le t < \tau ^x_{\mu _3}(0)$. Moreover, by subadditivity of the infimum,

$$\begin{aligned} I_3(t)= & {} \inf _{0\le s\le t}[U(s)-\mu _3s] = \inf _{0\le s\le t}\left[ a(U(s)-\mu _1s)+(1-a)(U(s)-\mu _2s)\right] \\\ge & {} a \inf _{0\le s\le t}[U(s)-\mu _1s] + (1-a) \inf _{0\le s\le t}[U(s)-\mu _2s]\\= & {} a I_1(t) + (1-a) I_2(t), \end{aligned}$$

and hence $D(t)\ge 0$ for $t \ge \tau ^x_{\mu _1}(0)$. Using the same argument, we deduce

$$\begin{aligned} ax - (1-a)I_2(t) + I_3(t) \ge a x - (1-a) I_2(t) + a I_1(t) + (1-a) I_2(t) = a(x + I_1(t)). \end{aligned}$$

In particular, for $t < \tau ^x_{\mu _1}(0)$ this value is strictly positive. As a result, $D(t)\ge 0$ for all $t\ge 0$. On top of that $D(t) > 0$ for $t\in [\tau ^x_{\mu _3}(0),\tau ^x_{\mu _1}(0))$. Accordingly, the latter implies strict positivity of $\mathbb {E}D(t)$, and therefore strict convexity of $\mathbb {E}Q_\mu (t)$, if the event $\{\tau ^x_{\mu _3}(0)\le t< \tau ^x_{\mu _1}(0)\}$ occurs with positive probability. That is,

$$\begin{aligned} P(D(t)>0)&\ge P\left( a(x+I_1(t))\mathbf 1 _{\{\tau ^x_{\mu _3}(0)\le t< \tau ^x_{\mu _1}(0)\}}> 0 \right) \nonumber \\&= P\left( x+ I_1(t)> 0 , \tau ^x_{\mu _3}(0)\le t< \tau ^x_{\mu _1}(0)\right) \nonumber \\&= P\left( x+ I_1(t)> 0 | \tau ^x_{\mu _3}(0)\le t< \tau ^x_{\mu _1}(0)\right) P\left( \tau ^x_{\mu _3}(0)\le t< \tau ^x_{\mu _1}(0)\right) \nonumber \\&= P\left( \tau ^x_{\mu _3}(0)\le t < \tau ^x_{\mu _1}(0)\right) = P(\tau ^x_{\mu _3}(0)\le t) - P(\tau ^x_{\mu _1}(0) \le t) > 0, \end{aligned}$$

(9.4)

by the stochastic dominance in (9.3). To ensure the strict inequality in (9.4), we have to enforce the condition

$$\begin{aligned} P(\tau ^x_{\mu _1}(0)<T) > 0. \end{aligned}$$

(9.5)

Remark

An example illustrating the need for this condition is the case in which U(t) is a compound Poisson process and $T< x/\mu _2 < x/\mu _1$. Then

$$\begin{aligned} Q_{\mu _k}(t) = x + U(t) - \mu _k t, \end{aligned}$$

for all $t\in [0,T]$, since $U(t)\ge 0$ and therefore $\tau ^x_{\mu _1}(0) > T$. Consequently, for all $a\in (0,1)$,

$$\begin{aligned} a\,Q_{\mu _1} + (1-a)\,Q_{\mu _2}(t) = Q_{\mu _3}(t), \end{aligned}$$

proving only convexity of $\mathbb {E}Q_{\mu }(t)$ and subsequently $\int _0^T \mathbb {E}[Q_\mu (t)]\,\mathrm{d}t$. In the case $\sigma >0$, the probability in (9.5) is necessarily positive.

The case $x=0$. By the fact that $\tau _{\mu }(0) = 0$ for all $\mu >0$, proving that $D(t)>0$ in the case $x=0$ reduces to showing that the probability of

$$\begin{aligned} D(t) = a I_1(t) + (1-a) I_2(t) - I_3(t)>0 \end{aligned}$$

happening is positive for all $t>0$. Define

$$\begin{aligned} t_0 := \inf \{ t> 0\, :\, U(t) > 0 \}, \end{aligned}$$

and

$$\begin{aligned} \tilde{\tau }_\mu := \inf \{ t > t_0\,: U(t) - \mu t \le 0\}. \end{aligned}$$

We note that $t_0$, as defined above, also defines the epoch of the start of a new excursion of the reflection $Q_\mu$ for all $\mu >0$. Namely,

$$\begin{aligned} U(s) \le 0 \quad\Rightarrow & {} \quad U(s) - \mu s \le -\mu s \qquad \text {for all }0\le s< t_0\\\Rightarrow & {} \inf _{0\le s< t_0} [U(s)-\mu s] \le -\mu t_0 \\\Rightarrow & {} U(t_0) - \mu t_0 - \inf _{0\le s < t_0} [U(s)-\mu s] \ge U(t_0) > 0. \end{aligned}$$

Then $Q_\mu (t_0-) = 0$ for all $\mu >0$. By virtue of the strong Markov Property, note that $Q_\mu (t_0+t) \,{\buildrel d \over =}\,Q_\mu (t)$. Hence, we assume without loss of generality that $t_0=0$. Again, we have a stochastic dominance relation similar to (9.3):

$$\begin{aligned} \tilde{\tau }_{\mu _2}< \tilde{\tau }_{\mu _3} < \tilde{\tau }_{\mu _1}, \end{aligned}$$

for all $\mu _1<\mu _3<\mu _2$. Then

$$\begin{aligned} D(t) \,{\buildrel d \over =}\,\left\{ \begin{array}{ll} 0, &{} \text {if } t< \tilde{\tau }_{\mu _2},\\ -(1-a)I_2(t), &{} \text {if } \tilde{\tau }_{\mu _2} \le t< \tilde{\tau }_{\mu _3},\\ (1-a)I_2(t) + I_3(t), &{} \text {if }\tilde{\tau }_{\mu _3} \le t < \tilde{\tau }_{\mu _1},\\ - a I_1(t) - (1-a) I_2(t) + I_3(t), &{} \text {if }t\ge \tilde{\tau }_{\mu _1}. \end{array} \right. \end{aligned}$$

Clearly, $D(t)\ge 0$ for all $t\ge 0$ and

$$\begin{aligned} -(1-a)I_2(t) + I_3(t) \ge a I_1(t) > 0, \end{aligned}$$

for $\tilde{\tau }_{\mu _3} \le t < \tilde{\tau }_{\mu _1}$. Hence, in a similar manner to (9.4),

$$\begin{aligned} P(D(t)>0)&\ge P\left( aI_1(t)\mathbf 1 _{\{\tilde{\tau }_{\mu _3} \le t< \tilde{\tau }_{\mu _1}\}}> 0 \right) \nonumber \\&= P\left( I_1(t)> 0 , \tilde{\tau }_{\mu _3} \le t< \tilde{\tau }_{\mu _1}\right) \nonumber \\&= P\left( I_1(t)> 0 | \tilde{\tau }_{\mu _3} \le t< \tilde{\tau }_{\mu _1}\right) P\left( \tilde{\tau }_{\mu _3} \le t< \tilde{\tau }_{\mu _1}\right) \nonumber \\&= P\left( \tilde{\tau }_{\mu _3} \le t < \tilde{\tau }_{\mu _1}\right) = P(\tilde{\tau }_{\mu _3}\le t) - P(\tilde{\tau }_{\mu _1} \le t) > 0, \end{aligned}$$

(9.6)

The last inequality is satisfied if $P(\tilde{\tau }_{\mu _1} < T) >0$, which is equivalent to $P( U(T) - \mu T \le 0 ) >0$, a condition that is clearly true for all our choices of U. In conclusion, for $x=0$, $\mathbb {E}[D(t)] >0$ and therefore $\mathbb {E}[Q_\mu (t)]$ is a strictly convex function of $\mu$. $\square$

1.2 Proof of Proposition 3

The proof of the proposition relies on the following auxiliary lemma, for which we include the proof for completeness.

Lemma 5

Consider the sequence for functions $f_n:\, [x_0,\infty ) \rightarrow \mathbb {R}$ and let $f: [x_0,\infty ) \rightarrow \mathbb {R}$ be the pointwise limit for some $x_0\in \mathbb {R}$. Assume f and $f_n$ are strictly convex for all n. Furthermore, let $f(y) \rightarrow \infty$ for both $y\rightarrow x_0^+$ and $y\rightarrow \infty$. If $x_n$ and x are the minimizers for $f_n$ and f, respectively, then $x_n\rightarrow x$ for $n\rightarrow \infty$.

Proof

We start by showing that the sequence $x_n$ is bounded. Fix $u_l, u_r$ such that $x_0<u_l< x < u_r$. We claim that there exists a $N\in \mathbb {N}$ such that $x_n\in [u_l,u_r]$ for all $n \ge N$. First, we prove the upper bound on $x_n$. For any strictly convex function h with minimizer $x_h$, the following statement holds true:

$$\begin{aligned} x_h < u_r \quad \Leftrightarrow \quad h \text { is strictly increasing at } u_r. \end{aligned}$$

(9.7)

The first implication follows from observing that $h(x_h) < h(y)$ for all $y> x^*$ and the definition of convexity:

$$\begin{aligned} 0<\frac{h(u_r)-h(x_h)}{u_r-x_h} \le \frac{h(u_r+\delta )-h(u_r)}{\delta }, \end{aligned}$$

for all $\delta >0$, so that $h(u_r)<h(u_r+\delta )$, i.e., h is increasing at $u_r$. The converse follows immediately by observing that $h(u_r) < h(u_r+\delta )$ for all $\delta >0$, so that $x_h < u_r$. Next, we show that $f_n$ must be increasing at $u_r$ for n sufficiently large. By pointwise convergence of $f_n$, we have

$$\begin{aligned} \lim _{n\rightarrow \infty } [f_n(u_r+\delta ) - f_n(u_r)] = f(u_r+\delta ) - f(u_r). \end{aligned}$$

Let $w_r:= f(u_r+\delta ) - f(u_r)>0$. Then

$$\begin{aligned} \exists N_r \in \mathbb {N}:\, \forall n\ge N_r:\, |[f_n(u_r+\delta ) - f_n(u_r)] - [f(u_r+\delta )-f(u_r)] | < w_r/2. \end{aligned}$$

Hence for $n\ge N_r$,

$$\begin{aligned} f(u_r+\delta )-f(u_r) - w_r/2< & {} f_n(u_r+\delta ) - f_n(u_r)< f(u_r+\delta )-f(u_r) + w_r/2\\ \Rightarrow 0< & {} w_r/2 < f_n(u_r+\delta ) - f_n(u_r). \end{aligned}$$

Hence, by (9.7), $x_n < u_r$ for sufficiently large n. Similarly, we argue

$$\begin{aligned} x_h > u_l \quad \Leftrightarrow \quad h \text { is strictly decreasing at } u_l, \end{aligned}$$

for any strictly convex function h with minimizer $x_h$. Note that $x_h > u_l$ implies $h(x_h) - h(u_l) < 0$ and for all $\delta >0$ we get by strict convexity

$$\begin{aligned} \frac{h(u_l)-h(u_l-\delta )}{\delta }< \frac{h(x_h)-h(u_l)}{x_h-u_l} < 0, \end{aligned}$$

by which $h(u_l-\delta )>h(u_l)$, i.e., h is decreasing in $u_l$. Moreover, if h is decreasing at $u_l$, then it is decreasing for all $y < u_l$, by arguments similar to the above. Therefore, $h(u_l-\delta )> h(u_l)$ for all $\delta >0$ and it must hold that $x_h>u_l$. Define $f(u_l) - f(u_l-\delta ) :=w_l < 0$, then, again by pointwise convergence, we have that

$$\begin{aligned} \exists N_l \in \mathbb {N}:\, \forall n\ge N_l:\, |[f_n(u_l) - f_n(u_l-\delta )] - [f(u_l)-f(u_l-\delta )] | < w_l, \end{aligned}$$

whereupon

$$\begin{aligned} f_n(u_l) - f_n(u_l-\delta )< f(u_l) - f(u_l-\delta ) + w_l = 2w_l < 0. \end{aligned}$$

Hence, for sufficiently large n, we also have $x_n > u_l$. Fix $N = \max \{N_l,N_r\}$. Then, for $n\ge N$, $x_n\in ( u_l,u_r)$. That is, the sequence $x_n$ is bounded. Therefore, by the theorem of Bolzano–Weierstrass, $x_n$ has to have a convergent subsequence. That is, there exists a sequence $n_k$ such that $n_k \rightarrow \infty$ and $x_{n_k}\rightarrow a$ as $k\rightarrow \infty$ for some $a \in [u_l,u_r]$. We prove that every subsequence must converge to x by contradiction. Suppose there exists a subsequence $n_k$ such that $x_{n_k}\rightarrow a\ne x$. Since, $x_n\in [u_l,u_r]$ for $n\ge N$, we may restrict our attention on the sequence of functions $\hat{f}_n:[u_l,u_r] \rightarrow \mathbb {R^+}$, consisting of the original function $f_n$ restricted to the domain $[u_l,u_r]$. To be precise, $x_n = \arg \min _y f_n(y) = \arg \min _y \hat{f}_n(y)$ for $n\ge N$. Because $\hat{f}_n$ and $\hat{f}$ are bounded, we furthermore have $\hat{f}_n \rightarrow \hat{f}$ uniformly.

Fix $\varepsilon >0$. By uniform convergence, there exists $K_{0} \in \mathbb {N}$ such that

$$\begin{aligned} | \hat{f}_{n_k}( y ) - \hat{f}( y)| < \varepsilon /2,\quad \forall k\ge K_0,\ y \in [u_l,u_r]. \end{aligned}$$

Also, because $\hat{f}$ is convex, it is continuous, so that there exists $\delta := \delta (\varepsilon )$ such that

$$\begin{aligned} |z-y|< \delta \quad \Rightarrow \quad |\hat{f}(z) - \hat{f}(y)| < \varepsilon /2. \end{aligned}$$

Let $K_1$ be such that $|x_{n_k}-a| < \delta$ for all $k\ge K_1$. Then for $k \ge K= \max \{K_0,K_1\}$ this implies

$$\begin{aligned} |f_{n_k}(x_{n_k}) - f(a)|&= |\hat{f}_{n_k}(x_{n_k}) - \hat{f}(a)| \\&\le |\hat{f}_{n_k}(x_{n_k}) - \hat{f}(x_{n_k})| + | \hat{f}(x_{n_k}) - f(a)| < \varepsilon /2 + \varepsilon /2 = \varepsilon . \end{aligned}$$

Hence, we conclude $\lim _{k\rightarrow \infty } \hat{f}_{n_k}(x_{n_k}) = f(a)$. Therefore,

$$\begin{aligned} \limsup _{n\rightarrow \infty } f_n(x_n) \ge f(a) > f(x), \end{aligned}$$

by minimality of x. However, $f_n(x_n) \le f_n(x)$, which implies $\limsup _{n\rightarrow \infty } f_n(x_n) \le \lim _{n\rightarrow \infty } f_n(x) = f(x)$, contradicting the strict inequality above. Hence, we deduce $x=a$. Consequently, every subsequence of $x_n$ converges to x and therefore $x_n\rightarrow x$ as $n\rightarrow \infty$. $\square$

Applying Lemma 5 to the functions $\varPi _T$ and $\varPi _\infty$ with $x_0=\lambda$, together with Lemma 4, we obtain the result immediately. $\square$

1.3 Proof of Proposition 4

Proof

Note that $\varPi _\infty$ is a smooth function. By the first optimality condition, $\varPi _\infty '(\mu _\infty ^\star ) = 0$. We first prove that also $\varPi _T(\mu )$ is differentiable with respect to $\mu$ for all $\mu \ge 0$. Recall (2.4), which defines the cost function as a combination of the accumulated expected transient queue length, and linear staffing costs. The latter term is clearly differentiable; hence, it remains to be proved that

$$\begin{aligned} C_T(\mu ) = \frac{1}{T}\int _0^\infty \mathbb {E}[Q_\mu (t)] \, \mathrm{d}t \end{aligned}$$

admits a derivative for all $\mu \ge 0$ with T fixed. This holds if and only if $\mathbb {E}[Q_\mu (t)]$ is differentiable for all $t\ge 0$. Let $Q(0)= x\ge 0$. Following (2.2),

$$\begin{aligned} \mathbb {E}[Q_\mu (t)]= & {} \mathbb {E}[X_\mu (t)] + \mathbb {E}\left[ \max \{ x, \sup _{s\in [0,t]}\{- X_\mu (s)\} \}\right] \\= & {} (\lambda -\mu )t+ \mathbb {E}\left[ \max \{ x, \sup _{s\in [0,t]}\{- X_\mu (s)\} \}\right] , \end{aligned}$$

where the first term is differentiable. Furthermore,

$$\begin{aligned} \mathbb {E}[ \max \{ x, \sup _{s\in [0,t]} \{ - X_\mu (s) \} \} ]&= x + \int _x^\infty P(\sup _{s\in [0,t]} \{ - X_\mu (s) \} > u )\, \mathrm{d}u \\&= x+\int _x^\infty P(\hat{\tau }^0(u) \le t )\, \mathrm{d}u, \end{aligned}$$

with $\hat{\tau }^0(u)$ as defined in (3.6).

Since $-X_\mu$ is a process with no positive jumps, we may apply Corollary VII3 of [8], which states that the following equivalence between measures holds:

$$\begin{aligned} s\,P( \hat{\tau }^0(u) \in \mathrm{d}s ) \mathrm{d}u = u\,P( -X_\mu (s) \in \mathrm{d}u ) \mathrm{d}s, \end{aligned}$$

(9.8)

so that

$$\begin{aligned} \int _{x}^\infty P(\hat{\tau }^0(u) \le t )\, \mathrm{d}u&= \int _{x}^\infty \int _{0}^t P( \hat{\tau }^0(u) \in \mathrm{d}s ) \mathrm{d}u\nonumber \\&= \int _{x}^\infty \int _{0}^t\,s^{-1} u\,\,P( -X_\mu (s) \in \mathrm{d}u ) \mathrm{d}s \end{aligned}$$

(9.9)

$$\begin{aligned}&= \int _{x}^\infty \int _{0}^t\,s^{-1} u\,P( X_\mu (s) \in \mathrm{d}u ) \mathrm{d}s\nonumber \\&= \int _{0}^t s^{-1} \mathbb {E}[ \max \{x, X_\mu (s)\} ] \, \mathrm{d}s \end{aligned}$$

(9.10)

$$\begin{aligned}&= \int _{0}^t \int _{x/s}^\infty P( X_\mu (s)/s> v ) \, \mathrm{d}v \, \mathrm{d}s\nonumber \\&= \int _{0}^t \int _{x/s}^\infty P( U(\lambda s)/s > v + \mu ) \, \mathrm{d}v \, \mathrm{d}s \end{aligned}$$

(9.11)

$$\begin{aligned}&= \int _{0}^t \int _{x/s+\mu }^\infty P( U(\lambda s)/s > w) \, \mathrm{d}w \, \mathrm{d}s, \end{aligned}$$

(9.12)

where the interchange of integrals is justified by Fubini’s theorem and this last form is differentiable with respect to $\mu$. Substituting Q(0) for x straightforwardly yields differentiability of the complete cost function $\varPi _T$ for all T.

Consequently, we invoke the first optimality condition for $\mu _T^\star$ to find

$$\begin{aligned} 0=\varPi _T'(\mu _T^\star )&= \varPi _\infty '(\mu _T^\star ) + \varPsi _T'(\mu _T^\star ) + O(1/T^2)\\&= \varPi _\infty '(\mu _\infty ^\star ) + \varPsi _T'(\mu _\infty ^\star ) + (\mu _T^\star -\mu _\infty ^\star )\left[ \varPi _\infty ''(\mu _\infty ^\star ) + \varPsi _T''(\mu _\infty ^\star ) \right] \\&\qquad + \frac{1}{2}(\mu _T-\mu _\infty ^\star )^2\left[ \varPi _T'''(\xi )+\varPsi _T'''(\xi ) \right] + O(1/T^2)\\&= \varPsi _T'(\mu _\infty ^\star ) + (\mu _T^\star -\mu _\infty ^\star )\left[ \varPi _\infty ''(\mu _\infty ^\star ) + \varPsi _T''(\mu _\infty ^\star ) \right] \\&\qquad + \frac{1}{2}(\mu _T-\mu _\infty ^\star )^2\left[ \varPi '''(\xi )+\varPsi _T'''(\xi )\right] + O(1/T^2), \end{aligned}$$

for some $\xi \in [\mu _T^\star ,\mu _\infty ^\star ]$. Rearranging this gives

$$\begin{aligned} \mu _T^\star -\mu _\infty ^\star= & {} \frac{-\varPsi _T'(\mu _\infty ^\star )}{\varPi _\infty ''(\mu _\infty ^\star )+\varPsi _T''(\mu _\infty ^\star ) + \frac{1}{2}(\mu _T^\star -\mu _\infty ^\star )(\varPi _\infty '''(\mu _T^\star )+\varPsi _T'''(\xi ))} + O(1/T)\\= & {} {-}\frac{\varPsi _T'(\mu _\infty ^\star )}{\varPi _\infty ''(\mu _\infty ^\star )} \left[ 1 - \frac{\varPsi _T''(\mu )}{\varPi _\infty ''(\mu _\infty ^\star )} - \frac{1}{2}(\mu _T^\star -\mu _\infty )\frac{\varPi _\infty '''(\mu _\infty ^\star )+\varPsi _T'''(\mu _\infty ^\star )}{\varPi _\infty ''(\mu _\infty ^\star )}\right] + O(1/T)\\= & {} {-}\frac{\varPsi _T'(\mu _\infty ^\star )}{\varPi _\infty ''(\mu _\infty ^\star )} [1 + o(1)] \end{aligned}$$

for $T\rightarrow \infty$, since both $\mu _T - \mu _\infty$ and $\varPsi _T''(\mu _\infty ^\star )$ are o(1). Let

$$\begin{aligned} \mu _\bullet := \lim _{T\rightarrow \infty } \frac{T \varPsi _T'(\mu _\infty ^\star )}{\varPi _\infty ''(\mu _\infty ^\star )}. \end{aligned}$$

By (3.9) we have

$$\begin{aligned} T \varPsi '_T(\mu ) = {-} \frac{\mathbb {E}[Q(0)^2]}{2(\mu -\lambda )^2} + \frac{\lambda u_3}{3(\mu -\lambda )^3} + \frac{3\lambda ^2u_2^2}{4(\mu -\lambda )^4}. \end{aligned}$$

Together with

$$\begin{aligned} \varPi _\infty ''(\mu ) = \frac{\lambda u_2}{(\mu -\lambda )^3} \end{aligned}$$

and (4.3) we obtain the expression for $\mu _\bullet$ in (4.4). $\square$

1.4 Proof of Proposition 5

Proof

We upper bound the optimality gap by using the decomposition in (4.1).

$$\begin{aligned} |\varPi _\infty ^\star - \varPi _T^\star |= & {} \left| \hat{\varPi }_T(\mu _\infty ) + \varDelta _T(\mu _\infty ^\star ) - \hat{\varPi }_T(\mu _T^\star ) - \varDelta _T(\mu _T^\star )\right| \nonumber \\\le & {} |\hat{\varPi }_T(\mu _\infty ^\star ) - \hat{\varPi }_T(\mu _T^\star )| + |\varDelta _T(\mu _\infty ^\star )| + |\varDelta _T(\mu _T^\star )|\nonumber \\= & {} |\hat{\varPi }_T(\mu _\infty ^\star ) - \hat{\varPi }_T(\mu _T^\star )| + O(1/T^2), \end{aligned}$$

(9.13)

since $\varDelta _T(\mu ) = O(1/T^2)$ by Proposition 2. Next, we find an upper bound for $|\hat{\varPi }_T(\gamma ) - \hat{\varPi }_T(\beta )|$, with $\hat{\varPi }_T(\cdot )$ as in (4.1), in terms of the difference between $\gamma$ and $\beta$. For simplicity, denote $\hat{\gamma } = \gamma - \lambda$ and $\hat{\beta } = \beta -\lambda$, implying $\hat{\gamma }-\hat{\beta }=\gamma -\beta$. Then, using (3.9), we get

$$\begin{aligned} |\hat{\varPi }_T(\mu _\infty ^\star ) - \hat{\varPi }_T(\mu _T^\star )|&= \left| \alpha (\hat{\gamma }-\hat{\beta }) +\left( \frac{\lambda u_2}{2} + \frac{\mathbb {E}[Q(0)^2]}{2T}\right) \left( \frac{1}{\hat{\gamma }}-\frac{1}{\hat{\beta }}\right) \right. \\&\qquad \left. -\frac{\lambda ^2 u_2^2}{4T}\left( \frac{1}{\hat{\gamma }^3}-\frac{1}{\hat{\beta }^3}\right) -\frac{\lambda u_3}{6T}\left( \frac{1}{\hat{\gamma }^2} - \frac{1}{\hat{\beta }^2}\right) \right| . \end{aligned}$$

Furthermore, we have

$$\begin{aligned} \frac{1}{\hat{\gamma }}- \frac{1}{\hat{\beta }}&= -\frac{\hat{\gamma }-\hat{\beta }}{\hat{\beta }^2} + \frac{(\hat{\gamma }-\hat{\beta })^2}{\hat{\beta }^3} + O\left( (\gamma -\beta )^3\right) ,\\ \frac{1}{\hat{\gamma }^2} - \frac{1}{\hat{\beta }^2}&= -\frac{2(\hat{\gamma }-\hat{\beta })}{\hat{\beta }^3} + \frac{3(\hat{\gamma }-\hat{\beta })^2}{\hat{\beta }^4} + O\left( (\gamma -\beta )^3\right) ,\\ \frac{1}{\hat{\gamma }^3} - \frac{1}{\hat{\beta }^3}&= -\frac{3(\hat{\gamma }-\hat{\beta })}{\hat{\beta }^4} + \frac{6(\hat{\gamma }-\hat{\beta })^2}{\hat{\beta }^5} + O\left( (\gamma -\beta )^3\right) .\\ \end{aligned}$$

Substituting these yields

$$\begin{aligned} |\hat{\varPi }_T(\gamma ) - \hat{\varPi }_T(\beta )|&= \left| (\gamma -\beta )\left[ \alpha - \frac{\lambda u_2}{2 \hat{\beta }^2} + \frac{1}{2T\hat{\beta }^2}\left( \mathbb {E}[Q(0)^2] + \frac{3\lambda ^2 u_2^2}{2\hat{\beta }^2} + \frac{2\lambda u_3}{3 \hat{\beta }}\right) \right] \right. \\&\quad \left. - (\gamma -\beta )^2\left[ \frac{\lambda u_2}{2 \hat{\beta }^3} + \frac{1}{2T\hat{\beta }^3}\left( \mathbb {E}[Q(0)^2] - \frac{3\lambda ^2 u_2^2}{\hat{\beta }^2} - \frac{\lambda u_3}{\hat{\beta }}\right) \right] \right| \\&\quad + O\left( (\gamma -\beta )^3\right) . \end{aligned}$$

Given that $\mu _T^\star = \mu _\infty ^\star + \mu _\bullet /T + o(1/T)$, we find

$$\begin{aligned} |\hat{\varPi }_T(\mu _\infty ^\star ) - \hat{\varPi }_T(\mu _T^\star )|&= \frac{|\mu _\bullet |}{T}\left( \alpha - \frac{\lambda u_2}{2(\mu _\infty ^\star -\lambda )^2}\right) + O(1/T^2)\\&= \frac{|\mu _\bullet |}{T}\left( \alpha - \frac{\lambda u_2}{2(\sqrt{\lambda u_2/2\alpha })^2}\right) + O(1/T^2) = O(1/T^2), \end{aligned}$$

which concludes the proof. $\square$