Distributionally Robust Observable Strategic Queues

1. Introduction

Imposing tolls to regulate queueing systems was first studied by Naor (1969). He considers a single-server first come, first served (FCFS) queue with stationary Poisson arrivals at a known rate λ. Service times are independent and identically and exponentially distributed with the rate μ. Customers are assumed to be risk neutral and homogenous from an economic perspective. Each customer receives a reward of $\$R$ upon service completion and incurs a cost of $\$C$ per unit of time spent in the system (including in service). In the observable model, every arriving customer inspects the queue length and decides whether to join (reneging is not allowed) or balk (i.e., not join the queue). This strategic decision making is the key factor differentiating this model from the classic M/M/1 queueing model.

Naor (1969) derives an optimal threshold strategy n. The customer joins the queue if and only if the system length is less than n. He computes this threshold value under three different control strategies: (1) individual optimization (n_e) where the customers act in isolation, aiming to maximize their own expected net benefit rate; (2) social optimization (n_s) where the objective is to maximize the long-run rate at which customers accrue net benefit; and (3) revenue maximization (n_r) where the agency imposes a toll on the customers joining the queue with the goal of maximizing its own revenue. The most important result by Naor (1969) is the relation $n_r\le n_s\le n_e$, which implies that the customers tend to join the system at a higher rate when left to themselves than is socially optimal. This is because customers do not consider the negative externalities they impose on customers who arrive later. The result also implies that the revenue-maximizing firms allow fewer customers to join their system than the socially optimal case.

Many authors have expanded on the seminal work by Naor (1969); a detailed review of these game-theoretic models is presented in a recent book by Hassin and Haviv (2003). Some of the other recent works (Burnetas and Economou 2007, Economou and Kanta 2008, Guo and Hassin 2011) involve deriving threshold strategies in a classic Naor setting with server shutdowns. Although Economou and Kanta (2008) study the system with server breakdowns and repairs, Burnetas and Economou (2007) analyze the system where the server shuts off when idle and incurs a setup time to resume. A slight variant of this model is given by Guo and Hassin (2011), where the server resumes only when the queue length exceeds a given critical length. Also, Guo and Zipkin (2007) explore the effects of three different levels of delay information and identify the specific cases that do and do not require such information to improve the performance. Haviv and Oz (2016) review the properties of several existing regulation schemes and devise a new mechanism where customers are given priority based on the queue length. Afèche and Ata (2013) study the observable M/M/1 queue with heterogenous customers, with some patient and some impatient of given proportions.

All the aforementioned works explore the Naor (1969) model by assuming deterministic arrival or service rates. However, in many real-world scenarios, customers may behave differently during different periods. Hence, there is merit in building a model that performs well under uncertainty. A possible approach is to consider distributional uncertainty on the customers’ interarrival times (Bandi et al. 2015). Unfortunately, such granular data are usually inaccessible or simply not stored in practice. In addition, even if the granular data are collected, the problem is still theoretically challenging as the M/M/1 structure no longer holds. In this case, calculating the long-run expected social benefit or revenue rates would be difficult, and there is no such study in the context of strategic queues. To avoid these shortcomings, some recent studies propose to address uncertainty by taking the arrival or service rate as a random variable (Liu and Hasenbein 2019, Hassin et al. 2023). For example, they assume that the customer arrival rate is a random variable for each weekday, and the queue manager seeks a strategy that maximizes the long-run social benefits or revenue. This modeling assumption is an expedient approach to account for customer arrival variation while considering the issues of data accessibility and model complexity. Compared with granular data, historical arrival rates are much easier to obtain; most restaurants can collect the historical daily arrival rates by referring to the recorded sales quantity in their accounting books, but only a tiny portion of restaurants keep records of the exact arrival time of each customer. In addition, by assuming that the arrival rate is a fixed random variable within each time slot, the model preserves the M/M/1 structure, which enables the use of many elegant results from the paper of Naor (1969).

Papers with related assumptions as in our model include Debo and Veeraraghavan (2014), who consider a system where the arriving customers cannot completely observe the service rate and value. They assume that the server belongs to one of two known types and that the service rate and prior probability for each type are known. Liu and Hasenbein (2019) study a stochastic extension of the Naor (1969) model by relaxing the assumption of a certain arrival rate. They assume that the arrival rate is drawn from a probability distribution that is known to the decision maker. Chen and Hasenbein (2020) further extend the stochastic model to the unobservable setting. They show that the social optimizer induces a lower expected arrival rate than the revenue maximizer in this setting. Hassin et al. (2023) also investigate the unobservable stochastic model from the perspective of strategic customers and demonstrate that the model exhibits a rate-biased arrivals see time averages property. Despite their conceptual appeal, all these works require that the arrival or service rate distribution is known precisely to decision makers, which may not be realistic in practice. In this paper, we extend the classical Naor model for observable systems by relaxing these assumptions, where we assume the arrival rate is uncertain and governed by an unknown underlying distribution, whereas the service rate is deterministic.

To this end, we consider an alternate modeling paradigm called the distributionally robust optimization (DRO) (Scarf 1957, Žáčková 1966, Shapiro and Kleywegt 2002). Unlike the traditional stochastic optimization model, DRO acknowledges the lack of full distributional information on the random arrival rate. Instead, the decision maker is assumed to have access to partial information, such as the moments and structural properties of the arrival rate distribution, or some limited historical observations. In this setting, the objective is to derive optimal threshold strategies that maximize the worst-case expected benefit rate, where the worst case is taken over an ambiguity set of all distributions consistent with the available information about the true distribution. Such max-min problems have been studied since the seminal work by Scarf (1957), but they have only received more attention with the advent of modern robust optimization techniques (Bertsimas and Sim 2004, Ben-Tal et al. 2009). Since then, a substantial body of literature has been devoted to studying well-known optimization problems under uncertainty in a distributionally robust setting; see Delage and Ye (2010), Li et al. (2014), Wiesemann et al. (2014), Hanasusanto et al. (2015), Shafieezadeh-Abadeh et al. (2015), and Ardestani-Jaafari and Delage (2021). Nevertheless, the distributionally robust framework has not been considered in the context of the classical Naor observable strategic queue model. The paper fills this gap in the literature.

We first study the distributionally robust queue model with a mean-absolute deviation (MAD) ambiguity set (Postek et al. 2018), where partial information about the distribution mean and MAD are known. Next, we extend our model to the data-driven setting, where queue system managers only have access to a finite number of independent and identically distributed training samples collected from historical observations. We construct a data-driven mean-absolute deviation (DD-MAD) ambiguity set that mitigates estimation errors from the empirical moment estimators. The resulting distributionally robust model with a data-driven ambiguity set admits a semidefinite programming (SDP) reformulation for the social optimization problem and a linear programming reformulation for the revenue maximization problem. To properly determine the robustness parameters, we establish a new distribution-free confidence interval for the empirical MAD. Although such confidence intervals exist for the empirical mean and variance (Delage and Ye 2010), to the best of our knowledge, none are available for the empirical MAD. Herrey (1965) derives the confidence interval for the empirical MAD under normal distribution data, whereas other works mostly focus on median-absolute deviation; see Bonett and Seier (2003), Abu-Shawiesh et al. (2018), and Arachchige and Prendergast (2019). Using this result, we further derive finite-sample guarantees for the data-driven MAD model, in which optimal values provide high-confidence lower bounds on the expected social benefit or revenue rate. We also benchmark our data-driven MAD ambiguity set with the popular Wasserstein ambiguity set (Pflug and Wozabal 2007, Esfahani and Kuhn 2018, Esfahani et al. 2018, Gao and Kleywegt 2023), which is widely used in the data-driven setting as it can offer attractive finite-sample guarantees. The results demonstrate that our proposed data-driven MAD model shares a similar guarantee as the Wasserstein model while generating a significantly more tractable reformulation.

The main contributions of this paper can be summarized as follows.

1. We propose a new model to tackle the uncertain arrival rate in the Naor strategic queue problem using the emerging DRO framework. The model does not impose any specific distributional assumption; instead, it optimizes in view of the worst-case distribution within a prescribed ambiguity set. Benefitting from this robustification framework, the model alleviates the overfitting issue and yields attractive out-of-sample performance.

