Second-order bounds for the M/M/ queue with random arrival rate
Introduction
The majority of stochastic models in queueing theory assume known arrival processes, which facilitates exact analysis. In fact, the dominant assumption is that potential customers arrive according to a Poisson process with a known intensity . In contrast, we interpret the arrival rate as an unknown parameter of which partial information is available. The parameter is then replaced by a random variable with some distribution representing what is believed or known about the market size. We primarily focus on the situation where we know the mean and variance of , with the connection to classical queueing theory when the variance is zero and, hence, the random variable is known to be always.
We set out to solve an optimization problem for determining tight bounds for the expected waiting time for all distributions of that comply with the partial information. For mean-variance information, we will show that these tight bounds are attained by two-point distributions, essentially saying that the worst-case market is one where the actual market size attains its maximum size or is well below the expected market size. To establish such tight bounds, we invoke primal-dual techniques for solving semi-infinite linear programs. Successfully applying these techniques will prove to depend crucially on the properties of the expected waiting time viewed as function of the arrival rate . For the main queueing system in this paper, the M/M/ queue, we will leverage that the expected waiting time is (i) increasing and (ii) convex in , with (iii) the first derivative with respect to (w.r.t) being a convex function of . All three properties will turn out to be crucial for solving the optimization problems in this paper, and the convex derivative will be particularly useful for the setting with mean-variance information.
We summarize the main contributions in this paper as follows. For the M/M/ queue with partially known arrival rate, we establish novel tight bounds. When mean-variance-support information is known, we show that the tight bound for the expected waiting time is attained by a two-point distribution. Our proof of these tight bounds combines two results from different areas. The first result stems from the area of optimization and says the semi-infinite linear program that describes the tight bound for the expectation of a function with a convex derivative is attained by a two-point distribution. The second result stems from queueing theory and shows that the expected waiting time as a function of the arrival rate is convex, and we establish that its derivative is also a convex function of the arrival rate. Combining these two results proves to be an effective strategy for finding distributionally robust bounds for queueing systems subject to parametric uncertainty. The tight bounds are leveraged for analysis of unobservable M/M/ queues that cater to rational users who decide to join or balk based on expected utility. We use the wait bounds to bound the equilibrium arrival rate, which in turn, leads to tractable maximin analyses for setting the revenue-maximizing price. In this way, we extend the classical strategic queueing literature with known model parameters as in [36, 13] to a setting where the arrival rate is only partially characterized in terms of range, mean, and variance.
The paper is structured as follows. In Section 1.1, we present a brief overview of related literature. Section 2 contains the proof of our main result, which provides bounds for the expected waiting time under uncertain arrival rates. Section 3 applies these bounds to a rational queueing model. Finally, in Section 4, we conclude and propose some further research directions.
Related literature
This research connects four themes in the queueing literature: random arrival rates, second-order bounds, parametric convexity, and rational queueing.
Queues with random arrival rate. Due to its mathematical tractability, the M/M/ queue is among the cornerstones of queueing theory and is widely applied in service operations management. To enhance practical relevance, several studies proposed to relax the assumption of Poisson arrivals, and instead model the arrival rate as a random variable [1, 23, 54, 47, 27, 5, 55, 4]. This gives rise to a Poisson mixture model for the arrival process, which can deal with forecast errors, overdispersion (compared with the natural fluctuations of a Poisson process), and unknown market size. The mathematical formulae of the M/M/ model, such as the Erlang-C formula for the delay probability, then still apply once the mixing distribution can be characterized. Jongbloed and Koole, [29] showed how to estimate this mixing distribution based on available data of arrivals.
Second-order bounds. For this M/M/ queue with random arrival rate, we derive performance bounds that do not require the full mixture distribution, but only utilize the first two moments of the arrival rate. To derive such second-order bounds, we need to solve a semi-infinite linear program. To do so, we first show that the expected stationary waiting time can be written as the expectation of a function of the random arrival rate that has a convex derivative. We then establish a general result, which provides tight second-order bounds as the solutions to semi-infinite linear programs for the expectation of functions with convex derivatives. Such semi-infinite linear programs arise naturally in studying second-order bounds, or more generally, moment problems. For example, [2, 10, 11] use second-moment information to compute tight bounds for distributionally robust stochastic programming. Second-order bounds also play a predominant role in the theory on performance bounds for queues, which started with the work of Kingman, [32]. Kingman derived bounds for the mean waiting time in the GI/G/1 queue, expressed in terms of the first two moments of the interarrival- and service-time distributions; however, these bounds are not in general tight.
Convexity and beyond. Numerous studies have taken advantage of the insight that moments of the waiting time can be viewed as the expected value of a convex function. Vasicek, [49] showed that the variance of waiting time under a general queueing discipline does not exceed that under the LCFS discipline, and Hajek, [18] established that among all arrival processes for an exponential server queue with specified arrival and service rates, the arrival process which minimizes the expected waiting time is the process with constant interarrival times. Weber, [51] revealed that the mean queueing time in the G/GI/ queue is nonincreasing and convex in the number of servers, , indicating that marginal analysis is optimal for determining server allocation among various service facilities so as to minimize the total expected queueing time. The proofs in [49, 18, 51] rely on convexity properties. Convexity of queueing performance metrics is a highly desired property in the parametric optimization of stochastic systems. A substantial body of research employs direct algebraic methods to establish convexity of the steady-state performance metrics as functions of the design parameters; see, for example, [15, 34, 28, 22, 19, 20, 52]. For a comprehensive discussion on the related notion of stochastic convexity, we refer the interested reader to [42, 44, 43]. However, as mentioned earlier, for the present study we require convexity of the objective function’s derivative. Several closely related concepts emerge in the extremal analysis of queueing systems. In this context, the need for a convex derivative is replaced by closely related sufficient conditions that warrant the application of the theory of complete Tchebychev systems [30]. For a detailed discussion of the relevance of these systems to queues, please refer to [12, 7, 17].
Rational queueing with limited information. A mature literature exists [25] that seeks to elucidate the rational decision-making processes of users who decide whether or not to access delay-sensitive services based on utility. Naor, [36] initiated this line of work for the observable M/M/1 queue with a customer utility that depends linearly on the price, value of service, and the expected waiting time. Edelson and Hilderbrand, [13] investigated similar questions for the unobservable M/M/1 queue. These works explore the environment with a deterministic arrival rate. Liu and Hasenbein, [35] expanded upon Naor’s model by assuming that the arrival rate is drawn from a known probability distribution accessible to the decision maker. This concept was further extended to the unobservable setting in [6]. Their research demonstrated that a socially optimal pricing strategy results in a lower expected arrival rate compared to a price set by a revenue-maximizing decision maker. Hassin et al., [26] examined the unobservable model with a random arrival rate from a distinct angle, where strategic customers base their decisions on a rate-biased time-average property when estimating their expected waiting time. Wang et al., [50] extended the traditional observable model to a distributionally robust setting by taking into account an uncertain arrival rate governed by an unknown underlying probability distribution.
Tight bounds for expected waiting time with limited market knowledge
In this section, we derive bounds for the expected waiting time with an uncertain arrival rate. Section 2.1 introduces the basic problem of finding a tight bound for the expected waiting time. Section 2.2 states the second-order bounds for general convex functions with a convex derivative. In Section 2.3, we apply these second-order bounds to the expected waiting time in the M/M/ queue. Section 2.4 incorporates unimodality information to refine the second-order bounds. In Section 2.5, we explore the extension wherein the service rate parameter is likewise a random parameter. Further, as an example of the broader applicability of our approach, we obtain second-order bounds for the expected waiting time in the M/G/1 queue with random arrival rate.
The model
Consider an M/M/ queue with Poisson arrivals with rate , unit mean exponential service requirements, and parallel servers, each operating with service rate . Let denote the expected waiting time experienced by a customer in this queueing system (including the time in service), which depends on the number of servers and the arrival rate . For notational convenience, we omit at first the functional dependence on the service rate parameter , but our findings remain valid for any given value of this parameter. Suppose that . Then the single-server M/M/1 queue has as expected waiting time . Observe that
and
Hence, is an increasing, convex function of with a convex derivative. We now present some general results for functions that have the same properties as . We present the results in general form, because we will apply these results later to with as well, to obtain the solutions for
(1) |
in which we focus on the case with the latter being shorthand notation for , the set of all distributions with mean , standard deviation and the support contained in the interval . That is,
where comprises all probability distributions with support contained in . We shall refer to such a set as an ambiguity set. Henceforward, the market size is viewed as a random variable with a distribution contained in . We further assume from this point onward that the system under consideration is stable for each realization of the random arrival rate; that is, .
Bounding the expectation of a function with convex derivative
We first develop a general bound for with being a function with a convex derivative, just like the function . We therefore solve where the constraints that define correspond to known first and second moments and support contained in . We thus need to solve
(2) | ||||
s.t. | ||||
a semi-infinite linear program with an (infinite-dimensional) decision variable, the probability distribution , which must reside within . Since there are three constraints (including ), it is known from general optimization theory (see, e.g., [40, Theorem 1], [45, Lemma 3.1], [46]) that among the optimal solutions of (2) there is a distribution with at most three points of support. We next prove a stronger result which says that among the optimal solutions to (2), there is a distribution with two points. The proof of this result uses the dual problem of (1), defined as finding the vector of dual variables that solves
(3) | ||||
s.t. |
and uses that the objective function has a convex derivative. Denote the optimal dual solution as .
Lemma 1.
Suppose that and . The probability distribution that solves (2) is supported on the points at which the function coincides with .
Proof.
According to Proposition 3.4 in [45], problems (2) and (3) have the same optimal value when the moment vector is an interior point of the set of feasible moment vectors, which is implied by the stated conditions and . Furthermore, [45, Proposition 3.4] asserts that, in this case, the optimal dual solution is attained if the common optimal value is finite. Since is continuous on the compact support , it readily follows that the optimal value is indeed bounded. Hence, the dual optimal solution is attained. Further, [45, Corollary 3.1] says that if the support is compact, the moment functions and the objective function are continuous, and the primal problem is consistent (i.e., the ambiguity set is nonempty), then there is no duality gap and the optimal primal solution is attained. Under the stated conditions on the moment vector , the primal problem is indeed guaranteed to admit a feasible distribution. Thus, it immediately follows from the continuity of and that the primal optimal solution is also attained. Therefore, the optimal solutions of both (2) and (3) are attained, and there is no duality gap. This, in turn, implies that complementary slackness holds [45, Proposition 2.1] and completes the proof. ∎
Proposition 1 (Extremal two-point distributions).
Let be a continuously differentiable function with its derivative being either convex or concave on the bounded interval . Suppose that and . Then for a random variable with distribution , the tight upper bound of is attained by a discrete distribution with two support points.
Proof.
It is evident that the optimal distribution cannot be a one-point distribution, since , so that is bound to intersect with more than once by Lemma 1. We first consider the case with being either strictly convex or strictly concave. It suffices to demonstrate that a candidate dual function coincides with at precisely two points, rather than three. Figure 1 illustrates the functions and , and their derivatives.