2. We prove that the revenue rate function is concave, whereas the social benefit rate function is either concave or unimodal under some mild prerequisites. We then show that these properties enable a closed-form solution for the worst-case expectation problem with an MAD ambiguity set. For the general cases, we derive an SDP reformulation for the social optimization problem and a linear programming reformulation for the revenue optimization problem.

3. We extend the distributionally robust model to the data-driven setting, where queue system managers only have access to a finite set of historical observations. To mitigate the adverse effect of the estimation errors from the empirical MAD, we robustify the ambiguity set by adding an extra layer of robustness to the empirical mean and MAD estimators. The data-driven MAD model admits an SDP reformulation for the social optimization problem and a linear programming reformulation for the revenue maximization problem. We then establish a distribution-free confidence interval for the empirical MAD and derive finite-sample guarantees for the distributionally robust model with a data-driven MAD ambiguity set. Compared with the Wasserstein ambiguity set, the data-driven MAD ambiguity set admits a more efficient reformulation of fixed complexity, where the number of constraints does not scale with the sample size.

The remainder of the paper is structured as follows. In Section 2, we propose the distributionally robust queue model and analyze the relationship between different thresholds under the distributionally robust setting. Section 3 presents tractable reformulations for the worst-case expectation problem with a classical MAD ambiguity set. Section 4 explores the distributionally robust model with a data-driven MAD ambiguity set and derives theoretical finite-sample guarantees. Finally, the out-of-sample performances of our distributionally robust models are assessed empirically in Section 5.

1.1. Notations

The set of all probability measures supported on Ξ is written as ${\mathcal{P}}_0(\mathrm{\Xi })≔\{\mu \in {\mathcal{M}}_{+}:{\int }_{\mathrm{\Xi }}\mu (d\xi )=1\}$, where ${\mathcal{M}}_{+}$ denotes the set of nonnegative Borel measures. All random variables are designated by tilde signs (e.g., $\tilde{\rho }$), whereas their realizations are denoted without tildes (e.g., ρ). We denote by ${\mathbb{E}}_{\mathbb{P}}[c(\tilde{\rho })]$ the expectation of a cost function with respect to the random variable $\tilde{\rho }$ under distribution $\mathbb{P}$. We define $\lfloor n\rfloor$ to be the largest integer less than or equal to n and $\Vert \mathit{x}{\Vert }_p$ to be the p-norm of a vector $\mathit{x}$. For any set Ξ, we let $\text{int}(\mathrm{\Xi })$ denote its interior. The cone of k × k positive semidefinite matrices is denoted by ${\mathbb{S}}_{+}^k$.

2. Distributionally Robust Strategic Queues Model

The extension of the Naor (1969) seminal queue model to the stochastic optimization setting with an uncertain arrival rate was first proposed by Liu and Hasenbein (2019), who consider an M/M/1 queue system with a random arrival rate $\tilde{\lambda }\sim {\mathbb{P}}^{\star }$ and a deterministic service rate μ. The queue system operates under a first come, first served discipline, and the true distribution of the uncertain arrival rate $\tilde{\lambda }$ is known by the system manager. Because the service rate μ is deterministic, without loss of generality, we consider the traffic intensity $\tilde{\rho }≔\frac{\tilde{\lambda }}{\mu }$ as the uncertain parameter throughout the remainder of the paper. The stochastic model aims to find an optimal threshold that maximizes the expected benefit rate: that is,

Here, $c_n(\tilde{\rho })$ is a general return function, which can be replaced with the social benefit rate function or the revenue rate function depending on the system manager’s objective.

In practice, the true distribution ${\mathbb{P}}^{\star }$ is never available to the system manager and typically has to be estimated using the empirical distribution generated from the historical observations. Although the empirical-based methods may work well on the observed data set, they often fail to achieve an acceptable out-of-sample performance because they do not consider any possible disturbances from the limited historical observations.

In this paper, we endeavor to address this fundamental shortcoming using ideas of DRO. The DRO approach does not impose any single distribution on the uncertain arrival rate. Instead, it constructs an ambiguity set $\mathcal{P}$ containing all plausible probability distributions that are consistent with the partial information as well as historical observations. In this setting, the objective is to derive an optimal threshold strategy $\widehat{n}$ that maximizes the worst-case expected benefit rate, where the worst case is taken over all distributions from within this ambiguity set: that is,

Because the model optimizes the expected benefit rate in view of the worst-case distribution, it mitigates overfitting to the observed samples and helps improve the performance in out-of-sample circumstances.

In this paper, we study the distributionally robust model from the perspective of an individual customer, a social optimizer, and a revenue maximizer. We first derive the results that hold for any generic ambiguity set $\mathcal{P}$.

2.1. Individual Optimization

We determine a pure threshold strategy in which each arriving customer decides to join or not join the queue based on the observed queue length, independent of the strategy adopted by other customers. A newly arrived customer makes a decision (to join or not join) based on the net gain $R-(i+1)C/\mu$, where i is the number of people currently in the queue, and will join the queue if it is nonnegative. Note that net gain is deterministic because it is independent of the random arrival rate. Thus, the optimal joining threshold for any arriving customer is given by

This result coincides with the original result of Naor (1969) (i.e., ${\widehat{n}}_e=n_e$) because the net gain of a newly arrived customer only depends on the current queue length and the service rate, which are all deterministic. On the other hand, as an individual optimizer, the customer can ignore the rates of future arrivals because they will not affect the time to service.

2.2. Social Optimization

We next analyze the distributionally robust threshold for a social optimizer. The social benefit rate for a realization of the traffic intensity ρ and a fixed threshold n is given by

One can verify that ${\lim }_{\rho \to 1}R\mu \frac{\rho (1-{\rho }^n)}{1-{\rho }^{n+1}}-C\left(\frac{\rho }{1-\rho }-\frac{(n+1){\rho }^{n+1}}{1-{\rho }^{n+1}}\right)=R\mu \frac{n}{n+1}-C\frac{n}{2}$, which indicates that the function $f_n(\rho )$ is continuous in ρ. The distributionally robust model determines an optimal threshold ${\widehat{n}}_s$ that maximizes the worst-case expected social benefit rate $Z_s(n)$: that is, ${\widehat{n}}_s\in {\text{arg}\hspace{0.17em}\text{max}}_{n\in {\mathbb{Z}}_{+}}{Z}_s(n)$, where

We first investigate the relationship between the optimal thresholds ${\widehat{n}}_e$ and ${\widehat{n}}_s$.

Proposition 1.

There exists an optimal threshold of the social optimizer less than or equal to the optimal threshold of an individual customer: that is, $\exists {\widehat{n}}_s\hspace{0.17em}s.t.\hspace{0.17em}{\widehat{n}}_s\le {\widehat{n}}_e.$

Proposition 1 enables decision makers to search for the best threshold from $\{1,\dots ,{\widehat{n}}_e\}$. We remark that the participation of customers is not affected by the distributionally robust setting because the queue adopts the FCFS discipline, so subsequent arrivals are immaterial after the customer has joined the queue. Thus, the socially optimal threshold is always achievable by reducing the individual threshold from n_e to ${\widehat{n}}_s$.

2.3. Revenue Optimization

We now consider a profit-maximizing firm that aims to maximize its expected revenue rate by imposing a toll t on every joining customer. In this setting, customers base their joining decision on this imposed toll t and evaluate the service completion only by R − t. Recall that customers join the queue if and only if the expected net gain is nonnegative. Therefore, determining an optimal toll t is equivalent to choosing a queue-length threshold n that maximizes the expected revenue rate, where $n=\lfloor \frac{(R-t)\mu }{C}\rfloor .$

The revenue rate for a realization of the traffic intensity and a fixed threshold n is given by

One can show that ${\lim }_{\rho \to 1}\frac{\rho (1-{\rho }^n)}{1-{\rho }^{n+1}}=\frac{n}{n+1}$, which indicates that $f_n(\rho )$ is continuous. The distributionally robust model determines an optimal threshold ${\widehat{n}}_r$ that maximizes the worst-case expected revenue rate $Z_r(n)$: that is, ${\widehat{n}}_r\in {\text{arg}\hspace{0.17em}\text{max}}_{n\in {\mathbb{Z}}_{+}}{Z}_r(n)$, where

Similarly, we first investigate the relationship between the optimal thresholds ${\widehat{n}}_e$ and ${\widehat{n}}_r$.

Proposition 2.

There exists an optimal threshold of the revenue maximizer less than or equal to the optimal threshold of an individual customer: that is, $\exists {\widehat{n}}_r\hspace{0.17em}s.t.\hspace{0.17em}{\widehat{n}}_r\le {\widehat{n}}_e.$

So far, we have presented the generic distributionally robust observable queue models for an individual customer, a social optimizer, and a revenue maximizer. However, we have not specified the ambiguity set for the social and revenue optimization problems. In the following sections, we will investigate different types of ambiguity sets and derive their tractable reformulations.

3. Distributionally Robust Strategic Queues with an MAD Ambiguity Set

We study the DRO model with an MAD ambiguity set. Suppose the support $[a,b]$, mean m, and MAD d of the random parameter $\tilde{\rho }$ are known to the decision makers. Then, we can construct an ambiguity set containing all possible distributions that are consistent with the partial information, defined as

We develop efficient solution schemes to find the optimal threshold strategies for a social optimizer and a revenue maximizer, given by ${\widehat{n}}_s$ and ${\widehat{n}}_r$, respectively, such that the worst-case expected benefit rates are maximized. In order to derive tractable reformulations for the distributionally robust models, we assume $m\in (a,b)$ and $d\in (0,\overline{d})$, where $\overline{d}≔\frac{2(m-a)(b-m)}{b-a}$ is the largest possible mean-absolute deviation attained by any distribution with the given support and mean.

3.1. Social Optimization

To determine an optimal joining threshold for a social optimizer, we compute the worst-case expected social benefit rate $Z_s(n)$ for every $n\in {\mathbb{Z}}_{+}$ satisfying $1\le n\le n_e$ and choose an ${\widehat{n}}_s$ such that ${\widehat{n}}_s\in {\text{arg}\hspace{0.17em}\text{max}}_{n\in {\mathbb{Z}}_{+}}{Z}_s(n)$. To this end, we show how to compute the worst-case expected social benefit rate for a fixed n. Suppose the distribution mean and MAD of $\tilde{\rho }$ are precisely known; then, the worst-case expected social benefit rate is given by the optimal value of the moment problem

where $\mathrm{\Xi }≔[a,b]$ is the support of $\tilde{\rho }$. The third constraint of (8) restricts the nonnegative measure ν to be a probability distribution, whereas the first and second constraints impose that the distribution’s MAD and mean are equal to d and m, respectively. The objective of the problem is to find a feasible distribution that minimizes the expected social benefit rate.

The semi-infinite linear optimization Problem (8) is hard to solve because it searches for the best decision from an infinite-dimensional space of probability measures. To derive a tractable reformulation, we focus on the dual problem. We first define $F(\rho )≔\alpha |\rho -m|+\beta \rho +\gamma$ and derive the dual problem as

Notice that $F(\rho )$ is a two-piece piecewise affine function majorized by $f_n(\rho )$. We know that if $f_n(\rho )$ is a piecewise affine function or a concave function, the semi-infinite constraint will reduce to a linear constraint because we only need to check the satisfaction of the constraint at points $\rho =a,m$, and b. However, the social benefit rate function is neither concave nor piecewise affine, making the problem difficult. To solve this optimization problem, we first investigate the properties of the social benefit rate function $f_n(\rho )$. For clarity of exposition, we relegate some of the proofs to Appendix B.

Lemma 1.

The social benefit rate function $f_n(\rho )$ has the following properties if $\frac{R\mu }{C}\ge n+1$.

1. $f_n(\rho )$ is strictly concave for $\rho \in [0,1]$.

2. $f_n(\rho )$ is either concave increasing or unimodal for $\rho \in [0,\infty )$.

3. The sign of the second derivative $f_n^{″}(\rho )$ changes at most once over $[0,\infty )$.

From Lemma 1, we know that the social benefit rate function has some appealing properties. Specifically, the function is either concave increasing or unimodal on the nonnegative axis, and when it is unimodal, the function changes from a concave function to a convex function at some point. The next lemma further asserts that the complementary slackness property holds for the primal and dual problems, which will later help us determine the worst-case distribution.

Lemma 2.

The optimal values of the primal-dual pair (8) and (9) coincide, and their optimal solutions ${\nu }^{\star }$ and $({\alpha }^{\star },{\beta }^{\star },{\gamma }^{\star })$, respectively, satisfy the complementary slackness condition

$$(f_n(\rho )-{\alpha }^{\star }|\rho -m|-{\beta }^{\star }\rho -{\gamma }^{\star }){\nu }^{\star }(\text{d}\rho )=0\hspace{1em}\hspace{1em}\forall \rho \in [a,b].$$

The proofs of the lemmas are relegated to Appendix B. Combining Lemmas 1 and 2, we are ready to show that Problem (8) can be solved analytically under certain conditions. Specifically, we divide this problem into three cases and derive an explicit expression of the worst-case distribution for each case.

Proposition 3.

Assume $m\in [0,1]$ and $\frac{R\mu }{C}\ge n+1$. Let $({\rho }_t,f_n({\rho }_t))$ be the tangent point on f_n for the line that passes through $(m,f_n(m))$. For any $n\ge 1$, we have one of the following three cases.

1. If $f_n(b)+f_n^{\prime }(b)(m-b)\ge f_n(m)$, then the extremal distribution that solves (4) is a three-point distribution supported on ${\rho }_1=a,\hspace{0.17em}{\rho }_2=m,\hspace{0.17em}{\rho }_3=b$, with corresponding probabilities

   $$p_1=\frac{d}{2(m-a)},\hspace{0.17em}{p}_2=1-\frac{d}{2(m-a)}-\frac{d}{2(b-m)},\hspace{0.17em}{p}_3=\frac{d}{2(b-m)}.$$

2. If $f_n(b)+f_n^{\prime }(b)(m-b)<f_n(m)$ and $d<d_0≔\frac{2(m-a)({\rho }_t-m)}{{\rho }_t-a}$, then the extremal distribution is a three-point distribution supported on ${\rho }_1=a,\hspace{0.17em}{\rho }_2=m,\hspace{0.17em}{\rho }_3={\rho }_t$, with probabilities

   $$p_1=\frac{d}{2(m-a)},\hspace{0.17em}{p}_2=1-\frac{d}{2(m-a)}-\frac{d}{2({\rho }_t-m)},\hspace{0.17em}{p}_3=\frac{d}{2({\rho }_t-m)}.$$

3. If $f_n(b)+f_n^{\prime }(b)(m-b)<f_n(m)$ and $d\ge d_0≔\frac{2(m-a)({\rho }_t-m)}{{\rho }_t-a}$, then the extremal distribution is a two-point distribution supported on ${\rho }_1=a,\hspace{0.17em}{\rho }_2=\frac{ad+2m(a-m)}{d+2(a-m)}$, with probabilities

   $$p_1=\frac{d}{2(m-a)},\hspace{0.17em}{p}_2=1-\frac{d}{2(m-a)}.$$

Figure 1 depicts the optimal two-piece piecewise affine function described in Proposition 3. We remark that the tangent point $({\rho }_t,f_n({\rho }_t))$ in Figure 1(b) can be determined efficiently by the bisection method. Specifically, we set $[l,u]=[m,b]$ as the initial search interval for the algorithm. In each iteration, we compute the derivative at the midpoint $\rho =\frac{u+l}{2}$, and we check whether it is the tangent point by calculating the difference between $f_n(m)$ and $f_n^{\prime }\left(\frac{u+l}{2}\right)\left(m-\frac{u+l}{2}\right)+f_n\left(\frac{u+l}{2}\right)$. If the difference is small enough, we terminate the algorithm; otherwise, we set $u=\frac{u+l}{2}$ if the difference is positive or set $l=\frac{u+l}{2}$ if the difference is negative, and then, we go back to the first step with the updated interval $[l,u]$.