We proceed by contradiction. Suppose that there exists a feasible such that it coincides with at three points. Let the first intersection point of the derivatives in Figure 1(b) correspond to a tangent point for the original functions. Due to dual feasibility, has to be greater than directly after the tangent point, and can only be tangent to once since is linear and is strictly convex. Hence, again by dual feasibility, there cannot exist a second intersection point in the interior of . Therefore, a candidate dual solution that coincides with thrice should intersect once in the interval and at both boundary points. Let be such that it is tangent to at some and intersects at . Assume now that also touches at . This implies that for to be dual feasible. However, as shown in Figure 1, this cannot hold, for ′ is strictly convex. Hence, we arrive at a contradiction. Likewise, for strictly concave , we construct such that it touches at , and it is easily shown that imposing this function to intersect also at leads to a similar contradiction.
To complete the proof, we must still consider the case where is either convex or concave, but not strictly so. This implies that is linear over certain intervals of its domain. Suppose we are provided with a dual feasible function . If does not align with ′ on any of these intervals, the preceding argument remains valid. However, if coincides with ′ on such an interval, and must match over the entire interval , for some , by dual feasibility. This condition results in a continuum of potential choices for the support of the extremal distribution, as dictated by complementary slackness. Nevertheless, this subsumes the situation of a two-point support. Therefore, the assertion presented in the claim holds true. An analogous argument applies for a concave , completing the proof of the claim. ∎
Proposition 1 partly overlaps with Theorem 5.1 in [2]. Proposition 1 stretches the conditions in [2, Theorem 5.1], as [2, Theorem 5.1] requires to be convex. On the other hand, [2, Theorem 5.1] allows the derivative to be convex on the interval and concave on for some . For the purpose of this paper, [2, Theorem 5.1] would have sufficed, but we present Proposition 1 as a slight generalization and to provide a self-contained exposition.
Proposition 1 defines a class of objective functions for which the mean-variance information set leads to two-point worst-case distributions. The defining feature of this class of functions is the convex derivative. Based on general optimization theory, we knew already that three points suffice for describing the worst-case distribution, but with Proposition 1 this is reduced to two points. It is thus enough to optimize over two-point distributions in the search for the extremal distribution that attains the tight bounds for . Kreĭn and Nudel’man, [33, Chapter 4] demonstrate that, when is continuously differentiable with a strictly convex derivative on , a unique two-point distribution solves the moment problem (2). The following lemma establishes a similar result.
Lemma 2.
Consider a function that is continuously differentiable on the bounded interval , and let be a random variable with distribution .
(i) If the derivative is convex on , the tight upper bound for is attained by a two-point distribution with values
and corresponding probabilities
.
(ii) If the derivative is concave on , the tight upper bound for is attained by a two-point distribution with values
with corresponding probabilities
.
Proof.
A two-point distribution with mean and variance 2 has support
(4) |
with probabilities , respectively. We thus need to solve with
and
due to the support being contained in . To show that is a maximizer, we will prove that is nonincreasing in . For differentiable , we have that
(5) |
In order for to be nonincreasing, we need
or
This inequality holds by convexity of ′. Thus, as is nonincreasing in , the maximum occurs at the lower bound of the feasible set (i.e., ). Substituting the optimal value ∗ into our parameterized expression for the two-point distribution then yields assertion (i).
For assertion (ii), note that for concave ′, is nondecreasing in and therefore maximized by . ∎
Notice that assertion (ii) yields the extremal distribution that solves , since determining the lower bound is tantamount to maximizing . Returning to the expected waiting time in the M/M/1 queue, we see that the tight bounds for are attained by the two-point distributions stated in Lemma 2. We shall prove this result for the multi-server system in the next subsection.
Let us mention two other streams of literature that also deal with tight bounds for stochastic models, and in particular three- and two-point worst-case distributions. The first stream stems from statistics. [16] presents general results for the maximization of integrals over distributions of which only the first two moments are given, and where for convex integrands the maximizing distribution turns out to be a three-point distribution. Among the many statistical applications of the results in [16] are bounding the expected sample maximum and the expected sample range, where for the latter, [24] claims that the extremal distribution can be reduced to a two-point distribution. The other line of research in which three- and two-point extremal distributions feature prominently deals with establishing tight bounds for the waiting time moments in the general GI/G/1 queue, or for more specific queues such as M/G/1 and G/M/1 queues; see [53, 12, 41, 37, 48] among others. In particular, a long-standing conjecture for the expected waiting time in GI/G/1 queue says that the tight bounds follow from the two-point distributions stated in Lemma 2. While recently some progress has been made, this conjecture is still open [9, 8].
Bounding the expected waiting time as function of the arrival rate
To apply the results in the previous section to the multi-server M/M/ queue with random arrival rate, we must show that the expected waiting time, as a function of the arrival rate, has a convex derivative. Therefore, we need to study the third derivative of with respect to . The expected waiting time is a classic formula in the queueing literature, for which a host of properties have been discovered, including convexity with respect to ; see, e.g., [15, 34]. But to the best of our knowledge, the property of a convex derivative with respect to has not been reported. However, this property was reported in [39] for the expected queue size defined as
(6) |
where is the Erlang-C formula, the probability that an arriving customer experiences delay in the M/M/ queue, and . Randhawa, [39] proved that as function of the arrival rate , indeed has a convex derivative. By Little’s law, , and one could try to see whether the convex derivative of implies the same property for . Unfortunately, we do not see a direct proof for this argument, in which we can use convexity of the derivative of the expected queue length to demonstrate convexity of the derivative of the expected waiting time . We therefore choose to provide a standalone proof for the fact that has a convex derivative, just like . The proof method we apply is largely inspired by that in [39]. In the following, we provide a proof sketch. The details are relegated to the Appendix. First, we determine an expression for the third derivative of the expected waiting time as a function of the server utilization . Then, we write this expression as a function of the Erlang-C formula and the number of servers . It then suffices to show that this expression is nonnegative for all values of and , which implies convexity of as a function of with the number of servers fixed.
Lemma 3 (Derivative w.r.t. is convex).
The expected waiting time as function of has a convex derivative.
Convexity of as a function of implies that the extremal distribution of is given by the extremal distribution in Lemma 2(i), thus yielding our main result.
Theorem 1 (Second-order bound for random arrival rate).
Consider an M/M/ queue with random arrival rate that follows a distribution belonging to the ambiguity set . The tight upper bound for the expected waiting time is attained by the two-point distribution with support
and corresponding probabilities
We demonstrate the effectiveness of our second-order bounds through a numerical experiment. The results are presented in Table 1. We assume that the “true” distribution of the random arrival rate follows a beta distribution with support and standard deviation fixed to . We consider the scenario with servers and vary the mean to obtain results for multiple utilization levels. The average utilization level is denoted by . As shown in Table 1, the second-order bounds proposed in Theorem 1 perform well for low- and medium-utilization regimes. However, for high utilization, the bounds deviate significantly from the true value when assuming a beta distribution for . Furthermore, the bounds widen as the upper bound of the support increases. The second column corresponds to the second-order lower bound, which can be obtained using the extremal distribution in Lemma 2(ii). Notice that this bound depends on , the lower bound of the support, rather than .
LB | UB | UB | UB | ||||
---|---|---|---|---|---|---|---|
0.2 | 1.0445 | 1.0449 | 1.0940 | 1.0449 | 1.1199 | 1.0449 | 1.4312 |
0.5 | 1.3389 | 1.3440 | 1.4655 | 1.3440 | 1.5315 | 1.3440 | 2.3236 |
0.7 | 1.9753 | 2.0094 | 2.3200 | 2.0096 | 2.5002 | 2.0098 | 4.6625 |
0.8 | 2.8099 | 2.9442 | 3.5518 | 2.9465 | 3.9446 | 2.9497 | 8.6515 |
0.9 | 5.3921 | 6.4285 | 7.6443 | 6.6019 | 9.0133 | 6.9387 | 25.0550 |
Sharper bounds for unimodal market size
Often more is known about the market-size distribution than just the mean and variance. For example, we may have some structural information regarding the distribution of the market size. In this section, we show how to use this structural information to improve the tight bounds, corresponding to the situation in which the service provider has the additional information that the market-size distribution is unimodal, which can be understood as there being one primary market segment that constitutes the bulk of the arrivals requesting service.
A random variable is said to be unimodal with mode if its distribution can be characterized as the mixture of a Dirac distribution and a distribution function that is a concave function on and a convex function on . Some examples of unimodal distributions are the normal, exponential, and beta probability distributions. In order to incorporate such unimodality, we apply Khinchine’s characterization theorem, which states that a random variable has a unimodal distribution with mode zero if, and only if, there exists a random variable that can be expressed as , where is a uniformly distributed random variable on independent of some random variable (see, e.g., [14, p. 158]). We adopt the following approach from [31] (see also [3]). First, we transform the problem from the unimodal random variable to the auxiliary variable . Then, we establish tight bounds for and use them to obtain bounds for . To be specific, if is unimodal with mode , then is also unimodal but with mode 0. Hence, according to Khinchine’s theorem, we have , where is uniformly distributed over . According to Kemperman’s approach, we can transform the moment problem in terms of the random variable into one in terms of the auxiliary random variable using that for any function , , where
(7) |
By taking , we can relate the moments of to those of through
so that . The mean and variance of , in terms of the mean and variance of , are given by
We can use the same integral relationship to find tight bounds for when is unimodal with mode . Specifically, this problem is equivalent to finding bounds for subject to the moment constraints on , where
We already have the solution for this transformed problem, since the (conditional) expectation operator (7) preserves the curvature properties. By substituting the transformed moments into the result stated in Lemma 2, we obtain the following result.
Theorem 2 (M/M/ queue with unimodality).
Consider an M/M/ queue with random arrival rate . Suppose that is unimodal with mode , mean , variance 2, and has support contained in . Let
Then, the tight upper bound for the expected waiting time is given by
where
Table 2 demonstrates the performance of the tight mean-variance bounds that incorporate unimodality information. Evidently, incorporating this structural information significantly sharpens the bounds compared to the standard mean-variance bounds without this information.
Var-UB | Unimod-UB | Var-UB | Unimod-UB | |||
---|---|---|---|---|---|---|
0.2 | 1.0449 | 1.1199 | 1.0583 | 1.0449 | 1.4312 | 1.0686 |
0.5 | 1.3440 | 1.5315 | 1.3879 | 1.3440 | 2.3236 | 1.4277 |
0.7 | 2.0096 | 2.5002 | 2.1502 | 2.0098 | 4.6625 | 2.3124 |
0.8 | 2.9465 | 3.9446 | 3.2477 | 2.9497 | 8.6515 | 3.7098 |
0.9 | 6.6019 | 9.0133 | 7.1468 | 6.9387 | 25.0550 | 9.3105 |
Broader robust applications
In practice, it might turn out that other parameters besides the arrival rate are uncertain, with only partial information available. It seems reasonable to assume that decision makers may not always have precise knowledge of service rates, as these rates might depend on specific server characteristics. Hence, it is of interest to establish also second-order bounds for the service rate at which the different servers operate. In this regard, we consider the scenario in which the service rate is replaced by a random variable . Convexity of the mean queue size of the M/M/ queue with respect to the system utilization is a classic result [15, 34]. The following result pertains to the multi-server setting, where we again will rely on the results in [39].
Lemma 4 (Derivative w.r.t. is concave).
The expected waiting time as function of has a concave derivative.
Proof.
Denote . To prove that is concave, we will use the result in [39] which proves convexity of as function of . Interpreting indeed as function of , define and , where is the variable and the parameters and are fixed to arbitrary constants. To see that the composite function is decreasingly convex, notice that
Little’s law indicates that , with a fixed constant, which completes the proof. ∎
Hence, is a decreasing, convex function of with a concave derivative. By Lemma 2(ii), the extremal distribution is given by a two-point distribution with support and respective probabilities
That is, the worst case is given by servers which either operate at their lowest service capacity or just above the mean service rate.
Moreover, it is possible to work with multiple uncertain parameters simultaneously. We next consider
which turns out to be solvable in closed form:
Theorem 3 (M/M/ queue with partially known service rate and market size).
Let the service rate follow a marginal distribution and the market size follow a marginal distribution . Suppose that . If and are independent, the sharpest possible upper bound for is attained by the product measure , with and as defined in Lemma 2(i) and (ii), respectively.
Proof.
The tower rule yields
where the expectation is well defined due to . Hence, fixing the distribution of gives the moment problem
For a fixed constant, has the required convexity properties, so is the solution to the moment problem. For to work also for , we need to check the signs of the first three derivatives. We interchange differentiation and expectation operations (which is allowed as is bounded), so that
As the expectation operator is nonnegative,
and likewise,
which proves that maximizes for any distribution of . Substituting the two-point distribution gives the moment problem
Based on the same line of reasoning, this expression is maximized by . ∎
Theorem 3 leverages that the extremal distributions for univariate settings are independent of the precise form of , and can thus be applied recursively when considering multiple uncertain parameters.
The applications of the second-order bound extend beyond the M/M/ queue and potentially include a vast range of uses. However, to apply the bounds, we need the “third-order property” of a convex derivative. As a concrete example, consider the Pollaczek-Khinchine formula for the expected waiting time in the M/G/1 queue given by
where denotes the variance of the service times. Now observe that
and moreover,
Hence, is increasingly convex in , and can be bounded using Lemma 2, yielding the following result.
Proposition 2 (M/G/1 queue with variance information).
Let the market size have a distribution that resides in the ambiguity set , and let denote the expected waiting time in the M/G/1 queue. Then, the tight upper bound for is attained by the two-point distribution stated in Theorem 1.
Hence, despite some potentially tedious algebra, the second-order bounds that incorporate variance information have the potential to be widely applicable to other stochastic systems. Notwithstanding, as a cornerstone of queueing theory, the M/M/ queue already boasts numerous pertinent applications. These include, but are not limited to, staffing optimization problems and rational queueing models, both of which find extensive usage in service operations management. We shall delve into the topic of rational queueing in the next section.
Rational queueing model
We now apply the results derived in the previous section to rational queueing models. Section 3.1 introduces the rational queueing model with a known market size. In Section 3.2, we extend this model to the setting in which customers are ambiguity-averse about the total market size. We provide a numerical example in Section 3.3.
Model with known market size
Consider a firm that sells delay-prone services to a market of rational delay-sensitive individuals. Individuals value service, but dislike waiting, and will only join when the net service value exceeds the wait costs. They cannot observe real-time queues and instead estimate their expected waiting time costs based on beliefs or information about the total arrival rate of all potential customers. The arrival rate or market size is measured in terms of the scale parameter , so that the expected time between arrivals of consecutive individuals is . All individual joining decisions together give an equilibrium arrival rate or market share e, which could be viewed as the product of times the probability that an individual decides to join. We then consider the firm as a price-setting monopolist seeking to maximize revenue. With the price for service, the firm can influence the net service value of individuals, and hence the joining probability. A low price yields small margins but high joining probability, while a high price increases the profit per customer but suppresses the market share. The firm should thus strike the optimal balance between profit per customer and market share. Such challenges have been thoroughly addressed in the rational queueing literature [25].
Assume that the firm operates according to an M/M/ queue, but that the arrival process consists of individuals that decide to join-or-not based on a rational-choice model. An individual’s utility in acquiring service is defined as
with the value of service, the price for service, the wait costs per time unit, and e the effective arrival rate of joining individuals. Here we assume that individuals cannot observe real-time queue lengths, but instead know how the expected waiting time varies with overall demand. Hence, the rational choices are based on a linear relation between net service value and expected waiting time costs, and an individual will only join when . This is a common set-up in the rational queueing literature with self-interested individuals that tend to overcrowd systems, because they ignore externalities—inconvenient side effects—that their decisions have on others. When an individual decides to join, this will increase congestion and wait for all customers. Standard game theory then predicts that all individual decisions together will lead to an equilibrium, expressed in the equilibrium arrival rate, the probability that an arbitrary individual joins the systems multiplied by the market size . Since not joining the queue gives zero utility, the equilibrium joining probability then solves . The arrival process of those who decide to join is then a Poisson process with rate . To optimize revenue, the firm needs to set the price that strikes the best balance between margin and market share . Notice that the joining probability is an implicit function that decreases with and can only be determined as the solution to when individuals know the arrival rate parameter . The firm then effectively solves the monopoly pricing problem, to attain maximal revenue. The optimal solution to this problem for a given market size is already well established [25]. Instead, in the next subsection, we choose to view the market size as a random variable of which only limited information is available.
Model with ambiguous market size
We interpret the market size as an unknown parameter about which the individuals that constitute the market form beliefs. The parameter is then replaced by a random variable with a distribution representing what individuals believe or know about the market size. We primarily focus on the situation where individuals only know the support, mean, and variance of . From the individual’s perspective, the variance of expresses uncertainty about the market size, and hence uncertainty about the expected waiting time. For this setting with partially-informed customers, or an ambiguously specified arrival process, we consider the revenue-maximizing firm as a maximin decision maker, first determining the worst-case market (a minimization problem), and then maximizing the revenue by selecting the best price for this worst-case market. The firm thus attempts to solve the maximin problem
(8) |
with the distribution of the market size and the ambiguity set that contains all distributions that satisfy the partial information, given by the mean, variance, and support. Minimax and maximin optimization problems such as (8) arise naturally in decision making under uncertainty. Our strategy for solving such problems will be to solve first the minimization problem , and then the maximization problem for this worst-case market. The technical challenge is to solve the minimization problem, which requires determining tight bounds for the expected waiting time for all distributions of that comply with the partial information. However, this part was already resolved in the previous section, where we derived an upper bound on the expected waiting time. A tight upper bound on the expected waiting time presents the worst-possible scenario in the rational choice model and yields the corresponding tight lower bound on the market share. Observe from the expression for the utility that solving is tantamount to solving . That is, the worst-possible expected market share arises from the worst-possible utility and hence worst-possible expected waiting time. To make this more precise, notice that the maximin problem (8) can also be written as an optimization problem in terms of the equilibrium joining strategy ; that is,
(9) | ||||
where the equivalence follows from the fact that the expected value is a known constant as it is part of the information contained in the ambiguity set. To see that (9) is equivalent to (8), we carefully go over the decision processes of the firm and the customers. First, the firm sets a price , and then the customers settle on a joining strategy with expected utility , while considering all possible distributions that possibly govern the market size . That is, the customers solve the equation , with the expected value operator appearing here because of the customers’ internal calculations. In equilibrium, the utility equation can then be rewritten as . Now, if the firm sets the price to , the expected revenue will equal this price times the market share, , yielding (9). In order to solve the utility equation, we must make additional assumptions regarding how customers experience waiting times in a mixed Poisson model.
The next question then is a rather philosophical one: What to assume for the wait expected by an arbitrary arriving customer? One natural choice is ; see, e.g., [5, 6]. Here, we assume nature picks a realization , and customers from time 0 onward arrive according to a Poisson process with this rate. Customers remain unaware of the specific universe (i.e., event ) they live in but hold a (Bayesian) prior belief about the probability of residing in one universe compared to another. Another option is to base the rational decision made upon arrival on the posterior distribution of , the updated version of the prior distribution of that accounts for size bias. That is, the customer uses the conditional distribution of given the realization of its own arrival event. Hassin et al., [26] call this phenomenon RASTA (Rate-biased ASTA), as a counterpart of PASTA. For any function , the posterior expectation of at arrival instants is , and when is nondecreasing and convex, then . As an application, consider the M/M/ queue. Then, the expected waiting time with a random arrival rate equals
The utility of a customer who joins the queue when all the other customers use strategy equals
and the best response of an individual customer is to join if and only if . When accounting for size bias, we thus need to solve
(10) |
to find the worst-case market share. Alternatively, assuming (P)ASTA instead of RASTA, we have to consider
(11) |
It is easy to show that (10) is also solved by the distribution in Lemma 2(ii), as it is known that is increasingly convex in (see [39]) and is a known constant since the mean is contained in the information set. For the RASTA arrival assumption, we assume the market consists of individuals that base their decision on the posterior distribution of . Then, the worst-case joining probability solves the equation
with , and as defined in Theorem 1. For PASTA arrivals, it is assumed the market consists of individuals that base their decision on the prior distribution of . Then, the worst-case joining probability is found by solving in the equation
with , and again as defined in Theorem 1. For conciseness, let denote this extremal distribution. We next determine the optimal price for PASTA arrivals. To accomplish this, we start by deriving an expression for the derivative of the expected waiting time with respect to , using the equations for and from [15]. By substituting into the expression for and applying Little’s law, the product rule, and the chain rule, we can determine
Proposition 3 (Optimal joining probability, PASTA).
Consider the M/M/ queue with PASTA arrivals. Suppose that . Then the joining probability that maximizes the revenue for the firm, in the worst-case market , is the unique solution in of
(12) |
Proof.
We shall show that the solution to (12) is the joining probability that maximizes
It suffices to consider as, due to the assumption, both and are not equilibrium strategies. It is fairly straightforward to see that the minimization problem is solved by , since the expected waiting time has the required third-order property. We next solve the maximization part. Taking the derivative with respect to , we obtain the first-order condition
which can be written as (12). Notice that we are allowed to interchange the differentiation and expectation operations as is bounded. Since is an increasing function of , as is (strictly) convex, this implies that (12) has a unique solution that is the global maximizer, as the objective function is a concave function of . ∎
We next determine the minimax price that induces the optimal joining strategy .
Proposition 4 (Optimal maximin price, PASTA).
For the M/M/ queue with PASTA arrivals, the optimal price for the firm that solves the maximin problem (8) is
(13) |
Numerical example
We next present numerical results for the M/M/2 model with random market size , , and PASTA arrivals. Figure 2 illustrates the equilibrium joining probabilities and revenues for various price values and different levels of variance 2. As expected, both the joining probabilities and revenues decrease with an increase in the dispersion of the market size. This observation aligns with the intuition that greater variance indicates more uncertainty about the market size and, as a result, greater uncertainty about the expected waiting time. Another noteworthy characteristic of the ambiguous model, as compared to the model with a known arrival rate, is that the equilibrium joining probability cannot exceed a certain level since, otherwise, the M/M/ system can become unstable.