We remark that the use of the MAD ambiguity set and its geometric interpretation is motivated by a recent work by van Eekelen et al. (2022), who analyze the worst-case performance of the GI/G/1 queue (Bhat 2008) under mean-dispersion constraints. The authors demonstrate that measuring the dispersion by MAD, instead of variance, significantly simplifies the analysis and enables a closed-form solution for the extremal distribution whenever the loss function is convex. Unfortunately, our problem is different as the social benefit rate function is neither convex nor concave. Nevertheless, by establishing some useful properties of the social benefit rate function and exploiting its geometric interpretation in the dual Problem (9), we are able to explicitly express the extremal distribution when $m\le 1$ and $\frac{R\mu }{C}\ge n+1$. Using this result, we can compute the worst-case expected social benefit rate $Z_s(n)$ efficiently.

Theorem 1.

Assume $m\in [0,1]$ and $\frac{R\mu }{C}\ge n+1$. Let $({\rho }_t,f_n({\rho }_t))$ be the tangent point on $f_n(\rho )$ for the line that passes through $(m,f_n(m))$. For any $n\ge 1$, we have the following three cases.

1. If $f_n(b)+f_n^{\prime }(b)(m-b)\ge f_n(m)$, then

   $$Z_s(n)=\frac{d}{2(m-a)}{f}_n(a)+\left(1-\frac{d}{2(m-a)}-\frac{d}{2(b-m)}\right)f_n(m)+\frac{d}{2(b-m)}{f}_n(b).$$

2. If $f_n(b)+f_n^{\prime }(b)(m-b)<f_n(m)$ and $d<d_0≔\frac{2(m-a)({\rho }_t-m)}{{\rho }_t-a}$, then

   $$Z_s(n)=\frac{d}{2(m-a)}{f}_n(a)+\left(1-\frac{d}{2(m-a)}-\frac{d}{2({\rho }_t-m)}\right)f_n(m)+\frac{d}{2({\rho }_t-m)}{f}_n({\rho }_t).$$

3. If $f_n(b)+f_n^{\prime }(b)(m-b)<f_n(m)$ and $d\ge d_0≔\frac{2(m-a)({\rho }_t-m)}{{\rho }_t-a}$, then

   $$Z_s(n)=\frac{d}{2(m-a)}{f}_n(a)+\left(1-\frac{d}{2(m-a)}\right)f_n\left(\frac{ad+2m(a-m)}{d+2(a-m)}\right).$$

Theorem 1 enables us to solve the worst-case expectation problem analytically under certain conditions. However, for the more general case, we are unable to solve it in a closed form. In the following theorem, we show that the worst-case expectation problem admits a semidefinite programming reformulation that can be solved in polynomial time using standard off-the-shelf solvers, such as SDPT3 (Toh et al. 1999) and MOSEK (ApS 2022).

Theorem 2.

For any $n\ge 1$, the worst-case expected social benefit rate $Z_s(n)$ coincides with the optimal value of the following semidefinite program:

$$\begin{array}{llll}\hfill & \mathit{\text{sup}}\hspace{0.17em}\hspace{0.17em}\hfill & \alpha d+\beta m+\gamma \hfill & \hfill \\ \hfill & \mathit{s}.\mathit{t}.\hspace{0.17em}\hspace{0.17em}\hfill & \alpha ,\beta ,\gamma \in \mathbb{R},\mathit{y},\mathit{z}\in {\mathbb{R}}^{n+3},\mathit{X},\mathit{X}\prime \in {\mathbb{S}}_{+}^{n+3}\hfill & \hfill \\ \hfill & \hfill & y_1=R\mu -C-y_0+y_{n+3},\hspace{0.17em}{y}_2=-R\mu -y_{n+3},\hfill & \hfill \\ \hfill & \hfill & y_3,\dots ,y_n=0,\hspace{0.17em}{y}_{n+1}=-R\mu +C(n+1)-y_0,\hfill & \hfill \\ \hfill & \hfill & y_{n+2}=R\mu -Cn+y_0-y_{n+3},\hspace{0.17em}\hfill & \hfill \\ \hfill & \hfill & y_0=\alpha m+\gamma ,y_{n+3}=-\alpha +\beta \hspace{0.17em}\hfill & \hfill \\ \hfill & \hfill & {\displaystyle \sum _{i+j=2l-1}{x}_{ij}}=0\hfill & \forall l\in [n+4]\hfill \\ \hfill & \hfill & {\displaystyle \sum _{i+j=2l}{x}_{ij}}={\displaystyle \sum _{q=0}^l{\displaystyle \sum _{r=q}^{n+3+q-l}{y}_r}}\left(\begin{array}{c}r\\ q\end{array}\right)\left(\begin{array}{c}n+3-r\\ l-q\end{array}\right)a^{r-q}{m}^q\hfill & \forall l\in [n+4]\cup \left\{0\right\}\hfill \\ \hfill & \hfill & z_1=R\mu -C-z_0+z_{n+3},\hspace{0.17em}{z}_2=-R\mu -z_{n+3},\hfill & \hfill \\ \hfill & \hfill & z_3,\dots ,z_n=0,\hspace{0.17em}{z}_{n+1}=-R\mu +C(n+1)-z_0,\hfill & \hfill \\ \hfill & \hfill & z_{n+2}=R\mu -Cn+z_0-z_{n+3}\hspace{0.17em}\hfill & \hfill \\ \hfill & \hfill & z_0=-\alpha m+\gamma ,z_{n+3}=\alpha +\beta \hspace{0.17em}\hfill & \hfill \\ \hfill & \hfill & {\displaystyle \sum _{i+j=2l-1}{x}_{ij}^{\prime }}=0\hfill & \forall l\in [n+4]\hfill \\ \hfill & \hfill & {\displaystyle \sum _{i+j=2l}{x}_{ij}^{\prime }}={\displaystyle \sum _{q=0}^l{\displaystyle \sum _{r=q}^{n+3+q-l}{y}_r^{\prime }}}\left(\begin{array}{c}r\\ q\end{array}\right)\left(\begin{array}{c}n+3-r\\ l-q\end{array}\right)m^{r-q}{b}^q\hfill & \forall l\in [n+4]\cup \left\{0\right\}.\hfill \end{array}$$(10)

The proof of this theorem relies on the following lemma, which expresses a univariate polynomial inequality in terms of semidefinite constraints.

Lemma 3

(Bertsimas and Popescu 2005, Proposition 3.1(f)). The polynomial $g(\rho )={\sum }_{r=0}^k{y}_r{\rho }^r$ satisfies $g(\rho )\ge 0$ for all $\rho \in [a,b]$ if and only if there exists a positive semidefinite matrix $X={[x_{ij}]}_{i,j=0,\dots ,k}\in {\mathbb{S}}_{+}^{k+1}$ such that

$$\begin{array}{ll}0=\hfill & {\displaystyle \sum _{i,j:i+j=2l-1}{x}_{ij}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}\hspace{0.5em}\forall l=1,\dots ,k\hfill \\ {\displaystyle \sum _{q=0}^l{\displaystyle \sum _{r=q}^{k+q-l}{y}_r}}\left(\begin{array}{c}r\\ q\end{array}\right)\left(\begin{array}{c}k-r\\ l-q\end{array}\right)a^{r-q}{b}^q=\hfill & {\displaystyle \sum _{i,j:i+j=2l}{x}_{ij}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall l=0,\dots ,k.\hfill \end{array}$$

Proof of Theorem 2.

Recall that the dual of ${\inf }_{\mathbb{P}\in \mathcal{P}}{\mathbb{E}}_{\mathbb{P}}[f_n(\tilde{\rho })]$ for $\tilde{\rho }$ supported on the interval $[a,b]$ is given by (see Problem (9))

$$\begin{array}{lll}\hfill & \underset{\alpha ,\beta ,\gamma \in \mathbb{R}}{\sup }\hfill & \alpha d+\beta m+\gamma \hfill \\ \hfill & \text{s}.\text{t}.\hfill & \alpha |\rho -m|+\beta \rho +\gamma \le f_n(\rho )\hspace{1em}\hspace{1em}\forall \rho \in [a,b].\hfill \end{array}$$

We can deal with the semi-infinite constraint separately for the cases $\rho \le m$ and $\rho \ge m$:

$$\begin{array}{ll}\hfill & \underset{\alpha ,\beta ,\gamma \in \mathbb{R}}{\sup }\hspace{0.17em}\hspace{0.5em}\alpha d+\beta m+\gamma \hfill \\ \hfill & \text{s}.\text{t}.\hspace{1em}\alpha (m-\rho )+\beta \rho +\gamma \le f_n(\rho )\hspace{1em}\hspace{1em}\forall \rho \in [a,m]\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}\alpha (\rho -m)+\beta \rho +\gamma \le f_n(\rho )\hspace{1em}\hspace{1em}\forall \rho \in [m,b].\hfill \end{array}$$