In conclusion, we have shown that our second-order bounds can aid in deriving the optimal maximin price for a firm seeking to maximize revenue. This firm serves a rational customer base that makes decisions based on limited information about the overall market size. Interestingly, this adverse market follows a computationally tractable two-point distribution, thus allowing for the application of established techniques designed for rational queueing models with stochastic arrival rates.
Conclusions and further research
We have presented tight second-order bounds for the expected waiting time in an M/M/ queueing system with an uncertain arrival rate. We have derived these bounds by using the fact that the expected waiting time can be written as the expectation of a convex function with a convex derivative. This convex derivative, in conjunction with mean-variance information, gave rise to a specific semi-infinite linear program with an explicit solution in terms of the two-point worst-case distribution.
Looking forward, it is worthwhile to search for other information sets that have, just as mean-variance, explicit and perhaps even two-point worst-case distributions. In Appendix B, we consider semi-variance and show that the worst-case distribution is either a three-point or a two-point distribution with some support points that can only be determined numerically. Hence, while variance leads to a closed-form two-point distribution, semi-variance leads to a less explicit worst-case scenario. The proof in Appendix B is based on solving the semi-infinite linear program with a primal-dual argument and again uses the convex derivative. Several recent papers have considered mean absolute deviation instead of variance, which leads for convex functions to an extremal three-point distribution that, as with variance, no longer depends on the objective function’s precise form; see, for example, [38, 48].
One might also consider data-driven information sets, such as the Wasserstein ambiguity set, which are equipped to handle limited data, as estimates of means and dispersion measures might not always be statistically accurate enough. Additionally, it is insightful to investigate to which other stochastic systems analogous types of distributionally robust techniques to control for parameter uncertainty can be applied, a topic we briefly touched upon in Section 2.5.
References
- Avramidis et al., [2004] Avramidis, A. N., Deslauriers, A., and L’Ecuyer, P. (2004). Modeling daily arrivals to a telephone call center. Management Science, 50(7):896–908.
- Birge and Dulá, [1991] Birge, J. R. and Dulá, J. H. (1991). Bounding separable recourse functions with limited distribution information. Annals of Operations Research, 30(1):277–298.
- Brockett and Cox Jr, [1985] Brockett, P. L. and Cox Jr, S. H. (1985). Insurance calculations using incomplete information. Scandinavian Actuarial Journal, 1985(2):94–108.
- Chen and Henderson, [2001] Chen, B. P. and Henderson, S. G. (2001). Two issues in setting call centre staffing levels. Annals of Operations Research, 108:175–192.
- Chen and Hasenbein, [2017] Chen, Y. and Hasenbein, J. J. (2017). Staffing large-scale service systems with distributional uncertainty. Queueing Systems, 87(1):55–79.
- Chen and Hasenbein, [2020] Chen, Y. and Hasenbein, J. J. (2020). Knowledge, congestion, and economics: Parameter uncertainty in Naor’s model. Queueing Systems, 96(1):83–99.
- Chen and Whitt, [2020] Chen, Y. and Whitt, W. (2020). Extremal models for the GI/GI/ waiting-time tail-probability decay rate. Operations Research Letters, 48(6):770–776.
- Chen and Whitt, [2021] Chen, Y. and Whitt, W. (2021). Extremal GI/GI/1 queues given two moments: Exploiting Tchebycheff systems. Queueing Systems, 97(1):101–124.
- Chen and Whitt, [2022] Chen, Y. and Whitt, W. (2022). Correction to: Extremal GI/GI/1 queues given two moments: Exploiting Tchebycheff systems. Queueing Systems, 102:553–556.
- Delage and Ye, [2010] Delage, E. and Ye, Y. (2010). Distributionally robust optimization under moment uncertainty with application to data-driven problems. Operations Research, 58(3):595–612.
- Dokov and Morton, [2005] Dokov, S. P. and Morton, D. P. (2005). Second-order lower bounds on the expectation of a convex function. Mathematics of Operations Research, 30(3):662–677.
- Eckberg Jr, [1977] Eckberg Jr, A. E. (1977). Sharp bounds on Laplace-Stieltjes transforms, with applications to various queueing problems. Mathematics of Operations Research, 2(2):135–142.
- Edelson and Hilderbrand, [1975] Edelson, N. and Hilderbrand, D. (1975). Congestion tolls for Poisson queuing processes. Econometrica: Journal of the Econometric Society, 43(1):81–92.
- Feller, [1971] Feller, W. (1971). An Introduction to Probability Theory and its Applications, Vol. II. John Wiley & Sons, New York, 2nd edition.
- Grassmann, [1983] Grassmann, W. (1983). The convexity of the mean queue size of the M/M/ queue with respect to the traffic intensity. Journal of Applied Probability, 20(4):916–919.
- Guljaš et al., [1998] Guljaš, B., Pearce, C. E., and Pečarić, J. (1998). Jensen’s inequality for distributions possessing higher moments, with application to sharp bounds for Laplace-Stieltjes transforms. The ANZIAM Journal, 40(1):80–85.
- Gupta and Osogami, [2011] Gupta, V. and Osogami, T. (2011). On Markov-Krein characterization of the mean waiting time in M/G/ and other queueing systems. Queueing Systems, 68:339–352.
- Hajek, [1983] Hajek, B. (1983). The proof of a folk theorem on queuing delay with applications to routing in networks. Journal of the ACM (JACM), 30(4):834–851.
- [19] Harel, A. (1990a). Convexity properties of the Erlang loss formula. Operations Research, 38(3):499–505.
- [20] Harel, A. (1990b). Convexity results for single-server queues and for multiserver queues with constant service times. Journal of Applied Probability, 27(2):465–468.
- Harel, [2010] Harel, A. (2010). Sharp and simple bounds for the Erlang delay and loss formulae. Queueing Systems, 64(2):119–143.
- Harel and Zipkin, [1987] Harel, A. and Zipkin, P. (1987). The convexity of a general performance measure for multiserver queues. Journal of Applied Probability, 24(3):725–736.
- Harrison and Zeevi, [2005] Harrison, J. M. and Zeevi, A. (2005). A method for staffing large call centers based on stochastic fluid models. Manufacturing & Service Operations Management, 7(1):20–36.
- Hartly and David, [1954] Hartly, H. and David, H. (1954). Universal bounds for mean range and extreme observations. The Annals of Mathematical Statistics, 25:85–99.
- Hassin, [2016] Hassin, R. (2016). Rational Queueing. CRC press, New York.
- Hassin et al., [2023] Hassin, R., Haviv, M., and Oz, B. (2023). Strategic behavior in queues with arrival rate uncertainty. European Journal of Operational Research, 309(1):217–224.
- Heemskerk et al., [2017] Heemskerk, M., van Leeuwaarden, J., and Mandjes, M. (2017). Scaling limits for infinite-server systems in a random environment. Stochastic Systems, 7(1):1–31.
- Jagers and Van Doorn, [1986] Jagers, A. and Van Doorn, E. A. (1986). On the continued Erlang loss function. Operations Research Letters, 5(1):43–46.
- Jongbloed and Koole, [2001] Jongbloed, G. and Koole, G. (2001). Managing uncertainty in call centres using poisson mixtures. Applied Stochastic Models in Business and Industry, 17(4):307–318.
- Karlin and Studden, [1966] Karlin, S. and Studden, W. J. (1966). Tchebycheff Systems: With Applications in Analysis and Statistics. John Wiley & Sons Interscience Publishers, New York.
- Kemperman, [1971] Kemperman, J. (1971). Moment problems with convexity conditions I. In Optimizing Methods in Statistics, pages 115–178. Academic Press, New York.
- Kingman, [1962] Kingman, J. F. C. (1962). Some inequalities for the queue GI/G/1. Biometrika, 49(3–4):315–324.
- Kreĭn and Nudel’man, [1977] Kreĭn, M. G. and Nudel’man, A. (1977). The Markov moment problem and extremal problems. American Mathematical Society, Providence, RI.
- Lee and Cohen, [1983] Lee, H. L. and Cohen, M. A. (1983). A note on the convexity of performance measures of M/M/ queueing systems. Journal of Applied Probability, 20(4):920–923.
- Liu and Hasenbein, [2019] Liu, C. and Hasenbein, J. J. (2019). Naor’s model with heterogeneous customers and arrival rate uncertainty. Operations Research Letters, 47(6):594–600.
- Naor, [1969] Naor, P. (1969). The regulation of queue size by levying tolls. Econometrica: Journal of the Econometric Society, 37(1):15–24.
- Pearce and Pečarić, [1995] Pearce, C. E. and Pečarić, J. (1995). An integral inequality for convex functions, with application to teletraffic congestion problems. Mathematics of Operations Research, 20(3):526–528.
- Postek et al., [2018] Postek, K., Ben-Tal, A., den Hertog, D., and Melenberg, B. (2018). Robust optimization with ambiguous stochastic constraints under mean and dispersion information. Operations Research, 66(3):814–833.
- Randhawa, [2016] Randhawa, R. S. (2016). Optimality gap of asymptotically derived prescriptions in queueing systems. Queueing Systems, 83(1):131–155.
- Rogosinski, [1958] Rogosinski, W. W. (1958). Moments of non-negative mass. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 245(1240):1–27.
- Rolski, [1972] Rolski, T. (1972). Some inequalities for GI/M/ queues. Applicationes Mathematicae, 1(13):42–47.
- Shaked and Shanthikumar, [1988] Shaked, M. and Shanthikumar, J. G. (1988). Stochastic convexity and its applications. Advances in Applied Probability, 20(2):427–446.
- Shaked and Shanthikumar, [2007] Shaked, M. and Shanthikumar, J. G. (2007). Stochastic Orders. Springer, New York.
- Shanthikumar and Yao, [1988] Shanthikumar, J. G. and Yao, D. D. (1988). Strong stochastic convexity and its applications in parametric optimization of queueing systems. In Proceedings of the 27th IEEE Conference on Decision and Control, pages 657–662. IEEE.
- Shapiro, [2001] Shapiro, A. (2001). On duality theory of conic linear problems. In Semi-Infinite programming, pages 135–165. Kluwer Academic Publishers, Dordrecht, The Netherlands.
- Smith, [1995] Smith, J. E. (1995). Generalized Chebychev inequalities: Theory and applications in decision analysis. Operations Research, 43(5):807–825.
- Steckley et al., [2009] Steckley, S. G., Henderson, S. G., and Mehrotra, V. (2009). Forecast errors in service systems. Probability in the Engineering and Informational Sciences, 23(2):305–332.
- van Eekelen et al., [2022] van Eekelen, W., den Hertog, D., and van Leeuwaarden, J. S. H. (2022). MAD dispersion measure makes extremal queue analysis simple. INFORMS Journal on Computing, 34(3):1681–1692.
- Vasicek, [1977] Vasicek, O. A. (1977). An inequality for the variance of waiting time under a general queuing discipline. Operations Research, 25(5):879–884.
- Wang et al., [2022] Wang, Y., Prasad, M. N., Hanasusanto, G. A., and Hasenbein, J. J. (2022). Distributionally robust observable strategic queues. arXiv preprint arXiv:2204.03221.
- Weber, [1980] Weber, R. R. (1980). Note—on the marginal benefit of adding servers to G/GI/ queues. Management Science, 26(9):946–951.
- Weber, [1983] Weber, R. R. (1983). A note on waiting times in single server queues. Operations Research, 31(5):950–951.
- Whitt, [1984] Whitt, W. (1984). On approximations for queues, I: Extremal distributions. AT&T Bell Laboratories Technical Journal, 63(1):115–138.
- Whitt, [2006] Whitt, W. (2006). Staffing a call center with uncertain arrival rate and absenteeism. Production and Operations Management, 15(1):88–102.
- Zan et al., [2018] Zan, J., Hasenbein, J. J., Morton, D. P., and Mehrotra, V. (2018). Staffing call centers under arrival-rate uncertainty with Bayesian updates. Operations Research Letters, 46(4):379–384.
Appendix A Proof of Lemma 2
Proof.
From the derivations in [15], we know that the derivative of can be written as
(14) |
Further differentiating both sides yields
(15) |
and
(16) | |||
By Little’s law, . The third derivative of the expected waiting time is then given by
(17) |
To prove the claim, we shall show that . Consecutively substituting (16), (15), (14) and finally (6) into (17), we obtain
(18) | ||||
Since is nonnegative, it suffices to show that
To simplify some of the terms of , it is convenient to work with bounds for in the remainder of the proof. We will use the simple bounds (i) and (ii) (see, e.g., [21]). It suffices to show for all that since the result is already demonstrated for the single-server queue. We proceed by showing that is nondecreasing for . Then, to complete the proof, it remains to show that for all . Now notice that
where
After some algebra, one sees that, for , , and clearly, . Since , has two positive roots. Denote the smaller root by . Now, to show , we demonstrate that , which is sufficient as is a convex quadratic function. We will instead prove this inequality for an upper bound on :
where the inequality follows from (ii) and . We next compute and , and show that is nonnegative. Demonstrating is sufficient since is a quadratic decreasing function on . Notice that
Since , it is sufficient to show by demonstrating nonnegativity of
After some tedious calculations, it follows from standard calculus that , for . So, for fixed , the auxiliary function is minimized at . Hence,
for all . Here the second-to-last inequality follows from minimizing . Therefore, or, equivalently, for . In the remainder of the proof, it thus suffices to concentrate on . Define
Observe that
From the well-known bound (i), it follows that for since . Further, it can be shown that
Since , for all . Hence, for ,
in which the final inequality follows from some straightforward calculus. This completes the proof. ∎
Appendix B Semi-variance as information set
We next consider a specific asymmetric dispersion measure that only measures dispersion above the mean. Let be the set of all distributions with mean , (upper) semivariance , where , and the support contained in the interval . Using primal-dual techniques, we establish the following result.
Theorem 4 (M/M/ queue with semivariance information).
Consider an M/M/ queue with random arrival rate that follows a distribution belonging to the ambiguity set . Suppose that and . Then, the tight upper bound for the expected waiting time corresponds with the maximum value of that results from the following two solutions:
-
(i)
the expected waiting time with the expectation taken over a three-point distribution with
where solves
-
(ii)
the expected waiting time with the expectation taken over the two-point distribution
Proof.
The primal problem for mean-(upper)semivariance information is given by
(19) | ||||||
s.t. |
and admits the following dual:
(20) | ||||||
s.t. |
The objective function is increasingly convex since it represents . Under the conditions imposed on the parameters of the ambiguity set, strong duality holds and the optimal values of the primal and dual problem coincide. In addition, since the common optimal value is finite, the dual optimal solution is attained [45, Proposition 3.4]. Moreover, as both the objective function and the dual function are continuous and the support is compact, the optimal primal solution is also attained [45, Corollary 3.1]. As a consequence, complementary slackness holds [45, Proposition 2.1]. We next use the complementary slackness property and structural properties of the functions and to determine which points should constitute the support of the extremal distribution. Since the primal problem has three constraints, we can restrict our search to a worst-case distribution with at most three support points [40]. Since the conditions in the theorem impose , we can readily exclude the case of a one-point worst-case distribution. Furthermore, we need both a point below and above the mean to satisfy the mean condition and to avoid the semivariance being equal to zero. For the support point below the mean, is the only dual-feasible option by convexity of , as is linear for all . We next consider two dual solutions that correspond with cases (i) and (ii) stated in the theorem, which assert a worst-case two- and three-point distribution, respectively.


(i) Fixing the first support point to , we next seek the other two points using another complementary slackness argument. The dual function can only be tangent to at one unique point on the interval . This follows from linearity of and convexity of , as illustrated in Figure 3, which rules out the possibility of a second tangent point. To see this, assume that there exists a second tangent point . To remain dual-feasible, the dual function has to satisfy for , and for . Therefore, has to intersect from below at , but since this already occurred at , this cannot happen a second time as otherwise has to be nonlinear or nonconvex. Hence, we arrive at a contradiction. The dual function can only coincide with one more time, at the upper bound of the support, . Therefore, if the worst case is given by a three-point distribution, then it has to admit this specific form as a result of complementary slackness. The probabilities, as functions of , follow from solving the moment constraints in (19). As is yet to be determined, the maximum value follows from solving the univariate optimization problem stated in the claim. There always exists a feasible solution to this optimization problem. Specifically, letting , we obtain the three-point distribution
of which is easily verified, using the conditions in the claim, that it is primal feasible for all possible parameter combinations.
(ii) We next consider the two-point solution. If the dual function only coincides with at and some , then the resulting two-point distribution can be directly derived from the moment conditions
which is a system of three equations with three unknowns with a unique solution, as stated in the claim. Further, from the assumptions in the claim, it follows that this distribution is always feasible in the primal. Taking the maximum expected value of the assertions (i) and (ii) completes the proof of the claim. ∎