Substituting the definition of $f_n(\rho )$ in (3) and applying algebraic reductions yield the following polynomial inequalities:

$$\begin{array}{ll}\hfill & -(\alpha m+\gamma ){\rho }^0+(R\mu -C-\beta +\gamma +\alpha m+\alpha )\rho +(-R\mu -\alpha +\beta ){\rho }^2+(-R\mu +Cn+C+\alpha m+\gamma ){\rho }^{n+1}\hfill \\ \hfill & \hspace{1em}\hspace{1em}+(R\mu -Cn-\alpha m-\alpha +\beta -\gamma ){\rho }^{n+2}+(\alpha -\beta ){\rho }^{n+3}\ge 0\hspace{1em}\hspace{1em}\hspace{0.17em}\forall \rho \in [a,m],\hspace{1em}\text{and}\hfill \\ \hfill & (\alpha m-\gamma ){\rho }^0+(R\mu -C-\alpha m-\alpha -\beta +\gamma )\rho +(-R\mu +\alpha +\beta ){\rho }^2+(-R\mu +Cn+C-\alpha m+\gamma ){\rho }^{n+1}\hfill \\ \hfill & \hspace{1em}\hspace{1em}+(R\mu -Cn+\alpha m+\alpha +\beta -\gamma ){\rho }^{n+2}-(\alpha +\beta ){\rho }^{n+3}\ge 0\hspace{1em}\hspace{1em}\forall \rho \in [m,b].\hfill \end{array}$$(11)

The inequalities are of the form $g_1(\rho )={\sum }_{r=0}^{n+3}{y}_r{\rho }^r\ge 0$ for $\rho \in [a,m]$ and $g_2(\rho )={\sum }_{r=0}^{n+3}{z}_r{\rho }^r\ge 0$ for $\rho \in [m,b]$, where $y=(y_1,\dots ,y_{n+3})$ and $z=(z_1,\dots ,z_{n+3})$ represent the coefficients of the respective polynomial inequalities. We now invoke the result of Lemma 3 with $k=n+3$ to express the inequalities in (11) as semidefinite constraints. The resulting semidefinite problem is equivalent to the original problem, which completes the proof. □

Remark 1.

In this subsection, we present two results. Theorem 1 provides a closed-form solution under certain prerequisites, whereas Theorem 2 derives an SDP reformulation for the general cases. It is worth noting that Theorem 1 requires the parameters to satisfy $n\le \frac{R\mu }{C}-1$. By Proposition 1, there exists an optimal threshold ${\widehat{n}}_s$ less than or equal to ${\widehat{n}}_e$ (i.e., $\exists {\widehat{n}}_s\le {\widehat{n}}_e=\lfloor \frac{R\mu }{C}\rfloor$). Thus, for a strategic queue with mean arrival rate $m\le 1$, Theorem 1 can be applied to compute the worst-case expected social benefit rate for the first $\lfloor \frac{R\mu }{C}\rfloor -1$ cases. This greatly speeds up the time to solve (4) because we only need to solve an SDP once for the remaining case $n=\lfloor \frac{R\mu }{C}\rfloor$. On the other hand, for a strategic queue with mean arrival rate m > 1, we cannot invoke Theorem 1 anymore, and we need to solve an SDP for each n satisfying $1\le n\le {\widehat{n}}_e,\hspace{0.17em}n\in {\mathbb{Z}}_{+}$.

3.2. Revenue Optimization

To determine an optimal joining threshold for a revenue maximizer, we compute the worst-case expected revenue rate $Z_r(n)$ for every $n\in {\mathbb{Z}}_{+},\hspace{0.17em}1\le n\le n_e$, and we choose an ${\widehat{n}}_r$ such that ${\widehat{n}}_r\in {\text{arg}\hspace{0.17em}\text{max}}_{n\in {\mathbb{Z}}_{+}}\{Z_r(n)\}$. To this end, we show how to compute the worst-case expected revenue for each n. Suppose the mean and MAD of the random parameter $\tilde{\rho }$ are known; then, the worst-case expected revenue rate is given by the following optimization problem:

To derive a tractable reformulation, we first investigate the property of the revenue rate function $r_n(\rho )$.

Lemma 4.

The revenue rate function $r_n(\rho )$ is concave for $\rho \in {\mathbb{R}}_{+}$.

Equipped with Lemma 4, we now show that the worst-case expectation Problem (12) admits a closed-form solution.

Theorem 3.

For any $n\ge 1$, the worst-case expected revenue rate can be derived as

$$Z_r(n)=\frac{d}{2(m-a)}{r}_n(a)+\left(1-\frac{d}{2(m-a)}-\frac{d}{2(b-m)}\right)r_n(m)+\frac{d}{2(b-m)}{r}_n(b).$$

To prove this theorem, we invoke a classical result that characterizes the worst-case distribution from the MAD ambiguity set for a concave loss function.

Lemma 5

(Ben-Tal and Hochman 1972, Theorem 3). Suppose $f(\rho )$ is a concave function and the ambiguity set is defined as $\mathcal{P}=\{\mathbb{P}\in {\mathcal{P}}_0([a,b]):\hspace{0.17em}{\mathbb{E}}_{\mathbb{P}}[\tilde{\rho }]=m,\hspace{0.17em}{\mathbb{E}}_{\mathbb{P}}[|\tilde{\rho }-m|]=d\}$. The extremal distribution that solves ${\inf }_{\mathbb{P}\in \mathcal{P}}{\mathbb{E}}_{\mathbb{P}}[f(\tilde{\rho })]$ is a three-point distribution supported on ${\rho }_1=a,\hspace{0.17em}{\rho }_2=m,\hspace{0.17em}{\rho }_3=b$ with probabilities

$$p_1=\frac{d}{2(m-a)},\hspace{0.17em}{p}_2=1-\frac{d}{2(m-a)}-\frac{d}{2(b-m)},\hspace{0.17em}{p}_3=\frac{d}{2(b-m)}.$$(13)

Proof of Theorem 3.

From Lemma 4, the revenue rate function $r_n(\rho )$ is concave. Therefore, applying Lemma 5 yields the result. □

Theorem 3 provides practical managerial insight for the decision maker. Observe that the extremal distribution is usually a discrete distribution supported on the mean and the lower and upper bounds of the support. Hence, instead of optimizing over the empirical distribution, the DRO scheme simplifies the problem into three cases: when the traffic intensity is extremely small ($\rho =a$), extremely large ($\rho =b$), or as expected ($\rho =m$). The weight for each scenario is determined by the MAD, which reflects the variation level of samples. This result aligns well with human intuition. To design a robust queue-regulating strategy, the decision maker may intuitively think about “how the queue behaves when the traffic intensity is extremely large, small, or as usual” and “what is the probability of these scenarios happening.” The closed-form solution provides an answer to these questions. For example, the quantities $p_1,p_2,p_3$ answer the question “What are the probabilities of these scenarios happening?,” whereas the quantities $r_n(a),r_n(b),r_n(m)$ answer the question “How does the queue behave when the traffic intensity is extremely large, small, or as expected?”

4. Extension to Data-Driven Problems

In this section, we apply the MAD ambiguity set to data-driven optimization problems. As we observed in the previous section, distributionally robust models with a moment ambiguity set necessitate decision makers to have access to exact values of the mean m or MAD d of the true unknown distribution, which may not be realistic in practice. A common approach is to construct such moment ambiguity sets by plugging in the point estimators generated from the historical samples. However, it is rarely the case that one can be entirely confident in these empirical estimators. For example, when the sample size is small, these empirical estimators might be far away from the true values; furthermore, some estimators, such as the empirical MAD, are even biased. In order to mitigate the adverse effects of the estimation errors, we design a data-driven MAD ambiguity set that contains the true underlying distribution with high confidence.

Unlike the setting in the previous section, here we assume that the queue system manager only has access to N independent and identically distributed samples of the traffic intensity given by ${\{{\widehat{\rho }}_i\}}_{i\in [N]}$, where ${\widehat{\rho }}_i={\widehat{\lambda }}_i/\mu$. In addition, we assume that decision makers have some prior knowledge or an educated estimate of the distribution support. Suppose the true mean and MAD of the underlying distribution are unknown and with high probabilities, belong to two confidence intervals $\mathcal{T}=[m_l,m_u]$ and $\mathcal{D}=[d_l,d_u]$ constructed using the samples. Then, the proposed data-driven distributionally robust model is formulated as

where the modified data-driven ambiguity set is defined as

One can verify that the results of Propositions 1 and 2 still hold, and we can obtain the optimal value of (14) by solving $\overline{Z}(n)$ for each $n\in {\mathbb{Z}}_{+}$ satisfying $1\le n\le n_e$ and select the one with the largest objective value.

We now derive the reformulations for the worst-case expected social benefit and revenue rates. To this end, we define the worst-case expected social benefit rate with the data-driven MAD ambiguity set by

and we define the worst-case expected revenue rate with the data-driven MAD ambiguity set by

The next theorem presents the reformulation of the worst-case expected social benefit rate. Some of the proofs of the results in this section are relegated to Appendix C.

Theorem 4.

For any $n\ge 1$, the worst-case expected social benefit rate ${\overline{Z}}_s(n)$ coincides with the optimal value of the following semidefinite program:

$$\begin{array}{llll}\hfill & \mathit{\text{sup}}\hspace{0.17em}\hspace{0.17em}\hfill & \gamma +{\theta }_1{d}_l-{\theta }_2{d}_u+{\theta }_3{m}_l-{\theta }_4{m}_u\hfill & \hfill \\ \hfill & \mathit{s}\mathit{.}\mathit{t}\mathit{.}\hspace{0.17em}\hspace{0.17em}\hfill & \gamma \in \mathbb{R},{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4\in {\mathbb{R}}_{+},\mathit{y},\mathit{z}\in {\mathbb{R}}^{n+3},\mathit{X},\mathit{X}\prime \in {\mathbb{S}}_{+}^{n+3}\hfill & \hfill \\ \hfill & \hfill & y_1=R\mu -C-y_0+y_{n+3},\hspace{0.17em}{y}_2=-R\mu -y_{n+3},\hfill & \hfill \\ \hfill & \hfill & y_3,\dots ,y_n=0,\hspace{0.17em}{y}_{n+1}=-R\mu +C(n+1)-y_0,\hfill & \hfill \\ \hfill & \hfill & y_{n+2}=R\mu -Cn+y_0-y_{n+3},\hspace{0.17em}\hfill & \hfill \\ \hfill & \hfill & y_0=({\theta }_1-{\theta }_2)\widehat{m}+\gamma ,y_{n+3}=-{\theta }_1+{\theta }_2+{\theta }_3-{\theta }_4\hspace{0.17em}\hfill & \hfill \\ \hfill & \hfill & {\displaystyle \sum _{i+j=2l-1}{x}_{ij}}=0\hfill & \forall l\in [n+4]\hfill \\ \hfill & \hfill & {\displaystyle \sum _{i+j=2l}{x}_{ij}}={\displaystyle \sum _{q=0}^l{\displaystyle \sum _{r=q}^{n+3+q-l}{y}_r}}\left(\begin{array}{c}r\\ q\end{array}\right)\left(\begin{array}{c}n+3-r\\ l-q\end{array}\right)a^{r-q}{\widehat{m}}^q\hfill & \forall l\in [n+4]\cup \left\{0\right\}\hfill \\ \hfill & \hfill & z_1=R\mu -C-z_0+z_{n+3},\hspace{0.17em}{z}_2=-R\mu -z_{n+3},\hfill & \hfill \\ \hfill & \hfill & z_3,\dots ,z_n=0,\hspace{0.17em}{z}_{n+1}=-R\mu +C(n+1)-z_0,\hfill & \hfill \\ \hfill & \hfill & z_{n+2}=R\mu -Cn+z_0-z_{n+3}\hspace{0.17em}\hfill & \hfill \\ \hfill & \hfill & z_0=-({\theta }_1-{\theta }_2)\widehat{m}+\gamma ,z_{n+3}={\theta }_1-{\theta }_2+{\theta }_3-{\theta }_4\hspace{0.17em}\hfill & \hfill \\ \hfill & \hfill & {\displaystyle \sum _{i+j=2l-1}{x}_{ij}^{\prime }}=0\hfill & \forall l\in [n+4]\hfill \\ \hfill & \hfill & {\displaystyle \sum _{i+j=2l}{x}_{ij}^{\prime }}={\displaystyle \sum _{q=0}^l{\displaystyle \sum _{r=q}^{n+3+q-l}{y}_r^{\prime }}}\left(\begin{array}{c}r\\ q\end{array}\right)\left(\begin{array}{c}n+3-r\\ l-q\end{array}\right){\widehat{m}}^{r-q}{b}^q\hfill & \forall l\in [n+4]\cup \left\{0\right\}.\hfill \end{array}$$

Note that when d_l = d_u and m_l = m_u, setting $\alpha ={\theta }_1-{\theta }_2$ and $\beta ={\theta }_3-{\theta }_4$ recovers the dual Problem (9) in the view of the primitive MAD ambiguity set, which corresponds to the case when we have absolute trust on the mean and MAD estimators.

The next theorem presents the reformulation of the worst-case expected revenue rate.

Theorem 5.

For any $n\ge 1$, the worst-case expected revenue rate ${\overline{Z}}_r(n)$ is equal to the optimal value of the following linear problem:

$$\begin{array}{l}\underset{\theta \in {\mathbb{R}}_{+}^4,\gamma \in \mathbb{R}}{\mathit{\text{sup}}}\hspace{1em}\gamma +{\theta }_1{d}_l-{\theta }_2{d}_u+{\theta }_3{m}_l-{\theta }_4{m}_u\hfill \\ s.t.\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}({\theta }_1-{\theta }_2)|a-\widehat{m}|+({\theta }_3-{\theta }_4)a+\gamma \le r_n(a)\hfill \\ ({\theta }_3-{\theta }_4)\widehat{m}+\gamma \le r_n(\widehat{m})\hfill \\ ({\theta }_1-{\theta }_2)|b-\widehat{m}|+({\theta }_3-{\theta }_4)b+\gamma \le r_n(b).\hfill \end{array}$$

Theorems 4 and 5 provide tractable reformulations for the social and revenue optimization problems. An advantage of the proposed data-driven model is that it can offer attractive finite-sample guarantees. Compared with the original MAD ambiguity set that imposes unique mean and MAD, the data-driven MAD ambiguity set allows these parameters to vary within the confidence intervals. In this way, we can assure that the set contains the true underlying distribution with a high probability, which immediately generates out-of-sample performance guarantees for the solution.

Theorem 6.

Let ${\{{\widehat{\rho }}_i\}}_{i\in [N]}$ be a set of N samples generated independently at random from ${\mathbb{P}}^{\star }$ and $v^{\star }$ denote the optimal value of (14). Define $\widehat{m}$ and $\widehat{d}$ as the empirical mean and MAD obtained from samples ${\{{\widehat{\rho }}_i\}}_{i\in [N]}$. By setting

$$\begin{array}{ll}\hfill & \mathcal{T}=\left[\widehat{m}-(b-a)\sqrt{\frac{\log 4/\delta }{2N}},\widehat{m}+(b-a)\sqrt{\frac{\log 4/\delta }{2N}}\right]\hfill \\ \hfill & \mathcal{D}=\left[\widehat{d}-(b-a)\sqrt{\frac{9 \log 4/\delta }{2N}},\widehat{d}+(b-a)\sqrt{\frac{9 \log 4/\delta }{2N}}\right],\hfill \end{array}$$(16)

we have

$$\text{Prob}(v^{\star }\le {\mathbb{E}}_{{\mathbb{P}}^{\star }}[c_{\widehat{n}}(\tilde{\rho })])\ge 1-\delta ,$$

where $\widehat{n}$ is the optimal threshold obtained from (14).

Proof.

The error of the empirical MAD estimate is given by

$$\begin{array}{ll}\hfill & \left|\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-\widehat{m}|\right]\right|\hfill \\ =\hfill & \max \left\{\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-\widehat{m}|\right],-\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|+\mathbb{E}\left[|\tilde{\rho }-\widehat{m}|\right]\right\}.\hfill \end{array}$$

We upper bound both terms inside the max operator. The first term is bounded by

$$\begin{array}{ll}\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-\widehat{m}|\right]\hfill & \le \hspace{0.17em}\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[||\tilde{\rho }-m|-|\widehat{m}-m||\right]\hfill \\ \hfill & \le \hspace{0.17em}\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-m|-|\widehat{m}-m|\right]\hfill \\ \hfill & \le \hspace{0.17em}\left|\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-m|\right]\right|+|\widehat{m}-m|,\hfill \end{array}$$

where the second inequality follows from reverse triangle inequality. Meanwhile, the second term is bounded by

$$\begin{array}{l}-\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|+\mathbb{E}\left[|\tilde{\rho }-\widehat{m}|\right]\le \hspace{0.17em}-\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|+\mathbb{E}\left[|\tilde{\rho }-m|+|\widehat{m}-m|\right]\hfill \\ \le \hspace{0.17em}\left|\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-m|\right]\right|+|\widehat{m}-m|.\hfill \end{array}$$

Because both of these two terms have the same upper bound, we have

$$\left|\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-\widehat{m}|\right]\right|\le \left|\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-m|\right]\right|+|\widehat{m}-m|.$$

As $\mathbb{E}[\widehat{m}]=\mathbb{E}[m]$ is an unbiased estimator, we can invoke the Hoeffding inequality to directly derive a confidence interval for the second term. However, the empirical MAD is biased (i.e., $\mathbb{E}[\frac{1}{N}{\sum }_{i=1}^N|{\rho }_i-\widehat{m}|]\ne \mathbb{E}\left[|\rho -m|\right]$), making the Hoeffding inequality not applicable. To derive a confidence interval for this term, we rewrite it as

$$\begin{array}{l}\left|\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-m|\right]\right|\hfill \\ =\max \left\{\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-m|\right],-\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|+\mathbb{E}\left[|\tilde{\rho }-m|\right]\right\}.\hfill \end{array}$$

We further upper bound the two terms inside the max operator. For the first term, we have

$$\begin{array}{ll}\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-m|\right]\hfill & \le \hspace{0.17em}\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-m|+|m-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-m|\right]\hfill \\ \hfill & \le \hspace{0.17em}\left|\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-m|-\mathbb{E}\left[|\tilde{\rho }-m|\right]\right|+|m-\widehat{m}|.\hfill \end{array}$$

For the second term, applying the reverse triangle inequality yields

$$\begin{array}{ll}\mathbb{E}\left[|\tilde{\rho }-m|\right]-\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|\hfill & \le \hspace{0.17em}\mathbb{E}\left[|\tilde{\rho }-m|\right]-\frac{1}{N}{\displaystyle \sum _{i=1}^N||{\widehat{\rho }}_i-m|-|\widehat{m}-m||}\hfill \\ \hfill & \le \hspace{0.17em}\mathbb{E}\left[|\tilde{\rho }-m|\right]-\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-m|+|\widehat{m}-m|\hfill \\ \hfill & \le \hspace{0.17em}\left|\mathbb{E}\left[|\tilde{\rho }-m|\right]-\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-m|\right|+|\widehat{m}-m|.\hfill \end{array}$$

Thus, we have

$$\left|\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-\widehat{m}|-\mathbb{E}\left[|\tilde{\rho }-m|\right]\right|\le \left|\mathbb{E}\left[|\tilde{\rho }-m|\right]-\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-m|\right|+2|\widehat{m}-m|.$$

Because both of these two terms are unbiased, we can apply the Hoeffding inequality and obtain

$$\begin{array}{ll}\hfill & \text{Prob}\left(\left|\mathbb{E}\left[|\tilde{\rho }-m|\right]-\frac{1}{N}{\displaystyle \sum _{i=1}^N|}{\widehat{\rho }}_i-m|\right|\ge r_1\right)\le 2 \exp \left(-\frac{2N{r}_1^2}{{(b-a)}^2}\right)\hspace{1em}\text{and}\hfill \\ \hfill & \text{Prob}\left(|\widehat{m}-m|\ge r_2\right)\le 2 \exp \left(-\frac{2N{r}_2^2}{{(b-a)}^2}\right).\hfill \end{array}$$

By applying the union bound and setting $r_1=r_2=r/3$, we arrive at the desired confidence intervals that the true mean m and MAD d satisfy

$$\begin{array}{l}\widehat{m}-(b-a)\sqrt{\frac{\log 4/\delta }{2N}}\le m\le \widehat{m}+(b-a)\sqrt{\frac{\log 4/\delta }{2N}}\hfill \\ \widehat{d}-(b-a)\sqrt{\frac{9 \log 4/\delta }{2N}}\le d\le \widehat{d}+(b-a)\sqrt{\frac{9 \log 4/\delta }{2N}}\hfill \end{array}$$

with probability at least $1-\delta$. Therefore, by setting the confidence interval $\mathcal{T}$ and $\mathcal{D}$ as in (16), we have

$$\text{Prob}\left({\mathcal{P}}_N^{\prime }\ni {\mathbb{P}}^{\star }\right)\ge 1-\delta ,$$

where ${\mathcal{P}}_N^{\prime }$ is the data-driven ambiguity set (15) constructed by N random samples drawn from the underlying distribution ${\mathbb{P}}^{\star }$. Because $v^{\star }≔{\inf }_{\mathbb{P}\in {\mathcal{P}}_N^{\prime }}{\mathbb{E}}_{\mathbb{P}}[c_n(\tilde{\rho })]$ and the probability that ${\mathcal{P}}_N^{\prime }$ contains the true distribution ${\mathbb{P}}^{\star }$ is at least $1-\delta$, we have

$$\text{Prob}(v^{\star }\le {\mathbb{E}}_{{\mathbb{P}}^{\star }}[c_{\widehat{n}}(\tilde{\rho })])\ge 1-\delta ,$$

which completes the proof. □

The theorem establishes that with judicious choices of the confidence interval lengths, the optimal value of the data-driven DRO model $v^{\star }$ provides a high-confidence lower bound on the expected benefit rate of the robust solution $\widehat{n}$ under the true underlying distribution ${\mathbb{P}}^{\star }$.

Remark 2.

An avid reader may be interested in employing the popular Wasserstein DRO model in the data-driven setting. Indeed, the model has been widely adopted because it can generate asymptotically consistent solutions and offer similarly attractive finite-sample guarantees. Unfortunately, the reformulation of this data-driven DRO model involves $\mathcal{O}(N)$ semidefinite constraints, which make the problem computationally intensive. For readers who are interested in the use of the Wasserstein ambiguity set, we provide a detailed discussion in Appendix B.

5. Numerical Experiments

In this section, we present the numerical experiments and examine the performance of different DRO policies. All optimization problems are implemented in MATLAB and solved by SDPT3 (Toh et al. 1999) via the YALMIP interface (Lofberg 2004). The experiments are run on a 2.2-GHz Intel Core i7 CPU laptop with 8 GB RAM.

We assess the out-of-sample performance of the data-driven policies for a social optimizer and a revenue maximizer through a fair out-of-sample experiment. We assume we have access to N independent samples ${\{{\widehat{\rho }}_i\}}_{i\in [N]}$ of the traffic intensity drawn from the true underlying distribution ${\mathbb{P}}^{\star }$, and we construct four ambiguity sets: an empirical MAD ambiguity set, an empirical variance ambiguity set, a DD-MAD ambiguity set, and a Wasserstein ambiguity set. The empirical MAD ambiguity set is defined in (7), where we directly substitute the empirical mean and MAD for m and d, respectively. The empirical variance model is another popular moment model that constructs its ambiguity set based on the empirical mean and variance (i.e., $\mathbb{E}[{(\tilde{\rho }-m)}^2]={\sigma }^2$). Because its formulation and derivation parallel those of the empirical MAD model, we omit its discussion for brevity. The DD-MAD ambiguity set is defined in (15), where rather than carelessly plugging in the empirical estimators, we construct a confidence interval around the empirical mean and MAD. The Wasserstein ambiguity set (Esfahani and Kuhn 2018, Gao and Kleywegt 2023) is a popular data-driven ambiguity set. However, its complexity scales with the number of samples, making the problem computationally intensive with large sample sizes. We derive the reformulation of the Wasserstein model in Appendix D. Once we constructed the different ambiguity sets, we then proceed to compute the distributionally robust thresholds that maximize the respective worst-case expected benefit rates. Finally, we compare the three solutions in a fair out-of-sample experiment relative to the sample average approximation (SAA) method, which naively assumes that the empirical distribution generated from the N samples is the true underlying distribution. The SAA method also represents the stochastic model (Liu and Hasenbein 2019) under the empirical distribution.

We conduct the out-of-sample trials for data sets containing $N=2,4,\dots ,10,20,40\dots ,100$ independent samples. We assume the arrival rate is generated by $\lambda =4\tilde{b}$, where $\tilde{b}\sim \text{Beta}(0.1,0.5)$. In addition, we assume the experienced decision maker has an educated guess for the distribution support as $[0,5]$. In each trial, we draw N independent training samples and obtain ${\{{\widehat{\rho }}_i\}}_{i\in [N]}$ from ${\mathbb{P}}^{\star }$. We then compute the optimal thresholds ${\widehat{n}}_d$, n_v,${\widehat{n}}_{dd}$, and ${\widehat{n}}_w$ for the MAD, variance, DD-MAD, and Wasserstein DRO models, respectively. We also compute the SAA threshold ${\widehat{n}}_{\mathit{\text{SAA}}}$ by solving the sample average approximation model. Based on the scaling rates derived in Theorem 6 and Esfahani and Kuhn (2018, theorem 3.4), the size of the confidence intervals in (14) is set to be $C_1/\sqrt{N}$, and the Wasserstein radius is set to be $C_2/\sqrt{N}$, where C₁ and C₂ are chosen from the set $\{5,1,0.5,0.1,0.05,0.01\}$ using a $k\mathit{\text{-fold}}\hspace{0.17em}\mathit{\text{cross}}\hspace{0.17em}\mathit{\text{validation}}$ procedure. Specifically, we partition the in-sample data ${\{{\widehat{\rho }}_i\}}_{i\in [N]}$ into $k=\min \{N,5\}$ folds and repeat the following procedure for each fold; the ith fold is taken as a validation data set, and the remaining k − 1 folds are merged to be a subtraining set. We repeat this process for each fold and choose the interval length that performs best in average. The reason why we do not directly plug in the theoretical values from Theorem 6 is that the bound holds for any underlying distributions, which can be overly conservative in practice. Finally, the out-of-sample expected benefit/revenue rate ${\mathbb{E}}_{{\mathbb{P}}^{\star }}[c_{\widehat{n}}(\tilde{\rho })]$ for each of the strategies is then estimated at high accuracy using 10,000 test samples from ${\mathbb{P}}^{\star }$.

Figure 2 depicts the out-of-sample performances of a social optimizer and a revenue optimizer under different DRO policies with R = 10, C = 1, and μ = 1. The expected values and 95 percentiles are computed from 50 independent trials. The y axis represents the improvements of the DRO policies relative to the SAA policy, whereas the x axis denotes the sample size. In the social optimization problem, the curve of the Wasserstein model terminates at n = 6 because the solver fails to converge when the sample size reaches eight. Meanwhile, the Wasserstein model dominates the SAA model uniformly across all sample sizes in the revenue maximization problem, whereas the MAD, variance, and DD-MAD models outperform the SAA model for small to medium sample sizes. This is because the Wasserstein ambiguity set converges to the true distribution as the number of samples grows, whereas the moment ambiguity sets fail to converge to the true distribution. We also find that the MAD model performs poorly when the sample size is small because the empirical MAD constitutes a biased estimator with significant estimation errors. On the other hand, the DD-MAD model—by optimizing in view of the most adverse mean and MAD—mitigates the detrimental effects of poor empirical estimations and generates high-quality policies. Moreover, we notice that the empirical MAD and variance models yield similar scores, suggesting that measuring dispersion by MAD or variance does not influence the performance of the model. Finally, we observe that the advantages of the DRO policies relative to the SAA method are generally more substantial in terms of the 95th percentiles of improvements. This underlines a major advantage of incorporating the DRO scheme as it reduces the likelihood of realizing inferior performance in the out-of-sample tests.

Table 1 reports the computation time of different models with the sample sizes varying from 2 to 100. We set the length of the confidence intervals and the radius of the Wasserstein ball to 0.1. In this experiment, the running time limit of SDPT3 is set to 600 seconds, and the number of iterations is set to 5,000. All computational times are averaged over 10 trials.

Table 1. Running Time (in Seconds) of Different Methods

The results in Table 1 indicate that the computational times of the MAD, variance, and DD-MAD models are size invariant because the number of constraints is independent of the number of samples. For the social optimization problem, the Wasserstein model is applicable to small-size instances. However, it encounters computational difficulties for moderate-size problem instances; when the sample size reaches 10, the SDP solver fails to converge within the time/iteration limit. Benefiting from the closed-form solution, the MAD model is more efficient than the variance and DD-MAD models in both the social and revenue optimization problems. For the revenue optimization problem, the DD-MAD and Wasserstein models admit a linear programming (LP) reformulation, whereas the variance model still leads to an SDP reformulation. This underlines a major advantage of using MAD at the moment ambiguity set as it can significantly improve the model’s efficiency.

The DD-MAD model is still size invariant, and its linear programming reformulation yields a much shorter computational time than the SDP reformulation for the social optimization problem. In addition, the Wasserstein model can be solved efficiently even for large sample sizes, benefiting from the linear programming reformulation.

Finally, we report the performance of the DRO models under different scaling parameters C in Figure 3. The expected improvements are computed with n = 5 samples from 50 independent trials. We observe that both models have large variations in performance with different scaling parameters. The Wasserstein and DD-MAD models yield unimodal curves in both the social and revenue optimization problems, implying a trade-off between performance and conservatism. Intuitively, including robustness can improve the out-of-sample performance, whereas being too conservative may also adversely affect the results. To achieve the best performance, one could set the size of the ambiguity set or confidence interval to the best radius. Unfortunately, we do not have access to this information. Although one can plug in the theoretical values obtained from concentration inequalities, these values are usually too conservative. In practice, decision makers can rely on a crossvalidation or bootstrap procedure to obtain a suitable size for the ambiguity set (Gotoh et al. 2021, Bates et al. 2023). From the figure, we further observe that the two models perform quite differently when the scaling parameter C is small; the DD-MAD model achieves significant improvement, whereas the Wasserstein model only yields a slight improvement. The reason is that the Wasserstein ambiguity set is centered on the empirical distribution. When the scaling parameter is small, all the distributions within this ambiguity set are close to the empirical distribution. Thus, it generates similar results to SAA and cannot achieve a large improvement. Conversely, the DD-MAD model converges to the empirical MAD model when C = 0. As illustrated in the first set of experiments, the MAD model outperforms the SAA method for small sample sizes. Hence, the DD-MAD method yields a substantial improvement when the scaling parameter is small.

6. Conclusion

This paper developed an extension of the Naor (1969) strategic queue model with uncertain arrival rates using the DRO framework. We showed that under the DRO setting, the optimal threshold of an individual optimizer coincides with the original result of Naor (1969), and there exist optimal thresholds of the social and revenue optimizers not larger than the optimal individual threshold. We then proved that the revenue rate function is concave, whereas the social benefit rate function is concave or unimodal under some mild conditions. These nice properties lead to a closed-form solution for the revenue maximization problem and an analytical solution for the social optimization problem.

Next, we considered the data-driven optimization setting, where decision makers only have access to limited historical samples. We proposed a data-driven MAD model by introducing an extra layer of robustness to the primitive MAD ambiguity set. As the model mitigates the detrimental estimation errors from the empirical mean and MAD, it achieves attractive performance in out-of-sample tests. We derived an SDP reformulation for the social optimization problem and a linear programming reformulation for the revenue maximization problem. We further established finite-sample guarantees for the data-driven model, which provide valuable guidance for choosing the robustness parameters in practice. Our experimental results demonstrate that a system manager who disregards ambiguities in the arrival rate distribution as well as errors from the empirical parameter estimations may incur large out-of-sample costs. Future work includes extending the DRO scheme to the unobservable strategic queues, where newly arrived customers cannot observe the current length of the queue system.