# Stochastic Approximation of Symmetric Nash Equilibria in Queueing Games

Liron Ravner  Department of Statistics, University of Haifa Ran I. Snitkovsky
Graduate School of Business, Columbia University  School of Data Science,
Shenzhen Research Institute for Big Data,

###### Abstract

We suggest a novel stochastic-approximation algorithm to compute a symmetric
Nash-equilibrium strategy in a general queueing game with a finite action
space. The algorithm involves a single simulation of the queueing process with
dynamic updating of the strategy at regeneration times. Under mild assumptions
on the utility function and on the regenerative structure of the queueing
process, the algorithm converges to a symmetric equilibrium strategy almost
surely. This yields a powerful tool that can be used to approximate
equilibrium strategies in a broad range of strategic queueing models in which
direct analysis is impracticable.

##  1 Introduction

This paper presents a simulation-based algorithm that computes, or learns, a
symmetric Nash equilibrium in a general class of queueing games. Our purpose
is to suggest a robust, easy-to-implement equilibrium approximation scheme,
mainly targeted at queueing games in which identifying equilibrium strategies
is difficult to pursue analytically. In particular, we consider service
systems to which rational customers arrive according to a renewal process,
each chooses between finitely many actions with the goal of maximizing their
expected utility. Customers in our model assume that the system operates in a
steady state, however, their assessment of the utility depends on other
customers’ actions, giving rise to a game-theoretic equilibrium problem:
Identifying a behavior profile that prescribes an optimal play for customers,
such that the induced steady state is consistent with each customer’s belief
regarding the behavior of the others.

Though the initiation of the literature about strategic customer behavior in
queues is often accredited to Naor [59], the most basic, yet non-degenerated
queueing game was first introduced by Edelson and Hildebrand [26], and was
coined the unobservable M/M/1 model. It discusses a single-server queue were
customer utility decreases linearly with the waiting time, and customers need
to decide whether to join the queue or balk, without observing the system
state. Since then, various different queueing games have been studied in the
literature, to model the behavior of callers in call centers, users of
communication networks, passengers in a transportation system, and more (see
overviews and exhaustive surveys in [37] and [34]).

Traditionally, this literature distinguishes between ‘unobservable’ and
‘observable’ models, though this conceptual distinction is rather obscure. By
unobservable, one usually means that the strategy of a customer cannot rely on
any state information but the premise that this state is drawn from the
stationary distribution. Conversely, observable (and partially-observable)
models are commonly used to describe games in which customers are endowed with
information about the state to which they arrive, hence, making their strategy
state dependent. In that sense, the observable framework arguably allows a
richer description of customer behavior, although technically, many observable
queueing games discussed in the literature can be reframed as unobservable
models via a suitable modification of the action set.

In most queueing games, observable or unobservable, customer expected utility
(ex ante) is a function of their own strategy, in which the expectation is
with respect to the stationary distribution arising from the cumulative
behavior of all players. When customers are homogeneous, the conventional
solution concept is the symmetric Nash equilibrium. Characterizing a symmetric
Nash equilibrium means finding a fixed point of the best-response function,
which is the function that maps a strategy profile to the best strategy to
play against it.

Unlike the Markovian setting in [26], for non-elementary queueing models, the
stationary state distribution, and hence also the expected utility, rarely
admit closed form expressions. Think, for example, of a GI/G/1 queue, with
customer utility depending on their delay. In this example, as well as in many
others, there is no explicit formula for the expected sojourn time of a
customer (let alone its distribution), hence, approaching an equilibrium
solution with only analytic tools is intractable. Tackling such problems often
brings the need for sophisticated numeric calculation and simulation schemes.

We present a stochastic-approximation (SA) algorithm that converges to an
equilibrium solution for a general class of models. Specifically, the
algorithm involves simulating the queueing process and updating the strategy
based on the realizations at carefully chosen regeneration times. The update
relies on estimating the deviation of the strategy from a best response, after
a smoothing transformation of the best-response function. Under mild
regularity assumptions it is shown that the algorithm converges almost surely
to a Nash equilibrium as the simulation length goes to infinity. The
convergence conditions are verified for several examples, and numerical
results are presented.

Simulation is a very popular technique among queueing theorists and
practitioners as a tool to approximate the performance of a model that cannot
be approached analytically. The literature studying the theory and
applications of simulation methods for queueing networks is extensive (see [7]
for a comprehensive overview), and algorithms can efficiently learn the
performance of quite intricate networks, given the set of primitives. However,
in the strategic queueing domain, simulation is hardly ever used, primarily
due to the difficulties that arise in the equilibrium search process: In
queueing games, the strategy profile governs the dynamics of the underlying
system. Thus, to identify its best response (or an approximation thereof), the
strategy should be fed to the simulator up front. For a strategy given a
priori, simulation can be used to verify, with high certainty, whether this
strategy approximately meets the equilibrium condition. But when the
equilibrium strategy is to be found, with no simple solution at hand, it
usually requires traversing through the set of possible strategies, and
performing a separate simulation at each iteration to check each strategy.
This procedure can be tedious and even impractical, especially when the
simulation time at each iteration is exceedingly long. Moreover, the outcome
of each iteration is subject to uncertainty, and so, if not carefully
implemented, such a process is not guaranteed to converge at all, or perhaps
worse – converge to an undesired limit.

It is possible that for these reasons the literature of strategic queueing is
often limited to stylized models, of which performance measures can be easily
expressed. In their own right, stylized models play an important role in
understanding the theoretical properties of queueing games. Yet, for the
purpose of departing from a merely theoretical framework, it is crucial to
also have the ability to compute the equilibrium outcome in elaborate systems.

We introduce a novel simulation scheme that jointly learns the system’s
statistical characteristics and customer best-response dynamics, such that the
process efficiently converges to a Nash equilibrium. Our framework allows for
customer strategies to depend on state information, thus, the results apply
for a wide range of applications. Still, the key results and intuition
supporting these results are more easily conveyed using the so-called
‘unobservable’ model. We first discuss our method in the context of
unobservable queueing games, and later complement the discussion by suggesting
a refinement for observable models as well. We further discuss several
practical and theoretical aspects of implementation.

Our method relies on the SA algorithm commonly known as the Robbins-Monro
algorithm [62], a fundamental building block underpinning various stochastic
optimization techniques, among which is the renowned Stochastic Gradient
Descent method. However, as opposed to standard optimization problems, we do
not seek the root of a gradient. In our framework, the equilibrium condition
needs to satisfy an indifference principle between randomly-chosen actions.
Under regularity conditions, an equilibrium is found by the SA algorithm in
the limit with probability one. In other words, a long enough (single)
simulation is guaranteed to converge to an equilibrium strategy. This idea
draws inspiration from our understanding of how equilibrium emerges in real
world systems: customers adapt their strategy based on observing past
empirical performance of the system, and eventually converge to an
equilibrium.

###  1.1 A motivating example – Two unobservable GI/G/1 queues in parallel

The main purpose of the paper is to describe a general solution technique that
can be applied to a broad range of applications, and therefore the fundamental
queueing game formulated in Section 2, and its extension in Section 5, are
relatively abstract. Yet to demonstrate the competency of our method we first
look at a special setup that is easily explained on one hand, but is
analytically intractable on the other hand.

Consider a network, termed the _system_ , consisting of two FCFS queues
(stations), indexed by $m\in\\{1,2\\}$, with dedicated servers working in
parallel, one for each queue. Service times at each station $m=1,2$ are
independent and follow a general light-tailed distribution $F_{m}$ with mean
$1/\mu_{m}$, such that w.l.o.g. $1/\mu_{1}\geq 1/\mu_{2}$. Potential customers
arrive at the system according to a renewal process with inter-arrival time
distribution $H$ and mean $1/\lambda$. The state of the system at an arbitrary
point in time is characterized by a vector $X=(X^{[1]},X^{[2]})$, with
$X^{[m]}$ being the workload at station $m$.

Upon arrival, potential customers make decisions based on a prior belief
regarding the mean waiting time at each queue (e.g., relying on historic
observations), yet, they cannot observe the current system state. Each
customer strategically chooses one of three possible actions: (1) Join station
$1$; (2) Join station $2$; or (3) Balk. Respectively, we denote the set of
actions for a customer by $\mathcal{A}=\\{1,2,3\\}$. A strategy
$\mathbf{p}=(p_{1},p_{2},p_{3})\in\Delta(\mathcal{A})$ is a distribution over
the action set with $p_{i}$ being the probability of taking action
$i\in\mathcal{A}$. We assume the utility from joining is linear in the waiting
time: Following the conventional notation in [37], let $R>0$ denote the
customer reward for service and $C>0$ be the customer waiting time cost. When
a customer finds the system at state $X$ and joins queue $m\in\\{1,2\\}$,
their utility is given by $v_{m}=R-C\cdot(X^{[m]}+Y_{m})$, where $Y_{m}\sim
F_{m}$ is a r.v. representing the customer’s service time at server $m$. The
utility from balking is normalized to 0, i.e., $v_{3}=0$.

Customers assume when they arrive, that the state $X$ admits its stationary
distribution (embedded at arrival instants), which is determined by the
population strategy $\mathbf{p}$. To avoid diverting the discussion from its
main purpose towards dealing with stability issues we assume that
$\lambda<\mu_{1}$, implying that $X$ indeed admits a stationary distribution
for every $\mathbf{p}$. A detailed discussion of the stability conditions for
generalizations of this model can be found in Section 4.1.

Equilibrium strategy. We are interested in characterizing an equilibrium
strategy
$\mathbf{p}^{e}=(p^{e}_{1},p^{e}_{2},p^{e}_{3})\in\Delta(\mathcal{A})$, such
that under the stationary state distribution induced by $\mathbf{p}^{e}$, each
customer’s choice is utility maximizing in expectation. In other words,
$\mathbf{p}^{e}$ is a strategy such that for each $i\in\mathcal{A}$

$p^{e}_{i}>0\Rightarrow
i\in\operatorname*{arg\,max}_{j\in\mathcal{A}}\mathrm{E}[v_{j}],$

where the expectation is taken jointly w.r.t the random variables
$Y_{1},Y_{2}$ and the steady-state distribution of $X$ induced by
$\mathbf{p}^{e}$. The following lemma, whose proof follows standard arguments
(see Section B.1 in the appendix), establishes existence and uniqueness of the
targeted equilibrium strategy.

######  Lemma 1.

Assume the inter-arrival distribution $H$ is continuous. For any pair of
distributions $(F_{1},F_{2})$ and parameters $R$ and $C$, there exists a
unique symmetric Nash equilibrium strategy $\mathbf{p}^{e}$.

However, except for some special cases, an explicit characterization of this
equilibrium strategy is not available, because the values of
$\mathrm{E}[v_{1}]$ and $\mathrm{E}[v_{2}]$ given strategy $\mathbf{p}$ are
inaccessible. Whether a strategy
$\mathbf{p}=(p_{1},p_{2},p_{3})\in\Delta(\mathcal{A})$ satisfies the
equilibrium condition is a question that cannot be answered accurately, but
rather approximately, based on one’s ability to approximate
$\mathrm{E}[v_{1}]$ and $\mathrm{E}[v_{2}]$ for the underlying strategy
$\mathbf{p}$. For example, one can determine whether $\mathbf{p}$ satisfies an
$\varepsilon$-equilibrium condition, with $\varepsilon$ depending on the
ability to bound the approximation error for $\mathrm{E}[v_{j}]$.

Stochastic approximation. For a given strategy $\mathbf{p}$, assume a customer
arrives at the system when both queues are empty. Then we can consider this
moment of arrival as an instant of system regeneration. Our method relies on
simulating regeneration cycles, i.e., simulating the queueing process between
two successive instants of system regeneration. During each cycle we assume
that the customer strategy $\mathbf{p}$ is fixed, and keep track of the system
state at every arrival instant. Letting the r.v. $L$ denote the number of
arrivals during a cycle (including the arrival that initiates the cycle), we
therefore generate $L\geq 1$ samples of the system state, $X_{1},\dots,X_{L}$,
with $X_{j}=(X_{j}^{[1]},X_{j}^{[2]})$, and combine them together by defining

$\mathbf{G}=\begin{pmatrix}G_{1}\\\ G_{2}\\\
G_{3}\end{pmatrix}=\sum_{j=1}^{L}\begin{pmatrix}R-C\cdot(X_{j}^{[1]}+1/\mu_{1})\\\
R-C\cdot(X_{j}^{[2]}+1/\mu_{2})\\\ 0\end{pmatrix}.$

In words, we sum up the expected utility over all arriving customers, for each
action they could have taken, had they observed the system state upon arrival.
The motivation for the construction of $\mathbf{G}$ above is that for any
$\mathbf{p}$, the term $G_{i}/\mathrm{E}[L]$, $i\in\\{1,2,3\\}$, can be
regarded as an unbiased point estimator for $\mathrm{E}[v_{i}]$. Yet the
values $\mathrm{E}[L]$ and $\mathrm{E}[v_{i}]$ themselves are assumed unknown,
and obviously are not given as inputs to the algorithm. This particular form
of the estimator $\mathbf{G}$ is tailored to correct for the length bias that
often arises when estimating performance measures of a stochastic processes
from random samples.

An application of our SA algorithm works as follows: Setting an arbitrary
initial strategy $\mathbf{p}^{(1)}$, at each iteration $n=1,2,\dots$, we
construct an estimator $\mathbf{G}^{(n)}$ assuming the strategy
$\mathbf{p}^{(n)}$ is fixed. Given a constant $\gamma_{0}>0$, we then apply
the following update rule:

$\mathbf{p}^{(n+1)}=\pi_{\Delta}\left(\mathbf{p}^{(n)}+\frac{\gamma_{0}}{n}\mathbf{G}^{(n)}\right),$

where $\pi_{\Delta}$ denotes the projection onto the standard (2-)simplex. In
Section 3.2 of the paper we prove that under mild regularity of the primitives
the proposed scheme converges (almost surely) to the true equilibrium,
$\mathbf{p}^{(n)}\>\to_{\rm as}\>\mathbf{p}^{e}$. The parameter $\gamma_{0}$
is measured in the reciprocal of utility units, and clearly, the choice of its
value will impact the performance of the algorithm, but we defer the
discussion of this issue to later sections.

Simulation results. Below are results from a simulation, in which $F_{1}$ is
$\rm{Beta}(10,10)+0.5$, $F_{2}$ is $\rm{Bernoulli}(0.1)\cdot 10$, and inter-
arrivals are distributed according to $\rm{Gamma}(0.1,11)$, implying
altogether that $\lambda=10/11<\mu_{1}=\mu_{2}=1$. The reward and cost
parameters are given by $R=5$ and $C=1$. To set up the simulation we
initialize $\mathbf{p}^{(1)}=(1/3,1/3,1/3)$ and $\gamma_{0}$=0.1, and run the
algorithm for $N=10^{6}$ iterations. The convergence of the sequence
$\\{\mathbf{p}^{(n)}\\}$ is depicted in Figure 1. After terminating at
iteration $N=10^{6}$, the algorithm produces an output
$\mathbf{p}^{(N)}=(0.525,0.330,0.145)$. We note that the coefficient of
variation of $F_{2}$ is significantly higher than that of $F_{1}$ ($3$
compared to $0.11$), and generally speaking, higher service-time variation
leads to longer waiting time and in turn to smaller utility. Thus, it is
expected that in equilibrium, the negative impact of the higher variance in
Server 2 will be effectively balanced by less arrivals, namely, that
$p^{e}_{1}>p^{e}_{2}$.

Figure 1: A stacked-area plot showing the convergence of the SA algorithm as
$n\to\infty$. The coordinates of
$\mathbf{p}^{(n)}=({p}^{(n)}_{1},{p}^{(n)}_{2},{p}^{(n)}_{3})$ are plotted vs.
$n$ on a square-root scale. The blue curve depicts ${p}^{(n)}_{1}$, and the
green curve depicts ${p}^{(n)}_{1}+{p}^{(n)}_{2}$, thus, the blue, green and
orange shaded areas correspond with ${p}^{(n)}_{1}$, ${p}^{(n)}_{2}$, and
${p}^{(n)}_{3}$, respectively.

Verification of approximate equilibrium. For validation purpose, we estimate
the mean virtual workloads for the two queues by simulating $2\times 10^{6}$
arrivals to each of the two GI/G/1 queues separately (assuming they are
initially empty). The choice of simulation length follows the framework laid
out in [74]. Note that traffic intensities are relatively low in these two
cases, and (recalling that $\mu_{1}=\mu_{2}=1$) are equal to $\lambda
p^{(N)}_{1}\approx 0.477$ in Server 1, and $\lambda p^{(N)}_{2}\approx 0.3$ in
Server 2. In the two cases below, inter-arrival times are distributed as
$S=\sum_{i=1}^{Z}U_{i}$ where $U_{i}\sim\rm{Gamma}(0.1,11)$, $Z$ is
geometrically distributed with a different parameter at each case, and
services are distributed as described:

  1. 1.

$Z\sim{\rm Geom}(p^{(N)}_{1})$ and services are distributed according to
$F_{1}$ (which is $\rm{Beta}(10,10)+0.5$);

  2. 2.

$Z\sim{\rm Geom}(p^{(N)}_{2})$ and services are distributed according to
$F_{2}$ (which is $\rm{Bernoulli}(0.1)\cdot 10$).

Under the algorithm’s output strategy $\mathbf{p}^{(N)}$, we obtain
(consistent) estimators for $\mathrm{E}[X^{[1]}]$ and $\mathrm{E}[X^{[2]}]$,
in Cases 1 and 2, respectively, by averaging the observed workloads at arrival
instants. To construct asymptotic confidence intervals for the mean workloads
we follow the CLT-based framework of [7, Ch. III]: We set the significance
level $\alpha=0.005$. Under Case 1 we obtain an average virtual workload equal
to 4.006 with an asymptotic $1-\alpha$ confidence interval $[3.996,4.017]$.
Under Case 2 we get a sample average equal to 4.005 with an asymptotic
$1-\alpha$ confidence interval $[3.989,4.021]$. This suggests that with more
than $(1-\alpha)^{2}(\approx 0.99)$ certainty, our algorithm’s output
$\mathbf{p}^{(N)}$ is such that
$(\mathrm{E}[v_{1}],\mathrm{E}[v_{2}])\in\mathcal{I}_{1}\times\mathcal{I}_{2}$,
where $\mathcal{I}_{1}=[-0.017,0.004]$ and $\mathcal{I}_{2}=[-0.021,0.011]$.
Hence, with high certainty (>99%) our approximated solution $\mathbf{p}^{(N)}$
satisfies the criterion for an $\epsilon$-approximate Nash equilibrium ([23])
for
$\varepsilon\leq\max_{(x,y)\in\mathcal{I}_{1}\times\mathcal{I}_{2}}|x-y|=0.028$.

###  1.2 Outline and main results

As discussed, our main result is the derivation of a simulation-based
approximation technique to compute equilibria. The method is robust in the
sense that it can be applied to a wide range of strategic queueing models.
Below we summarize the organization of the paper and its results:

Theoretical results. In Section 2 we introduce a unified framework for
unobservable queueing games with homogeneous customers and finitely many
actions. As a byproduct we derive sufficient conditions for the existence of a
symmetric Nash equilibrium in this class of games (Lemma 3), which we show how
to reduce down to a fixed-point problem of a real continuous function. This
naturally yields an iterative (deterministic) mechanism to approximate
equilibrium solutions, had customer stationary utility been known. In Section
3 we explain how this stationary utility can be estimated based on samples
drawn from simulated regeneration cycles. The proposed estimator serves in
place of the unknown original utility, in combination with the iterative
fixed-point scheme from Section 2, to form the SA algorithm. The SA algorithm
is shown, in Theorem 5, to converge almost surely to an equilibrium under mild
regularity conditions. Our main convergence result is accompanied with a
thorough discussion of the assumptions and convergence-rate results in
Sections 3.3 and 3.4. In Section 5 we then extend the method to account for
state-dependent strategies in order to deal with observable models.

Examples and implementation. All along Sections 2 and 3 we use the canonical
unobservable M/G/1 queue as an instructive tool to exemplify and generate
intuition from the general model. In Section 4 we discuss relevant
applications, namely multiple GI/G/1 queues in parallel (Section 4.1), and
selective routing between queues with different buffer capacity (Section 4.2).
We verify the sufficient convergence conditions for these applications, and
use them to present refinements aimed at improving algorithm’s run time. In
Section 5 we discuss the observable GI/G/1 model as an example of observable
queueing games and provide results from numeric experiments.

###  1.3 Related literature

Fundamentals of queueing games. The game of joining or balking from an
unobservable queue was first introduced by Edelson and Hilderbrand [26]. They
consider an M/M/1 queue with strategic customers wishing to maximize a utility
function comprised of a fixed reward and a linear cost for waiting. As a
simple formula is readily available for the expected waiting time in an M/M/1
queue, the equilibrium joining probability can be computed explicitly. This
model easily extends to general service times [9]. Numerous other queueing
games have been studied over the years and we refer interested readers to the
books of Hassin and Haviv [37] and Hassin [34] for complete overview and
surveys. Recently, Haviv and Oz [42] suggest a unified approach to formulate
unobservable queueing games by analyzing the interactions between the (random)
set of customers who visit the system during the same busy period. Their
formulation is related to our work in that our SA algorithm updates the
strategy in between busy periods. Broadly speaking, both works rely on the
understanding that, given the strategy, all the statistical properties of the
underlying queueing process are encapsulated in a single regeneration cycle.

Stochastic approximations in queueing literature. A comprehensive overview of
the theory of stochastic approximations, which was initiated by the seminal
work of Robbins and Monro [62], can be found in Kushner and Yin [49] and in
Borkar [12]. An adaptation of the Robbins-Monro algorithm for optimization
purposes was first presented by Kiefer and Wolfowitz [48], who laid the
groundwork for the rapid development of stochastic optimization techniques.
Among their notably-wide range of applications, these techniques have been
applied to problems of optimizing the steady-state performance of queueing
systems, to which we give emphasis next.

The computation of an optimal service rate for a single-server queue using
stochastic approximations was initially studied by Suri and Zazanis [68] where
the gradient estimation problem is addressed via the method of Infinite
Perturbation Analysis. Empirical experiments for the M/M/1 queue (building on
the machinery developed in [68]) were carried by Suri and Leung [69]. Fu [27]
later extended this theory to the GI/G/1 queue, exploiting the regenerative
structure underlying it. Further convergence properties are given for the
GI/G/1 model in L’Ecuyer and Glynn [50], for different estimators of the
gradient such as finite-difference and likelihood-ratio estimators. A
numerical study of the techniques discussed in [50] for the M/M/1 case can be
found in L’Ecuyer et. al [51]. A review of optimization techniques using
simulation including applications to GI/G/1 queues is provided by Fu [28], and
a unified framework for stochastic optimization of the steady-state
performance in regenerative systems is established in Andradóttir [5]. In all
these settings the objective function is continuously differentiable in the
parameter and the focus is on estimating the gradient towards the ultimate
goal of computing its root. Our setting is different in that we do not focus
on optimizing a given objective, but rather on identifying Nash equilibrium
solutions – a problem that is not dealt with in this line of literature.

Recently, Chen et. al [18] proposed an SA algorithm for jointly optimizing
price and capacity in a GI/G/1 queue, and formulated conditions for
convergence to the optimal solution. Chen et. al [18] further perform thorough
analysis of the regret, drawing connections between stochastic approximations
and reinforcement learning in the context of revenue management in queues.
Reinforcement learning has been employed in recent studies to approach optimal
control policies as a means of optimizing steady-state performance in queueing
networks, by Liu et. al [52] and by Dai and Gluzman [21].

Simulation and learning in queueing games. In the area of strategic queueing,
literature involving simulation-based methods to approach equilibrium
solutions is scarce. Altman and Shimkin [3] explicitly analyze an observable
processor-sharing queue with strategic customers, exponential services, and
linear waiting cost. They suggest an ad-hoc simulation-based learning
algorithm, to demonstrate the convergence of customers’ iterative decision-
making process to a symmetric Nash equilibrium. This was extended to a model
with heterogeneous customers in Ben Shahar et al. [10]. Buche and Kushner [13]
verify convergence of the algorithm presented in [3] to a Nash equilibrium and
explain how it extends to non-linear cost functions. In the context of
customers strategically timing their arrivals to a transient queueing
processes, Sakuma et. al [65] construct a heuristic dynamic algorithm to
approximate the Nash equilibrium arrival strategy. We note, however, that the
framework of strategic arrival scheduling in queues is significantly different
than ours – in the former, the focus is on transient queueing processes, and
additionally, customer action space is usually assumed continuous (see more
details in [43]).

Indirectly related to our work is the study of best-response dynamics;
iterative updating of strategies by maximizing utility given the previous
strategy in the iteration. Best-response algorithms are known to converge in
some queueing games, such as S-modular games as defined by Yao [77]. For
example, the S-modular framework was applied to a problem of decentralized
control of a wireless network by Altman and Altman [2]. However, in most
settings the known conditions required for convergence of best-response
dynamics are not satisfied, and typically cannot even be tested due to the
intractability of the stationary performance measures of the system.

Potential applications. In the queueing game we introduce, customers are
modeled as short-lived entities who arrive according to a renewal process to a
general, regenerative system, and choose one out of finitely-many actions,
with the objective of maximizing their utilities. Ever since Naor’s [59]
seminal work, this modeling structure has been accepted as a standard approach
to study customer behavior in service systems, and has been applied widely in
the operations-management literature. While in [59] the existence of a
dominant strategy deems the equilibrium analysis trivial, here we focus our
attention on cases where the equilibrium search problem is non-degenerated.
Below is a noninclusive list of few representative examples of such models.

Among the classical problem themes in the study of strategic queueing are
join-or-balk decisions in observable and unobservable queues [26, 15, 25, 31,
47]; decentralized selective routing in queueing networks [9, 33, 61, 38];
paying to reduce wait via priority or premium service classes [1, 35, 72, 20];
and provision and acquisition of information in queues [32, 75, 40, 39, 41].
Concrete applications that motivate the study of such models traditionally
include communication and computer information services [57, 3, 54, 46] and
transportation networks [56, 55]. Some very recent business applications that
are increasingly emerging in this line of literature also extend across ride-
sharing platforms [70, 45, 44, 16]; food delivery and curbside pickup in
restaurants [17, 67]; as well as omnichannel services [8, 63, 30, 16].

In many of these papers, solving for customer equilibrium is a preliminary,
yet a crucial step towards answering questions concerning pricing, policy
selection, capacity planning and the like. Due to the intricacy of the
equilibrium problem, the majority of these examples aim to capture the main
characteristics of the equilibrium behavior through parsimonious models, on
which direct analysis can be carried out. Two exceptions here are [75], who
adopt a mean-field approach, and [54] who resort to heavy-traffic
approximations (in the Halfin-Whitt regime). In [54], the authors further
highlight the need for exhaustive simulations in tackling the pre-limit analog
of their model. The framework we study is fairly general, covers both
observable and unobservable games, and is easily extended to multiple
(finitely-many) customer types. It is rich enough to capture the customer-
behavior model employed in the examples listed above (with the restriction to
finite heterogeneity types when considered). Thus, our results potentially
offer a useful method to approach many of these examples’ variants and
extensions that are not amenable to direct analysis.

###  1.4 Notation

In the paper, the domain of interest over which we define our problem is
$\mathbb{R}^{k}$, $k>1$. Thus, we designate vectors in $\mathbb{R}^{k}$ by
bold letters, and so we do for functions whose image is in $\mathbb{R}^{k}$.
We denote by $\mathbf{e}_{i},i=1,\dots,k$, the standard-basis unit vector with
$1$ in its $i$-th coordinate, and denote by $\mathbf{e}$ the all-1
$k$-dimensional vector. By default, we define vectors as column vectors and
use ′ to denote their transpose. The non-negative real half line is denoted by
$\mathbb{R}_{+}$. Given a function $g$ with a domain in $\mathbb{R}$ we denote
its limit from the left at $t$ by $g(t-)=\lim_{s\uparrow t}g(s)$. For a real
vector space, $\|\cdot\|$ is used by default to refer to the $L^{2}$-norm, and
for $p\in[1,\infty)\cup\\{\infty\\}$, $\|\cdot\|_{p}$ denotes the
$L^{p}$-norm. We further use $\|\cdot\|_{0}$ to denote the $L^{0}$-“norm”,
i.e., the number of non-zero elements. For any non-empty set $\mathcal{S}$ and
a point $x$ both defined in a real vector space we denote by
$\pi_{\mathcal{S}}(x)$ the projection of $x$ onto $\mathcal{S}$, namely,
$\pi_{\mathcal{S}}(x)=\operatorname*{arg\,max}_{y\in\mathcal{S}}\|x-y\|$. With
a slight abuse of notation, when $\pi_{\mathcal{S}}(x)$ is a singleton we will
refer to it as a point in $\mathcal{S}$. For $k$ being the dimension of the
problem, we denote by $\Delta$ the $(k-1)$-simplex (i.e., the unit simplex
with $k$ vertices), and its relative interior by $\Delta^{\mathrm{o}}$;
$\Delta^{\mathrm{o}}=\\{\mathbf{x}+\theta\mathbf{y}\mid\mathbf{x,y}\in\Delta,0<\theta<1\\}$.
Given a probability space we denote by $\mathit{1}(A)$ the r.v. representing
the indicator of an event $A$. For a random variable $X$, $\sigma(X)$ is the
$\sigma$-algebra generated by $X$. The symbols $\>\to_{\rm as}\>$,
$\>\to_{\mathrm{P}}\>$ and $\>\to_{\rm d}\>$ indicate convergence of a
sequence of random variables in almost sure, in probability, and in
distribution, respectively. Almost-sure convergence of a random sequence
$X_{n}$ to a set $\mathcal{S}$, denoted by $X_{n}\>\to_{\rm as}\>\mathcal{S}$,
implies that $\inf_{x\in\mathcal{S}}\|X_{n}-x\|\>\to_{\rm as}\>0$ as
$n\to\infty$. When comparing random variables, $\leq_{\rm st}$ symbolizes
inequality in first-order stochastic dominance.

##  2 Model and preliminaries

We begin by introducing a general formulation of an unobservable queueing
game. By unobservable, we mean that each customer commits to a specific action
(which is possibly chosen randomly) prior to observing any information about
the system’s state. However, we treat the concept of an action here in rather
general abstraction, which allows us to cover many observable models as well
(see discussion in 3.5). Later, in Section 5, we will allow customers to rely
on state-information in their decisions, thereby extending our current
formulation, providing a more refined treatment of observable queueing games.
To ease the exposition, we assume throughout the paper that customers are
homogeneous, however all our results can be generalized to deal with finitely
many heterogeneity types of customers through only incremental changes.

The model. Consider a service system in which the state at time $t$ represents
some (possibly multi-dimensional) buffer content, taking values in a state
space $\mathcal{X}\subset\mathbb{R}^{d}_{+}$. Potential customers arrive at
the system according to a renewal process with inter-arrivals
$\\{A_{n}\\}_{n\geq 1}$ (i.e., $A_{n}$ are iid), where we interpret $A_{n}$ as
the time between the arrivals of the $(n-1)$-st and the $n$-th customers. The
arrival epoch of the $n$-th customer, $n\geq 1$, is therefore given by
$T_{n}=\sum_{i=1}^{n}A_{i}$.

Every arriving customer chooses one out of $k\geq 2$ possible actions, with
the action set denoted by $\mathcal{A}=\\{a_{1},\dots,a_{k}\\}$. Thus, the
$(k-1)$-simplex $\Delta$ represents the set of possible strategies (namely,
distributions over $\mathcal{A}$). For any strategy $\mathbf{p}\in\Delta$ we
denote its $i$-th coordinate by $p_{i}$ which is the probability assigned to
action $a_{i}$.

We denote by $\\{X(t;\mathbf{p})\\}_{t\geq 0}$ the stochastic process
representing the state of the system at time $t$ when the strategy employed by
customers is $\mathbf{p}$. By convention, $\\{X(t;\mathbf{p})\\}_{t\geq 0}$ is
assumed to be right-continuous with left-hand limits. The evolution of
$X(t;\mathbf{p})$ depends on the strategy $\mathbf{p}$ played by the
customers, however, in the general setup we suppress the specifics of how new
arrivals and their corresponding actions change the state of the buffer
content because that depends on the specific application.

######  Example 1.

For the sake of generating intuition, one can think of the canonical
Unobsrvable M/G/1 queue (see [9]), in which customers choose whether or not to
join a single-server queue, thus $k=2$. For concreteness, we assume that
$\lambda<\mu$, where $\lambda$ is the Poisson arrival rate and $1/\mu$ is the
mean service time. In this example, a strategy is characterized by
$\mathbf{p}=(p_{1},p_{2})=(p,1-p)\in\Delta$ for some $p\in[0,1]$, prescribing
the joining and balking probabilities, $p$ and $1-p$, respectively. The
content process $X(t;\mathbf{p})$ corresponds to the workload (or virtual
waiting time), hence when a customer arrives at the system, with probability
$p$ they join the queue and add a random job size (with mean $1/\mu$) to the
workload. $\diamond$

Let $X_{n}(\mathbf{p})=X(T_{n}-;\mathbf{p})$ denote the system state just
before the $n$-th arrival given that customers adopt the strategy
$\mathbf{p}$, and further assume that $X_{1}(\mathbf{p})=0^{d}$. In other
words the system starts empty, so that the first customer arrives (at time
$T_{1}$) to an empty system. The strategy $\mathbf{p}\in\Delta$ defines a
probability measure $\mathbb{P}_{\mathbf{p}}$ for
$\\{X_{n}(\mathbf{p})\\}_{n\geq 1}$. Expectation with respect to the measure
$\mathbb{P}_{\mathbf{p}}$ is denoted by $\mathrm{E}_{\mathbf{p}}$. Define the
_cycle-length_ r.v. $L(\mathbf{p})=\inf\\{n\geq 1\mid
X_{n+1}(\mathbf{p})=0^{d}\\}$. With the assumption $X_{1}(\mathbf{p})=0^{d}$,
$L(\mathbf{p})$ is a r.v. describing the number of arrivals during a typical
regenerative cycle, provided customers play according to $\mathbf{p}$. Note
that separate cycle lengths are iid, that is, for any $n_{1},n_{2},\dots$
satisfying $X_{n_{i}}(\mathbf{p})=0^{d}$, $i=1,2,\dots$, the sequence
$\inf\\{n\geq 1\mid X_{n_{i}+n+1}(\mathbf{p})=0^{d}\\}$ for $i=1,2,\dots$
consists of iid random variables. Let
$\ell(\mathbf{p})=\mathrm{E}_{\mathbf{p}}L(\mathbf{p})$ denote the mean cycle
length when the strategy is given by $\mathbf{p}$, and let
$\ell^{r}(\mathbf{p})=\mathrm{E}_{\mathbf{p}}L^{r}(\mathbf{p})$ denote its
$r$’th moment. By embedding our system at arrival epochs we impose that
$L(\mathbf{p})\geq 1$ with probability 1, therefore $\ell(\mathbf{p})\geq 1$
for every $\mathbf{p}\in\Delta$.

Considering the strategy $\mathbf{p}$ as given, it is known that if the first
moment of the cycle length is finite, i.e., if $\ell(\mathbf{p})<\infty$, then
$\\{X_{n}(\mathbf{p})\\}_{n\geq 1}$ is positive Harris recurrent; for further
details see [6, Ch. VI]. This implies $X_{n}(\mathbf{p})\>\to_{\rm
d}\>X(\mathbf{p})$ as $n\to\infty$, where $X(\mathbf{p})$ is a random variable
corresponding to the stationary distribution at arrival times. Note that
$X(\mathbf{p})$ may be different from the stationary, _time-averaged_
distribution of $X(t;\mathbf{p})$. The conditions for convergence of the
algorithm presented below demand that the cycle length has a finite second
moment, implying that its first moment, $\ell(\mathbf{p})$, is also finite.

The utility of a customer depends on the system state at the time of their
arrival, their chosen action $a_{i}$ and possibly some random outcome.
Formally, given a _realization_ of the state $x\in\mathcal{X}$ and a
realization of the random outcome $y$, we let $v_{i}(x,y)$, $i=1,\dots,k$
denote the value associated with action $a_{i}$. We define

$\mathbf{v}(x,y)=\big{(}v_{1}(x,y),\dots,v_{k}(x,y)\big{)}.$ (1)

The function $\mathbf{v}$ is considered a model’s primitive. Note that the
random outcome may depend on the action chosen, in which case $y$ can be
modeled as a vector $y\in\mathbb{R}^{k}$ with each coordinate corresponding to
an action $a_{i}\in\mathcal{A}$. The dimension of $y$ is in fact irrelevant to
our analysis and for the sake of simplicity $y$ can be thought of as single
valued. The function $\mathbf{v}$, as well as the distribution of the random
outcome, are allowed to depend on the strategy $\mathbf{p}$, but for brevity
we suppress this dependence in the notation.

Our general interest is in treating the state and the random outcome as
(possibly dependent) random variables, say $X$ and $Y$. The value vector
$\mathbf{v}(X,Y)$ then constitutes a $k$-dimensional random variable. We
introduce the vector of expected stationary values as

$\mathbf{u}(\mathbf{p})=\mathrm{E}_{\mathbf{p}}\bigg{[}\mathbf{v}\big{(}X(\mathbf{p}),Y\big{)}\bigg{]},$
(2)

where we recall that $X(\mathbf{p})$ is the random stationary state given
$\mathbf{p}$. We assume throughout that a stationary distribution exists for
every $\mathbf{p}\in\Delta$, and furthermore that the function
$\mathbf{v}(X,Y)$ is integrable with respect to the measure corresponding to
that stationary distribution. This ensures that $\mathbf{u(p)}$ is well
defined given $\mathbf{p}$, thus, it describes a function
$\mathbf{u}:\Delta\to\mathbb{R}^{k}$; For each $i$, $u_{i}(\mathbf{p})$ takes
the interpretation of the mean utility of some “controlled” customer who
arrives at a stationary system and is told to play $a_{i}$. To clarify, we
highlight that $\mathbf{u(p)}$ does _not_ represent the (single-valued)
average utility obtained when customers employ the strategy $\mathbf{p}$ – the
latter can be expressed as $\mathbf{u(p)}^{\prime}\mathbf{p}$.

######  Example 1.

(Continued) Consider once again the example of the unobservable M/G/1, with
arrival and service rates $\lambda$ and $\mu$, respectively, satisfying
$\lambda<\mu$. As explained, a strategy is given by $\mathbf{p}=(p,1-p)$ where
$p\in[0,1]$ depicts the joining probability, thus, we interpret
$X(\mathbf{p})$ as the stationary virtual workload at arrival epochs
associated with customers joining at rate $p\lambda$, and $Y$ is the service
time of a tagged arrival. The possible actions are joining ($a_{1}$) or
balking ($a_{2}$), hence $v_{1}(x,y)=R-C\cdot(x+y)$ for $R$ and $C$ being the
reward from service and cost for unit of delay, and $v_{2}(x,y)=0$. In
addition, letting $w(\cdot)$ denote the mean virtual workload in the system as
a function of the arrival rate, we have
$\mathrm{E}_{\mathbf{p}}[X(\mathbf{p})]=w(p\lambda)$ and therefore
$\mathbf{u}(\mathbf{p})=(R-C\cdot(w(p\lambda)+\mu^{-1}),0)$. $\diamond$

As in Example 1, it is oftentimes the case that $Y$ and $X(\mathbf{p})$ are
independent. However, in general the distribution of $Y$ may depend on the
observed state, $x$, as well as on the strategy of others, $\mathbf{p}$. One
classic example is a processor-sharing system in which the waiting time
depends on the strategy of future arrivals (see [3]).

Equilibrium strategy. Next, we define an equilibrium, the desired solution
concept at the center of attention in this work.

######  Definition 1.

The best-response set for a strategy $\mathbf{p}\in\Delta$ is the set

$\mathcal{BR}(\mathbf{p})=\operatorname*{arg\,max}_{\mathbf{q}\in\Delta}\mathbf{u(p)}^{\prime}\mathbf{q}.$

For a strategy $\mathbf{p}\in\Delta$, each of the elements of
$\mathcal{BR}(\mathbf{p})$ is commonly termed a _best-response strategy_ (or
_best response_ in short) for $\mathbf{p}$. Thus, $\mathbf{p}$ is a symmetric
equilibrium if it is a best response to itself, or synonymously, if it is a
fixed point of $\mathcal{BR}$, when the latter is viewed as a set-valued map,
$\mathcal{BR}:\Delta\to 2^{\Delta}$. Hence,

######  Definition 2.

A Symmetric Nash Equilibrium strategy is a strategy $\mathbf{p}\in\Delta$ such
that

$\mathbf{p}\in\mathcal{BR}(\mathbf{p})=\operatorname*{arg\,max}_{\mathbf{q}\in\Delta}\mathbf{u(p)}^{\prime}\mathbf{q}.$

Nonetheless, describing an equilibrium as a fixed point of the best-response
function $\mathcal{BR}$ does not yield any simple method to compute it in
general. This is because $\mathcal{BR}$ is a set-valued mapping in a
continuous space, and is naturally hard to work with. Moreover, even if
$\mathcal{BR}$ is a singleton almost everywhere in the domain $\Delta$ and
continuous at every such point, it is often the case that the fixed point
$\mathbf{p}^{e}\in\mathcal{BR}(\mathbf{p}^{e})$ is a point of jump
discontinuity, in the sense that the limits approaching $\mathbf{p}^{e}$ from
different sides do not agree.

######  Example 1.

(Continued) In the unobservable M/G/1, recall that a strategy is given by
$\mathbf{p}=(p,1-p)$ with $p\in[0,1]$ and that
$\mathrm{E}_{\mathbf{p}}[X(\mathbf{p})]=w(p\lambda)$ which is the mean virtual
workload when the arrival rate is $p\lambda$. Suppose that $R,C$ and $\mu$ are
such that $0<R/C-1/\mu<w(\lambda)$. Then there exists a unique equilibrium
joining probability $p^{e}\in(0,1)$, which is characterized by the unique
solution to the equation $R-C\cdot(w(p^{e}\lambda)+1/\mu)=0$, and the
corresponding equilibrium strategy is $\mathbf{p}^{e}=(p^{e},1-p^{e})$. The
best-response function (plotted in Figure 2) for a strategy
$\mathbf{p}=(p,1-p)$ takes the form:

$\mathcal{BR}(\mathbf{p})=\begin{cases}\mathbf{e}_{1}\quad&\mbox{if
}p<p^{e};\\\ \Delta&\mbox{if }p=p^{e};\\\ \mathbf{e}_{2}&\mbox{if
}p>p^{e}.\end{cases}$

Denoting $\tilde{\mathbf{e}}=\mathbf{e}_{1}-\mathbf{e}_{2}$, it can be seen
that for any $\epsilon>0$,

$\|\mathcal{BR}(\mathbf{p}^{e}+\epsilon\tilde{\mathbf{e}})-\mathcal{BR}(\mathbf{p}^{e}-\epsilon\tilde{\mathbf{e}})\|=\|\mathbf{e}_{1}-\mathbf{e}_{2}\|=\|\tilde{\mathbf{e}}\|=\sqrt{2}.$

$\diamond$

Example 1 shows that the set-valued mapping $\mathcal{BR}$ in general is not
lower-hemicontiuous. Thus, even in the fundamental unobservable M/G/1 setting,
it is not obvious how to approach the equilibrium strategy based on fixed-
point iterations. To overcome this difficulty, we define next a modified
version of the best-response function, which is a vector-valued function,
$\mathbf{f}:\Delta\to\Delta$, expressed as:

$\mathbf{f}(\mathbf{p})=\pi_{\Delta}\big{(}\mathbf{p}+\mathbf{u(p)}\big{)},$
(3)

recalling that $\pi_{\Delta}(x)$ denotes the projection of $x$ onto the
simplex $\Delta$. This vector-valued function can replace the raw definition
of $\mathcal{BR}$ for our purpose, and in addition, under appropriate
conditions on the primitives, also possesses desirable attributes like
continuity and smoothness.

The intuition behind the definition of $\mathbf{f}$ is that
$\mathbf{f}(\mathbf{p})$ maps $\mathbf{p}$ to a strategy obtained by deviating
from $\mathbf{p}$ in the direction of one of its (actual) best responses. This
can be noticed observing that $\mathbf{p}+\mathbf{u(p)}$ is the unique
solution of the following (convex) optimization problem:

$\max_{\mathbf{q}\in\mathbb{R}^{k}}\left\\{\mathbf{u(p)}^{\prime}\mathbf{q}-\frac{1}{2}\|\mathbf{p-q}\|^{2}\right\\}.$

Whenever the best response for $\mathbf{p}$ is not unique, $\mathbf{f(p)}$
“pushes” $\mathbf{p}$ in the direction of the best response that is the
closest to $\mathbf{p}$ in the standard euclidean sense (note that the
projection of $\mathbf{p}$ onto $\mathcal{BR}(\mathbf{p})$ is unique because
the latter is a convex polyhedron in $\mathbb{R}^{k}$). It is intuitive
therefore that a strategy $\mathbf{p}$ is a symmetric equilibrium if and only
if it is a fixed point of $\mathbf{f}$. This statement is rigorously
formulated in the next lemma, whose proof appears in Appendix A.1:

######  Lemma 2.

A strategy $\mathbf{p}\in\Delta$ is a Symmetric Nash Equilibrium, i.e.,
$\mathbf{p}\in\mathcal{BR}(\mathbf{p})$, if and only if it satisfies
$\mathbf{p}=\mathbf{f}(\mathbf{p})$.

$p$$45^{\circ}$0101$p^{e}$ first component
---
of $\mathcal{BR}(\mathbf{p})$
first component
---
of $\mathbf{f}(\mathbf{p})$
Figure 2: The first component of the best-response function
$\mathcal{BR}(\mathbf{p})$ (red), and the first component of the best-response
surrogate $\mathbf{f(p)}$ (blue), for the unobservable M/G/1 model. At $p^{e}$
customers are indifferent between joining and balking, thus, every joining
probability $p\in[0,1]$ constitutes a best-response strategy,
$(p,1-p)\in\mathcal{BR}(\mathbf{p}^{e})$. The red curve’s left and right
limits at $p^{e}$ do not coincide. The blue curve is a single-valued
continuous function over $[0,1]$ and is smooth at $p^{e}$.

In Figure 2, both the best-response function, $\mathcal{BR}$, and its
surrogate, $\mathbf{f}$, are illustrated for the unobservable M/G/1 model of
Example 1. It can be seen that both functions admit a unique fixed point at
the equilibrium strategy. However, $\mathcal{BR}$ is discontinuous (in the
sense of lower-hemicontinuity) at that point, while $\mathbf{f(p)}$ is indeed
continuous.

######  Remark 1.

The result of Lemma 2 would still prevail if we extended our definition of
$\mathbf{f(p)}$ to
$\pi_{\Delta}\big{(}\mathbf{p}+h(\mathbf{p})\mathbf{u(p)}\big{)}$, for any
positive real function $h:\Delta\to\mathbb{R}$. In using our suggested
equilibrium-approximation scheme, it is useful to choose a function $h$ that
is bounded away from 0.

Iterative equilibrium-approximation scheme. Equipped with Lemma 2 we focus our
interest on characterizing a solution to $\mathbf{p}=\mathbf{f}(\mathbf{p})$.
To ensure that such a fixed point exists, we assume that $\mathbf{u}$ is
continuous, in which case we have, as an implication of the maximum theorem,
that $\mathbf{f}$ is continuous as well. We can then state the following
lemma.

######  Lemma 3.

Assume that $\mathbf{u(p)}$ is continuous for all $\mathbf{p}\in\Delta$. Then
a symmetric equilibrium strategy $\mathbf{p}^{e}\in\Delta$ exists, and this
strategy satisfies $\mathbf{p}^{e}=\mathbf{f}(\mathbf{p}^{e})$.

Lemma 3 follows directly from Lemma 2 and the continuity of $\mathbf{f}$ by
applying Brouwer’s fixed-point theorem. As a consequence of Lemma 3, when
$\mathbf{u(p)}$ is accessible and can be evaluated directly, then a simple
fixed-point iteration algorithm can be used to approximate an equilibrium; To
this end, we think of $\mathbf{u(p)}$ as the vector that determines the
direction of the step at each iteration, giving rise to the following
(deterministic) adaptive-iteration scheme:

${\mathbf{p}}^{(n+1)}=\pi_{\Delta}\left({\mathbf{p}}^{(n)}+\gamma_{n}\mathbf{u}({\mathbf{p}}^{(n)})\right),$
(4)

where $\\{\gamma_{n}\\}_{n\geq 1}$ is a positive sequence called the step-size
sequence, and satisfies $\sum_{n=1}^{\infty}\gamma_{n}=\infty$. By choosing
$\\{\gamma_{n}\\}_{n\geq 1}$ such that $\gamma_{n}\to 0$, assuming
$\mathbf{u}(\cdot)$ is continuous, any limit point of the algorithm in (4)
corresponds to an equilibrium strategy. Yet, the latter statement alone does
not make any conclusion about whether or not convergence takes place, and in
some cases the iterate may indeed end up oscillating between non-equilibrium
strategies. Further regulating assumptions on $\mathbf{u}(\cdot)$ are needed
(for example, Assumption A3 below) to guarantee convergence to an equilibrium.
General conditions under which (4) converges is a rather technical issue that
exceeds the scope of our discussion here; some relevant information is
provided in Appendix A.2.

######  Remark 2.

Consider a point $\mathbf{p}^{e}$ in the (relative) interior of $\Delta$,
$\mathbf{p}^{e}\in\Delta^{\mathrm{o}}$. Define

$\mathbf{z}=\pi_{\Delta}\left(\mathbf{p}^{e}+\mathbf{u}(\mathbf{p}^{e})\right)-(\mathbf{p}^{e}+\mathbf{u}(\mathbf{p}^{e})),$

which is the shortest vector needed to bring
$\mathbf{p}^{e}+\mathbf{u}(\mathbf{p}^{e})$ back into $\Delta$. A necessary
condition for the algorithm in (4) to converge to $\mathbf{p}^{e}$ is
$\mathbf{z}=-\mathbf{u}(\mathbf{p}^{e})$, i.e., the shortest trajectory from
$\mathbf{p}^{e}+\mathbf{u}(\mathbf{p}^{e})$ back to the simplex has to be
along the opposite direction of $\mathbf{u}(\mathbf{p}^{e})$. When
$\mathbf{p}^{e}$ is an interior point, the latter condition means that
$\mathbf{u(p}^{e})$ has to be orthogonal to the hyperplane spanned by
$\Delta$. In other words, all entries of $\mathbf{u(p}^{e})$ have to be equal.
This observation falls in line with the well-known _indifference principle_ :
If $\mathbf{p}$ is a strategy such that when customers play $\mathbf{p}$, an
individual customer is indifferent to choosing any action in $\mathcal{A}$,
then $\mathbf{p}$ induces equilibrium. The latter argument and its intuitive
interpretation can be extended to any strategy
$\mathbf{p}=(p_{1},\dots,p_{k})$ that lies in a facet of $\Delta$, with its
corresponding subset of supported actions, $\\{a_{i}\in\mathcal{A}\mid i\text{
\rm{ s.t.} }p_{i}>0\\}$.

Clearly, the deterministic algorithm just described is only useful when the
utility function $\mathbf{u(p)}$ can be computed. In many settings, however,
the entries of $\mathbf{u(p)}$ may not be given in closed form, and even
numerical approximations can be difficult to obtain. Therefore we devise a
method of approaching the desired equilibrium using an SA algorithm.

##  3 Simulation and stochastic approximation

The stochastic approximation procedure is based on the Robbins-Monro
algorithm, which aims to mimic the deterministic fixed-point iteration of (4),
replacing the progression direction $\mathbf{u(p)}$ with a noisy estimator
based on a realization of the queueing process. However, estimating
$\mathbf{u(p)}$ directly, e.g by taking a sample average over a (finitely-
long) simulation path, often results in a biased estimator. Using such biased
estimators in the SA scheme, one risks the algorithm will converge to an
undesirable limit. In light of Remark 1, we thus replace $\mathbf{u(p)}$ in
(4) with $\mathbf{g(p)}=\ell(\mathbf{p})\mathbf{u(p)}$, which we can easily
estimate, without bias, based on realizations of regeneration cycles. A
typical iteration of the SA algorithm then runs as follows. At iteration $n$,
assume that the strategy is $\mathbf{p}^{(n)}$, and that we are provided with
an estimator $\mathbf{G}^{(n)}$ for $\mathbf{g}(\mathbf{p}^{(n)})$. Recall
that $\gamma_{n}$ is the $n$-th step size. The update at the next iteration of
the SA scheme is done by

$\mathbf{p}^{(n+1)}=\pi_{\Delta}\left(\mathbf{p}^{(n)}+\gamma_{n}\mathbf{G}^{(n)}\right),\
n\geq 1,$ (5)

where $\mathbf{p}^{(1)}\in\Delta$ is an arbitrary initial strategy. If
$\mathbf{G}^{(n)}$ is unbiased for $\mathbf{g}(\mathbf{p}^{(n)})$, then under
mild regularity conditions the algorithm will converge to a fixed point of
$\mathbf{f}$ (see Theorem 5), i.e., to a Nash equilibrium. The remainder of
this section is dedicated to the construction of an unbiased estimator and to
presenting the necessary assumptions and steps for our main convergence
result.

###  3.1 Unbiased estimator

Assume that the strategy $\mathbf{p}$ is given. To avoid using cumbersome
notation we suppress in some places the dependency on $\mathbf{p}$, however
all the distributions of the random variables to be defined next should be
regarded as functions of $\mathbf{p}$. To construct the estimator, we simulate
a single regenerative cycle starting with the arrival of a customer at an
empty system; $X_{1}(\mathbf{p})=0^{d}$. As before, we denote by
$L=\inf\\{n\geq 1\mid X_{n+1}(\mathbf{p})=0\\}$ the number of arrivals during
that cycle. For each arrival $j=1,\dots,L$, let $X_{j}$ be the state and
$Y_{j}$ be the corresponding random outcome. For the $j$-th customer, the
vector $\mathbf{v}(X_{j},Y_{j})$ as defined in (1) reflects the potential
utility that each action will yield to that customer. We clarify that
calculating $\mathbf{v}(X_{j},Y_{j})$ for each customer $j$ requires to assess
their realized value for _every_ action in $\mathcal{A}$ _as if_ they chose
it, yet in the simulation process, the action taken by the $j$-th customer is
fully governed by $\mathbf{p}$.

Denote the expectation of the utility conditional on observing state $X$ by

$\mathbf{\overline{v}}(X)=\mathrm{E}_{Y}\left[\mathbf{v}\big{(}X,Y\big{)}\>\big{|}\>X\right].$
(6)

Throughout the paper we assume, for all $\mathbf{p}\in\Delta$, that
$\mathbf{\overline{v}}$ is integrable with respect to the distribution of
$X(\mathbf{p})$. Then the vector of expected stationary utilities,
$\mathbf{u}(\mathbf{p})$, can be written as
$\mathbf{u}(\mathbf{p})=\mathrm{E}_{\mathbf{p}}\big{[}\mathbf{\overline{v}}(X(\mathbf{p}))\big{]}$.
Define

$\mathbf{G}=\sum_{j=1}^{L}\mathbf{\overline{v}}(X_{j}).$ (7)

The following lemma, which is a direct result of a classical formula of
regenerative processes (see, for example, [19, Prop. A.3]), states that
$\mathbf{G}$ is an unbiased estimator for $\mathbf{g(p)}$.

######  Lemma 4.

Suppose that $\ell^{2}(\mathbf{p})<\infty$ for all $\mathbf{p}\in\Delta$, then

$\mathrm{E}_{\mathbf{p}}\mathbf{G}=\ell(\mathbf{p})\mathbf{u(p)}=\mathbf{g(p)}.$
(8)

Lemma 4 explains the rational behind the choice of an estimate for
$\mathbf{g(p)}$ over an estimate for $\mathbf{u(p)}$. As discussed, one may
attempt to estimate $\mathbf{u(p)}$ directly, for example by considering the
naive estimator

$\mathbf{Z}=\frac{1}{L}\sum_{j=1}^{L}\mathbf{\overline{v}}(X_{j}),$

however, such an estimator is generally biased. To see why, we provide an
example.

######  Example 1.

(Continued) Recall the unobservable M/G/1 model with service rate $\mu$,
reward $R$, and delay cost $C$, where we used $\mathbf{p}=(p,1-p)$ to denote a
strategy and $X_{j}$ to describe the virtual workload at the $j$-th arrival,
so that $\mathrm{E}_{\mathbf{p}}[X(\mathbf{p})]=w(p\lambda)$ is the mean
virtual workload when the arrival rate is $p\lambda$. In particular, the
estimator $\mathbf{Z}$ when constructed from a single regeneration cycle
(i.e., one busy period) is a 2-dimensional vector,
$\mathbf{Z}=(R-C\cdot(\bar{X}+1/\mu),0)$, where
$\bar{X}=(1/L)\sum_{j=1}^{L}X_{j}$ is the average workload at arrival instants
over all customers of that cycle. It is known that $\bar{X}$ is a biased
estimator of $w(p\lambda)$, i.e., that $\mathrm{E}_{\mathbf{p}}[\bar{X}]\neq
w(p\lambda)$ (see, for example, [19]), and so, $\mathbf{Z}$ is biased for
$\mathbf{u(p)}$. The fact that $\bar{X}$ is biased can be explained as a
ramification of the well-known _length bias_ : The average virtual workload
over long busy periods will naturally tend to be larger than that over short
busy periods. Yet, during a long busy period more customers will typically
arrive at the system than in a short one. By averaging the sampled workloads
one erroneously normalizes the length of the busy period, and as a result,
gives disproportional importance to short busy periods in the estimation,
resulting in an underestimate of the mean workload. By contrast, the estimator
$\mathbf{G}$ presented above, for the unobservable M/G/1 model takes the form

$\mathbf{G}=\sum_{j=1}^{L}\mathbf{\overline{v}}(X_{j})=\sum_{j=1}^{L}\begin{pmatrix}R-C\cdot(X_{j}+1/\mu)\\\
0\end{pmatrix}.$

This estimator is indeed unbiased for $\mathbf{g(p)}$, and in particular,
$\mathrm{E}_{\mathbf{p}}\left(\sum_{j=1}^{L}X_{j}\right)=\ell(\mathbf{p})w(\lambda
p)$. $\diamond$

It is important to point out the difference between using the realized value
vector for customer $j$, $\mathbf{v}(X_{j},Y_{j})$, and using the conditional
expectation $\mathbf{\overline{v}}(X_{j})$, as prescribed in (7). The latter
has two main advantages. Firstly, using the expected value instead of the
realization reduces the variance of the estimator. Secondly,
$\mathbf{\overline{v}}$ is merely a function of the state, and hence is more
convenient to work with. Yet, in some settings it may not be trivial to
express this function explicitly.

In principle, one can directly use $\mathbf{v}(X_{j},Y_{j})$ in Equation (7)
instead of $\mathbf{\overline{v}}(X_{j})$. However, since
$\mathbf{v}(X_{j},Y_{j})$ is a function of both the observed state $X_{j}$ and
the random outcome $Y_{j}$, in order to apply the result of Lemma 4, one must
make assumptions about the dependence between $X_{j}$ and $Y_{j}$. For
example, the following assumption is sufficient: At the $j$-th arrival with
the observed state $x_{j}$, the random outcome conditioned on the state
observed, $Y_{j}\mid X_{j}=x_{j}$, should be independent of $X_{j^{\prime}}$
and $Y_{j^{\prime}}$ for all $j^{\prime}\neq j$. This assumption may not hold
in systems with elaborate service policies such as priority queues or
processor sharing: For instance, in processor sharing, the waiting of the
$j$-th arriving customer indeed depends on the service demand (as well as the
action) of later-arriving customers. Furthermore, the value
$\mathbf{v}(X_{j},Y_{j})$ is not a byproduct of a single busy period
simulation because it must be computed for every possible action of a
customer, including those not taken in the realization of the simulation. A
possible solution is, for each $j=1\ldots,L$, to obtain an independent sample
of $\mathbf{v}(x_{j},Y_{j})$ using subordinate simulations, independently of
the main simulation path. To generate a sample of any component
$v_{i}(x_{j},Y_{j})$ of the vector $\mathbf{v}(x_{j},Y_{j})$, one can run an
independent subordinate simulation that starts with the arrival of a customer
at state $x_{j}$ who takes action $a_{i}$, and terminate when the utility of
that customer is realized. When possible, the same subordinate simulation can
be used to evaluate as many components of the vector $\mathbf{v}(x_{j},Y_{j})$
as needed, as long as this is done independently of the main simulation path.

###  3.2 Convergence and main result

In the approximation process, given a strategy $\mathbf{p}$, we make use of
the estimator $\mathbf{G}$ at each iteration of the algorithm as a “proxy” for
the true value $\mathbf{g(p)}$. Following the framework of [49, Ch. 5] we
present conditions for this algorithm to converge almost surely to a set of
equilibrium points. In particular, we assume that the variance of the unbiased
estimator $\mathbf{G}$ is finite for any strategy $\mathbf{p}$. This is
combined together with the assumptions required for existence of an
equilibrium, as discussed in Section 2, to obtain the convergence result
described in Theorem 5. The proof of that theorem appears in Appendix A.2.

To state our assumptions, we first introduce the hyperplane
$\mathcal{H}=\\{\mathbf{x}\in\mathbb{R}^{k}\mid\mathbf{e}^{\prime}\mathbf{x}=0\\}$
which is a $(k-1)$-dimensional subspace of $\mathbb{R}^{k}$, parallel to
$\Delta$. We further note that the orthogonal projection onto $\mathcal{H}$,
$\pi_{\mathcal{H}}:\mathbb{R}^{k}\to\mathcal{H}$, is a linear transformation,
and can be represented by a symmetric matrix $\mathbf{H}\in\mathbb{R}^{k\times
k}$, namely, $\pi_{\mathcal{H}}(\mathbf{x})=\mathbf{Hx}$.

######  Assumption A1.

The second moment of the cycle length $L(\mathbf{p})$ is bounded for every
strategy; $\ell^{2}(\mathbf{p})<\infty$ for all $\mathbf{p}\in\Delta$.

######  Assumption A2.

The second moment of the estimator $\mathbf{G}$ is uniformly bounded on
$\Delta$;
$\sup_{\mathbf{p}\in\Delta}\mathrm{E}_{\mathbf{p}}\|\mathbf{G}\|^{2}<\infty$.

######  Assumption A3.

The projected utility function $\mathbf{Hu}(\cdot)$ is continuous, and there
exists a differentiable function $u^{*}:\Delta\to\mathbb{R}$ such that
$-\mathbf{Hu}(\mathbf{p})$ is the gradient of $u^{*}$ at $\mathbf{p}$, i.e.,
$-\mathbf{Hu}(\mathbf{p})=\nabla u^{*}(\mathbf{p})$.

######  Assumption A4.

The step-size sequence $\\{\gamma_{n}\\}_{n\geq 1}$ satisfies

$\sum_{n=1}^{\infty}\gamma_{n}=\infty,\quad\sum_{n=1}^{\infty}\gamma_{n}^{2}<\infty,$
(9)

######  Theorem 5.

Suppose Assumptions A1–A4 are satisfied. As $n\to\infty$,
$\mathbf{p}^{(n)}\>\to_{\rm as}\>\mathcal{S}^{e}$, where
$\mathcal{S}^{e}\subseteq\Delta$ is a set of equilibrium strategies, namely,
$\mathbf{f}(\mathbf{p})=\mathbf{p}$ for every $\mathbf{p}\in\mathcal{S}^{e}$.

######  Corollary 6.

Under the conditions of Theorem 5, if the set of equilibrium strategies is
finite, then as $n\to\infty$, $\mathbf{p}^{(n)}\>\to_{\rm
as}\>\mathbf{p}^{e}$, where $\mathbf{p}^{e}$ is an equilibrium.

###  3.3 Discussion of assumptions

Assumption A1 is required for positive recurrence of the underlying
regenerative process, which in turn ensures the existence of a stationary
distribution for any strategy. The step-size condition of Assumption A4 is
standard and enables the application of the law of large numbers for
martingale differences, which underlies the convergence result in Theorem 5.
Assumptions A2–A3 impose regularity conditions on the utility function which
are specific to the framework of queueing games, and therefore deserve more
attention, as discussed below.

A straightforward generalization of Lemma 3 shows that the continuity of
$\mathbf{Hu}(\cdot)$ implies that an equilibrium exits, in accordance with the
first part of Assumption A3. Stated in different words, the second part of the
assumption, which relates to the existence of some potential function $u^{*}$,
means that $\mathbf{Hu}(\cdot)$ should be a conservative vector field on
$\Delta$. Still, the function $u^{*}$ need not be specified. Its existence is
needed to ensure convergence, and prevents the iterate from circulating a
perpetual loop of off-equilibrium strategies. Importantly, without this
assumption there is no guarantee that the deterministic algorithm converges to
an equilibrium, that is, even when $\mathbf{u(p)}$ is known. For a continuous
utility $\mathbf{u}(\cdot)$, Assumption A3 holds, for example, when $k=2$, or
alternatively when $k>2$ and $\mathbf{u}$ itself is a conservative vector
field on $\Delta$. Further elaboration on Assumption A3 is provided in
Appendix A.2.

Assumption A2 imposes that the variance of the estimator $\mathbf{G}$ is
finite. Finite variance of the increments is quite a common assumption in the
framework of stochastic approximation. However, verifying it for specific
queueing games with only limited knowledge about the stationary distribution
can be challenging. In this paper we provide specific examples in which we
rely on queueing-analytic tools to verify this assumption. Example 1 below
provides an analogous condition to Assumption A2 for the unobservable M/G/1
queue, which is formulated in terms of moments of the service-time
distribution.

######  Example 1.

(Continued) Consider once more the unobservable M/G/1 model described in
Sections 2 and 3: Denote the generic service time random variable by $Y$ with
$\mathrm{E}Y=1/\mu$. As before, we assume that $\lambda<\mu$. We further
assume $\mathrm{E}Y^{4}<\infty$, which clearly implies
$\mathrm{E}Y^{2}<\infty$. In this case Assumptions A1–A2 can be verified by
applying known properties, which is done in Section 4.1. In the particular
unobservable M/G/1 model, $\mathrm{E}Y^{4}<\infty$ is also a necessary
condition for Assumption A2; Recall that $G=\sum_{i=1}^{L}(R-C(X_{i}+1/\mu))$,
where $X_{i}$ is the virtual workload observed by customer $i$ and $Y_{i}$ is
their service time. Hence, $G$ is a linear function of the sum of waiting
times during a busy period. An explicit formula of the variance of the
cumulative workload during a busy period was given in [22] as a function of
$\lambda$ and $\mathrm{E}Y^{j}$ for $J=1,2,3,4$. Therefore, for the strategy
$\mathbf{p}=\mathbf{e}_{1}$ such that all customers join, $\lambda<\mu$ and
$\mathrm{E}Y^{4}<\infty$ imply

$\mathrm{E}_{\mathbf{e}_{1}}\|\mathbf{G}\|^{2}=\mathrm{E}_{\mathbf{e}_{1}}\Bigg{\|}\sum_{j=1}^{L}\mathbf{\overline{v}}(X_{j}(\mathbf{e}_{1}))\Bigg{\|}^{2}<\infty.$

Furthermore, by standard coupling arguments it can be verified that for every
strategy $\mathbf{p}\in\Delta$,

$\mathrm{E}_{\mathbf{p}}\|\mathbf{G}\|^{2}\leq\mathrm{E}_{\mathbf{e}_{1}}\|\mathbf{G}\|^{2}.$

We conclude that $\lambda<\mu$ together with $\mathrm{E}Y^{4}<\infty$ form a
necessary and sufficient condition for Assumption A2. Theorem 5 then implies
that under these conditions the SA algorithm converges almost surely to the
unique equilibrium. $\diamond$

Assumptions A1–A2 demand that the conditions are satisfied for all strategies
in the simplex, but this can in some cases be relaxed with a small
modification of the algorithm. If there exists a known closed convex set of
strategies $\mathcal{S}\subseteq\Delta$ that contains an equilibrium, and
conditions A1–A2 are satisfied for all $\mathbf{p}\in\mathcal{S}$, then the
projection $\pi_{\mathcal{S}}$ can be used in (5) and all our results carry
over. In practice, such a set can often be found because stability conditions
are known for many models even without explicit derivation of the stationary
distribution. However, in absence of a-priori knowledge of the stability
region, the algorithm must be modified to incorporate ‘stability-enforcing’
dynamics. Appendix C presents a heuristic modification of the algorithm that
terminates busy cycles if they get too long. Of course, this introduces a bias
in the estimation step. To eliminate this bias, the modified algorithm
increases the termination threshold gradually, with the goal of making the
bias asymptotically negligible once the algorithm is absorbed in the stability
region.

That said, we acknowledge that the task of detecting stability regions using
simulation is challenging in general. In fact, the problem of determining
whether a queueing system is stable under a general schedueling policy was
shown in [29] to be undecidable even in a simple system configuration. A
simulated-annealing method to identify stability regions in some practical
setups is suggested in [53]. Some ad-hoc solutions have been discussed in
recent papers that involve learning and control algorithms for queueing
systems (see [21] and [18]).

###  3.4 Rate of convergence and approximate equilibrium

An important question from a practical point of view is how many iterations
are needed to obtain a ‘good’ approximation. This question has been the
subject of interest for a large body of works in the literature the stochastic
gradient-descent algorithms, where it is often implied by appropriate
conditions on the primitives that the iterate $\mathbf{p}^{(n)}$ converges
(almost surely) to a unique point $\mathbf{p}^{e}$. Then, the above question
is addressed via studying different notions of convergence of the scaled
approximation error, $(\mathbf{p}^{(n)}-\mathbf{p}^{e})/\sqrt{\gamma_{n}}$.
For a specific choice of $\\{\gamma_{n}\\}$, [60] established
$L^{2}$-convergence of the scaled error, by imposing strong convexity and
smoothness on the potential function ${u}^{*}$. Under these assumptions, the
equilibrium exists uniquely (as it is the minimum point of $u^{*}$), and the
vector $-\mathbf{Hu}$ is the gradient for which we have noisy observations. In
what follows we leverage this framework to bound our algorithm’s equilibrium-
approximation error with high probability.

Rather than studying the scaled error per se, we find it more meaningful in
our setup to understand, given the iteration index $n$, how ‘close’ the
strategy $\mathbf{p}^{(n)}$ is to meeting the equilibrium condition. We
therefore adopt the concept of _$\epsilon$ -approximate Nash equilibrium_ (see
[23]). Recall that a strategy $\mathbf{p}$ is a (symmetric)
$\epsilon$-approximate Nash equilibrium, or an $\epsilon$-equilibrium for
short, if no deviation from this strategy can increase a customer’s expected
utility by more than $\epsilon$:

$\mathbf{u}(\mathbf{p})^{\prime}\mathbf{p}\geq\max_{\mathbf{q}\in\Delta}\mathbf{u}(\mathbf{p})^{\prime}\mathbf{q}-\epsilon.$
(10)

Proposition 7 below relates the convergence of the scaled errors to
$\epsilon$-equilibria. The proof is provided in Appendix A.3.

######  Proposition 7.

Suppose that Assumptions A1–A3 hold. Suppose also that $\mathbf{u}$ is
(locally) Lipschitz on $\Delta$, and that $u^{*}$ in A3 is strongly convex on
$\Delta$, i.e., there exists a constant $C>0$ such that

$u^{*}(\mathbf{p})\geq
u^{*}(\mathbf{q})+\mathbf{Hu}(\mathbf{q})^{\prime}(\mathbf{p}-\mathbf{q})+\frac{C}{2}\|\mathbf{p}-\mathbf{q}\|^{2},\
\forall\mathbf{p,q}\in\Delta.$ (11)

Let $\gamma_{n}=\eta/n$, where $\eta>1/(2C)$. Then there exists a constant
$M$, such that for all $n\geq 1$ and $\delta\in(0,1)$,

$\mathrm{P}\left(\mathbf{u}(\mathbf{p}^{(n)})^{\prime}\mathbf{p}^{(n)}\geq\max_{\mathbf{q}\in\Delta}\mathbf{u}(\mathbf{p}^{(n)})^{\prime}\mathbf{q}-\sqrt{n^{-\delta}}\right)\geq
1-M\cdot n^{\delta-1};$ (12)

that is, for all $n$ and $\delta\in(0,1)$, the $n$-th iterate
$\mathbf{p}^{(n)}$ is a $\sqrt{n^{-\delta}}$-equilibrium with probability (at
least) $1-M\cdot n^{\delta-1}$.

Notice that the choice of the step size $\gamma_{n}=\eta/n$ in Proposition 7
satisfies A4. Following [60], under the conditions of Proposition 7, this
choice of $\gamma_{n}$ implies that
$\mathrm{E}\|\mathbf{p}^{(n)}-\mathbf{p}^{e}\|^{2}=\mathrm{o}(1/n)$.
Proposition 7 shows this also yields an $\epsilon_{n}$-equilibrium with high
probability ($1-\mathrm{o}(n^{\delta-1})$), where
$\epsilon_{n}=\sqrt{n^{-\delta}}$. For example, with $\delta=1/2$, the result
ensures that an $\epsilon$-equilibrium is reached, with
$1-\mathrm{o}(\epsilon^{2})$ certainty, after $\epsilon^{-4}$ iterations.

There are various other conditions in the literature for convergence of the
scaled error, including almost-sure convergence rates (e.g., [66] for
unconstrained stochastic optimization) and weak convergence to the normal
distribution (e.g., [12, Ch. 8.4]). All of these results require assumptions
on the utility vector $\mathbf{u}$ and the penitential function $u^{*}$ that
ensure the algorithm is “well behaved”, so to speak, in the neighborhood of
equilibrium points. In Proposition 7 this takes the form of (11) and the
corresponding step-size constant $\eta>1/(2C)$. Though in applications these
conditions can be hard to verify, our result does shed light on the obtainable
rates when some limited structure of the stationary expected utility is
available. For a detailed discussion of the convergence rate for constrained
and unconstrained stochastic-approximation algorithms, see [14] and [49].

###  3.5 A direct extension to observable queueing games

In the context of queueing games, it is natural to think of an ‘action’ as
representing a possible terminal outcome of a customer decision process, for
example, ‘join’ or ‘balk’, as in Example 1. However, in our framework, the
term ‘action’ is in fact synonymous with ‘pure strategy’ (this interpretation
slightly deviates from the terminology used in extensive-form games, see [58,
Ch. 3]). Thus, a single action can potentially represent a complex behavior
rule, which, for instance, can prescribe joining at some queue-lengths while
balking at others.

As a brief example, consider an observable GI/G/1/$N$ queue: Each customer
observes the number $n$ of customers in the queue, and assuming it is not full
($0\leq n\leq N-1$), they choose whether to join or balk. Then a pure strategy
is an element in $\\{0,1\\}^{N}$ determining which queue lengths to join, and
we associate it with a single action in the (finite) set $\mathcal{A}$.

Observable queueing games with finitely many pure strategies are covered by
the formulation in section 2 and hence our results can indeed be applied to
such games. However it can be noticed, particularly in the example above, that
the dimension of the strategy space grows exponentially with the number of
signals potentially observed. When the number of signals is large, the
approach discussed in this section may not be practical. In Section 5 we
discuss an extension of the algorithm that traverses a strategy space whose
dimension is proportional (as opposed to exponential) to the number of
signals.

##  4 Applications, implementation and refinements

This section presents two applications of the SA algorithm. The first,
discussed in Section 4.1, is a general game of choosing a queue from a pool of
parallel single-server queues, or balking from the system altogether. This
model, inspired by [9], generalizes the canonical M/G/1 model of Example 1 and
also the motivating example of two parallel queues in Section 1.1. The second
example, in Section 4.2 is a generalization of the model introduced in [38],
considering a system comprised of two parallel servers; one with an infinite
buffer queue and one with no queueing buffer, and customers have the option to
sense the no-buffer server at a cost before joining the infinite buffer queue.
For both applications the sufficient convergence conditions of Theorem 5 are
verified in detail, and some numerical results are reported. The related
proofs appear in Appendix A. Using the second example we focus on reducing the
variance of the estimator and dynamically choosing the step size. The
modifications are illustrated by considering an example model.

###  4.1 Unobservable GI/G/1 queues in parallel

The main purpose of the section is to rigorously verify Assumption A1–A2,
which imply the convergence of the SA algorithm. The proof of Proposition 8 is
provided in Appendix B.2 and relies on a coupling construction, which we
believe can be found useful in verifying the assumptions for further
generalizations and variations of the model here.

For some integer $k>1$, consider a system of $k-1$ parallel GI/G/1 queues
where customers can choose which queue to join, with the possibility of
balking, i.e., not joining any queue at all. Arrivals to the system are
generated according to a renewal process with a continuous inter-arrival
distribution $H$ and mean $1/\lambda$. Each arriving customer chooses, without
observing the system state, whether at all to join a queue, and if so, which
of the $k-1$ queues to join. Denote the generic service time variable in
station $m\in\\{1,\ldots,k-1\\}$ by $Y_{m}\sim F_{m}$, and its mean by
$1/\mu_{m}$. A strategy $\mathbf{p}=(p_{1},\dots,p_{k})\in\Delta$, for each
$m\in\\{1,\ldots,k-1\\}$ expresses the probability $p_{m}$ of joining queue
$m$, with $p_{k}$ being the probability of balking. Arriving at a realization
of the state $x=(x^{[1]},\dots,x^{[k-1]})$, the net value for a customer
joining queue $m=1,\dots,k-1$, is given by
${v}_{m}(x,y_{m})=\nu_{m}(x^{[m]}+y_{m})$, where $y_{m}$ is a realization of
$Y_{m}$. The value from balking is normalized to zero; ${v}_{k}(x)=0$.

Clearly, this is a generalization of the motivating example presented in
Section 1.1 and the unobservable M/G/1 model discussed in Example 1. The next
proposition presents sufficient conditions for convergence of the SA algorithm
to a Nash equilibrium.

######  Proposition 8.

Assume that for every $m\in\\{1,,\ldots,k-1\\}$: (i) $\lambda<\mu_{m}$; (ii)
$\mathrm{E}[Y_{m}^{4}]<\infty$; and (iii) the function
$\nu_{m}:\mathbb{R}_{+}\to\mathbb{R}$ is Lipschitz continuous. Then for any
initial strategy $\mathbf{p}^{(1)}\in\Delta$ and step-size sequence
$\\{\gamma_{n}\\}_{n\geq 1}$ satisfying Assumption A4, the SA algorithm
defined in (5) converges almost surely to a set of equilibrium strategies, in
the sense of Theorem 5.

###  4.2 Selective routing in queues with different buffer capacities

Consider a queueing network comprised of two servers, Server 1 that has no
waiting room, and Server 2 that offers unlimited queuing capacity. Customers
arrive at the system following a Poisson process with rate $\lambda$, with
i.i.d service-length requirements distributed according to $F$. Let $Y\sim F$
represent some arbitrary service duration, and let $\mathrm{E}[Y]=1/\mu$.

Upon arrival a customer cannot observe the system’s state, and is given the
option to _sense_ , that is, to make a costly attempt to obtain service by
Server 1, at a fixed cost $c_{s}$. A customer that chooses to sense and finds
Server 1 idle, immediately commences service, otherwise they are
instantaneously sent to the end of the infinite-capacity queue and waits to be
served by Server 2. The sensing option is offered only upon the customer’s
arrival, and the cost is incurred regardless of whether Server 1 is found
idling or not. Customers who do not choose the sensing option do not incur a
sensing cost and are sent directly to the end of the queue served by Server 2.
All customers queueing for Server 2 incur a linear waiting cost at rate
$c_{w}$ per unit time.

Let $\mathcal{A}=\\{a_{1},a_{2}\\}$, where $a_{1}$ represents the action
“sense” and $a_{2}$ the action “not sense”. A strategy in this game is given
by $\mathbf{p}=(p,1-p)$, where $p>0$ is the probability of choosing $a_{1}$.
It is shown in [38] that when customers’ service demand is exponentially
distributed, the game admits a unique equilibrium strategy which can be
efficiently approximated. Yet even for the special case of exponential
services, a simple expression for the equilibrium strategy is not available,
because the solution in [38] for the stationary distribution given a strategy
$\mathbf{p}$ makes use of Cardano’s formula. It is worth mentioning in this
context that in general, the queue served by Server 2 cannot be treated as a
standard M/G/1, because the effective arrival process to Server 2 is not
renewal, as it depends on the state of Server 1.

To approach the equilibrium using simulation we define the state as a pair
$X=(X^{[1]},X^{[2]})$ where each component $X^{[m]}$ stands for the workload
at Server $m\in\\{1,2\\}$. Given a realization $x=(x^{[1]},x^{[2]})$ of $X$,
the components of the vector function $\overline{\mathbf{v}}(x)$ take the form

$\begin{split}\overline{v}_{1}(x)&=-c_{s}-c_{w}\cdot\mathit{1}(x^{[1]}>0)\cdot
x^{[2]},\\\ \overline{v}_{2}(x)&=-c_{w}\cdot x^{[2]}.\end{split}$

######  Proposition 9.

Assume in the model of Section 4.2 that $\lambda<\mu$ and
$\mathrm{E}[Y^{4}]<\infty$. Then for any initial strategy
$\mathbf{p}^{(1)}\in\Delta$ and step-size sequence $\\{\gamma_{n}\\}_{n\geq
1}$ satisfying Assumption A4, the SA algorithm defined in (5) converges almost
surely to a set of equilibrium strategies, in the sense of Theorem 5.

The simplicity of the model here sets up a convenient platform to demonstrate
some practical heuristics that can enhance the algorithm’s performance, which
is what we focus on next.

###  4.3 Variance reduction using control variates

The variance reduction technique discussed here is a fundamental control-
variate method based on [64], that exploits our knowledge of the strategy
$\mathbf{p}$ at the beginning of each iteration, as well as some model
parameters used to construct the simulation. Adapted to our framework, the
goal is to modify the estimator $\mathbf{G}$ in a way that does not introduce
bias, but reduces its generalized variance. The underlying principle is the
following: Suppose one can characterize an $l$-dimensional random vector
$\mathbf{C}$, referred to as the control vector, with zero mean (entrywise),
such that it is also strongly correlated with $\mathbf{G}$. Then an
appropriate matrix of control coefficients
$\mathbf{\Psi}\in\mathbb{R}^{l\times k}$ can be found, for which the new
estimator $\mathbf{G}-\mathbf{\Psi}^{\prime}\mathbf{C}$ is unbiased for the
original unknown quantity, but has smaller generalized variance than that of
$\mathbf{G}$. Mathematically, this means that we would like to choose
$\mathbf{\Psi}$ such that $\mathbf{G}-\mathbf{\Psi}^{\prime}\mathbf{C}$
satisfies

$\mathrm{E}_{\mathbf{p}}[\mathbf{G}-\mathbf{\Psi}^{\prime}\mathbf{C}]=\mathrm{E}_{\mathbf{p}}[\mathbf{G}]=\mathbf{g(p)},\quad\text{and}\quad\mathrm{det}(\mathrm{Var}_{\mathbf{p}}[\mathbf{G}-\mathbf{\Psi}^{\prime}\mathbf{C}])<\mathrm{det}(\mathrm{Var}_{\mathbf{p}}[\mathbf{G}]).$

Observe that unbiasedness is obtained irrespective of the choice of
$\mathbf{\Psi}$. It is shown in [64] that the optimal choice for
$\mathbf{\Psi}$ (in terms of minimal variance) is given by
$\mathbf{\Psi}^{*}=\mathbf{\Sigma}_{\mathbf{GC}}\mathbf{\Sigma}^{-1}_{\mathbf{CC}}$,
where $\mathbf{\Sigma}_{\mathbf{CC}}=\mathrm{Var}_{\mathbf{p}}[\mathbf{C}]$ is
the variance-covariance matrix of the control vector $\mathbf{C}$, and
$\mathbf{\Sigma}_{\mathbf{GC}}=\mathrm{Cov}_{\mathbf{p}}(\mathbf{G,C})$ is the
cross-covariance matrix of $\mathbf{G}$ with $\mathbf{C}$. In most practical
applications, both $\mathbf{\Sigma}_{\mathbf{CC}}$ and
$\mathbf{\Sigma}_{\mathbf{GC}}$ are unknown, and so $\mathbf{\Psi}^{*}$ has to
be estimated from the data. Moreover, in the setting here the optimal
$\mathbf{\Psi}^{*}$ varies with the strategy $\mathbf{p}$. Luckily, our
empirical results show that significant variance reduction is obtained even
when using a rough approximation of $\mathbf{\Psi}^{*}$ which is based on
naive sampling. Our guiding intuition is that as the iteration index $n$
increases, the change in the strategy $\mathbf{p}^{(n)}$ becomes smaller.
Given a constant strategy, it is known that the estimator for
$\mathbf{\Psi}^{*}$ constructed from the sample covariance and cross-
covariance matrices is likelihood maximizing (see [4], Section 4.3.1).

At iteration $n>1$, let $\mathbf{G}^{(n)}$ be the estimator defined in (7),
and $\mathbf{C}^{(n)}$ be the control vector. We use a naive choice for the
matrix of control coefficients $\hat{\mathbf{\Psi}}^{(n)}$ by defining (note
that $\mathbf{C}^{(1)},\dots\mathbf{C}^{(n)}$ have mean zero):

$\hat{\mathbf{\Sigma}}^{(n)}_{\mathbf{CC}}=\frac{1}{n-1}\sum_{m=1}^{n}\mathbf{C}^{(m)}\left.\mathbf{C}^{(m)}\right.^{\prime},\quad\hat{\mathbf{\Sigma}}^{(n)}_{\mathbf{GC}}=\frac{1}{n-1}\sum_{m=1}^{n}\mathbf{G}^{(m)}\left.\mathbf{C}^{(m)}\right.^{\prime},\text{
and}\quad\hat{\mathbf{\Psi}}^{(n)}=\left.{\hat{\mathbf{\Sigma}}^{(n)}_{\mathbf{GC}}}\right.^{-1}\hat{\mathbf{\Sigma}}^{(n)}_{\mathbf{CC}}.$

Thus, at every iteration $n>1$ we replace $\mathbf{G}^{(n)}$ in Equation (5)
by
$\mathbf{G}^{(n)}-\left.\hat{\mathbf{\Psi}}^{(n)}\right.^{\prime}\mathbf{C}^{(n)}$.

######  Example 2.

Consider the model presented in Section 4.2. We choose a 3-dimensional control
vector $\mathbf{C}=(C_{1},C_{2},C_{3})$. As prescribed in Section 3.2, given a
strategy $\mathbf{p}=(p,1-p)$, we generate $L$ observations
$X_{1},\dots,X_{L}$ based on a single regeneration cycle, and construct the
original estimator $\mathbf{G}=\sum_{j=1}^{L}\mathbf{\overline{v}}(X_{j})$
introduced in (7). Let $D_{1},\dots,D_{L}$, be a sequence of _sensing
indicators_ , namely for $j=1,\dots,L$, $D_{j}\sim\text{Bernoulli}(p)$ is the
indicator of the event “the $j$-th customer chooses action $a_{1}$ (sense)”.
Let $Y_{j}$ be the service demand of that $j$-th customer. Our control vector
of choice, $\mathbf{C}$, is given by:

$\mathbf{C}=\begin{pmatrix}C_{1}\\\ C_{2}\\\
C_{3}\end{pmatrix}=\sum_{j=1}^{L}\begin{pmatrix}D_{j}-p\\\ Y_{j}-1/\mu\\\
\mathit{1}(X_{j}^{[1]}>0)-\frac{\lambda p}{\mu+\lambda p}\end{pmatrix}.$

Because the underlying process is regenerative, we have for every
$\mathbf{p}$,

$\mathrm{E}_{\mathbf{p}}[\mathbf{C}]=\ell(\mathbf{p})\cdot\begin{pmatrix}\mathrm{E}_{\mathbf{p}}[D]-p\\\
\mathrm{E}_{\mathbf{p}}[Y]-1/\mu\\\
\mathrm{E}_{\mathbf{p}}[\mathit{1}(X^{[1]}>0)]-\frac{\lambda p}{\mu+\lambda
p}\end{pmatrix}=\mathbf{0},$

where $D\sim\text{Bernoulli}(p)$ is the sensing indicator of an arbitrary
customer, $Y\sim F$ is an arbitrary service demand, and $X=(X^{[1]},X^{[2]})$
is the stationary system state, given strategy $\mathbf{p}=(p,1-p)$. It is
therefore immediate that $\mathrm{E}_{\mathbf{p}}[D]=p$ and
$\mathrm{E}_{\mathbf{p}}[Y]=1/\mu$, and because Server 1 can essentially be
modeled as an M/G/1/1 loss-system with arrival rate $\lambda p$ and service
rate $\mu$, we also have in stationarity that
$\mathrm{E}_{\mathbf{p}}[\mathit{1}(X^{[1]}>0)]=\frac{\lambda p}{\mu+\lambda
p}$. This specific choice of the control vector $\mathbf{C}$ is led by the
understanding that customers’ actions, their realized service demand, as well
as the state of Server 1 observed by the different customers during the
simulated cycle, are all affecting the utility sum $\mathbf{G}$ in that cycle.
$\diamond$

Figure 3 depicts the performance of the algorithm incorporating control
variates and dynamic step sizes (described below) for the model described in
Section 4.2. Further information about the efficiency and selection of control
coefficients for the method of control variates in vectored-output simulations
can be found in [64] and [76].

###  4.4 Dynamic step size selection

Next we consider the selection of the step size sequence,
$\\{\gamma_{n}\\}_{n\geq 1}$. Our main idea underlying the SA algorithm is
that instead of estimating $\mathbf{u(p)}$, to which unbiased estimators in
general are not known, one can estimate $\ell(\mathbf{p})\mathbf{u(p)}$ with
no bias, based on a single regeneration cycle, and use this estimator to
determine the progression of the algorithm. Yet, as noted in Remark 1, the
function $\ell(\mathbf{p})$ chosen is not exclusive; If
$h:\Delta\to\mathbb{R}_{+}$ is a real continuous function such that
$h(\mathbf{p})\geq a$ for some constant $a>0$, and in addition if one is
endowed with a finite-variance unbiased estimator for
$h(\mathbf{p})\mathbf{u(p)}$, then this estimator can be used instead of the
one suggested in Section 3, and the convergence result in Theorem 5 carries
over. The exact form of $h(\mathbf{p})$ need not be expressed or specified,
neither should it be introduced as input to the algorithm.

In theory, if $\ell(\mathbf{p})$ can be computed, then the estimator
$\mathbf{G}/\ell(\mathbf{p})$ can be taken to directly estimate
$\mathbf{u(p)}$ without bias. This is essentially equivalent to dynamically
choosing the step size at iteration $n$ to be $\gamma_{n}/\ell(\mathbf{p})$
instead of the original $\gamma_{n}$. The dynamic choice of strategy-dependent
step size has the advantage of “neutralizing”, in some sense, the unwanted
impact of the cycle lengths on the algorithm’s progression. Of course, in
applications, exact knowledge of $\ell(\mathbf{p})$ is an unrealistically-
strong requirement. However, even a coarse approximation
$\tilde{\ell}(\mathbf{p})\approx\ell(\mathbf{p})$ can impose a similar
neutralizing effect, leading to substantial efficiency improvement. Thus,
instead of using the estimator $\mathbf{G}$ to estimate
$\ell(\mathbf{p})\mathbf{u(p)}$, given an approximation
$\tilde{\ell}(\mathbf{p})$ we can use the estimator
$\mathbf{G}/\tilde{\ell}(\mathbf{p})$ to estimate
$(\ell(\mathbf{p})/\tilde{\ell}(\mathbf{p}))\mathbf{u(p)}$, or equivalently,
divide the step size by $\tilde{\ell}(\mathbf{p})$. We believe this method to
be effective especially in heavily loaded systems, in which the mean cycle
length is typically large, and in addition, many approximation techniques tend
to become more accurate under this regime.

######  Remark 3.

Natural approximations for ${\ell}(\mathbf{p})$ can be derived from available
approximations for mean busy periods in queues. For instance, many
approximations, in addition to some special-case exact expressions, are
available for the mean busy period in the GI/G/1 queue (see, for instance,
[11]). The length of an arbitrary regeneration cycle $L$ of a GI/G/1 can be
interpreted as the number of service completions during the corresponding busy
period $B$. Thus, the relation between $L$ and $B$ is expressed by
$B=\sum_{j=1}^{L}Y_{j}$, where $Y_{1},\dots,Y_{L}$ are the service lengths
during the cycle. This relation implies $\mathrm{E}[L]=\mathrm{E}[B]/E[Y]$,
where $Y$ is the length of an arbitrary service.

######  Example 2.

(Continued) In the model of Section 4.2, a back-of-the-envelope approximation
for $\ell(\mathbf{p})$ is obtained by (erroneously) treating the two nodes,
Server 1 and 2, as two independent Markovian systems with service rate $\mu$
in both systems. Thus, given $\mathbf{p}=(p,1-p)$, we have one M/M/1/1 system
with arrival rate $\lambda^{[1]}=p\lambda$, and another independent M/M/1
system with arrival rate
$\lambda^{[2]}=\lambda\cdot(1-p+\lambda^{[1]}/(\lambda^{[1]}+\mu))$. The
stationary probability of finding each server idling is
$\pi_{0}^{[1]}=\mu/(\mu+\lambda^{[1]})$ in the M/M/1/1 system and
$\pi_{0}^{[2]}=1-\lambda^{[2]}/\mu$ in the M/M/1 system, and by independence
the joint distribution of the events is $\pi_{0}^{[1]}\pi_{0}^{[2]}$. By the
PASTA principle, this coincides with the probability that a customer finds
both servers idle upon arrival, therefore the mean number of arrivals between
two subsequent such incidents is given by
$\tilde{\ell}(\mathbf{p})=1/(\pi_{0}^{[1]}\pi_{0}^{[2]})$, and we shall use
this quantity to approximate the true mean cycle length $\ell(\mathbf{p})$.
$\diamond$

In the following numerical example, we employ variance reduction and dynamic
step-size selection on the model described in Section 4.2. Services are
exponentially distributed, $\mu$ and $c_{w}$ are normalized to 1,
$\lambda=0.99$ and $c_{s}=5$. The initial strategy for the SA algorithm is
$\mathbf{p}^{(1)}=(1/2,1/2)$. In Figure 3 we compare, over the first
$N=5\times 10^{4}$ iterations, the performance of the crude implementation of
the algorithm with that of the refined version. For the sequence of
strategies, $\\{\mathbf{p}^{(n)}\\}_{n\geq
1}=\\{({p}^{(n)},1-{p}^{(n)})\\}_{n\geq 1}$, the figure depicts the
convergence of the sequence ${p}^{(n)}$ to the (unique) equilibrium
probability ${p}^{e}$, that corresponds with the unique equilibrium strategy
$\mathbf{p}^{e}=(p^{e},1-p^{e})$. The estimator used at iteration $n$ in the
refined version (orange) is
$\mathbf{G}^{(n)}-\left.\hat{\mathbf{\Psi}}^{(n)}\right.^{\prime}\mathbf{C}^{(n)}$,
and the step size at iteration $n$ is $(1/n)\tilde{\ell}(\mathbf{p}^{(n)})$.
In the crude version (blue), we keep our original estimator
$\mathbf{G}^{(n)}$, with comparable, non-dynamic step size given by
$(1/n)\tilde{\ell}(\mathbf{p}^{(1)})$ for each $n$. For reference, we plot the
correct equilibrium (red dashed horizontal line), which, for exponential
services can be approached numerically up to any arbitrary precision level
using the machinery developed in [38].

Figure 3: Convergence of the SA algorithm as $n\to\infty$. The first
coordinate of $\mathbf{p}^{(n)}=({p}^{(n)},1-{p}^{(n)})$, i.e., the sensing
probability at stage $n$, is plotted vs. $n$. The blue curve corresponds with
the crude algorithm, and the orange corresponds with the refined version. The
red dashed line depicts the correct equilibrium sensing probability, $p^{e}$,
computed numerically in 8-digit precision.

##  5 Games with state information

In this section, we explain how to amend our SA algorithm to cope with so-
called observable models. More precisely, we assume that customers, at the
moment they arrive to the system and make their decision, have access to an
information signal, which is a function of the system state. Unlike the
approach described in 3.5, the one we take here does not suffer from the same
issue of dimensionality.

From a mathematical perspective, the previous formulation described in Section
2 can be seen as a special case of the model below when the function mapping
states to signals is constant (i.e., non-informative). Although, from a
nomenclature perspective, this understanding is somewhat inconsistent with how
the terms ’observable’ and ’unobservable’ are used in the literature to
describe two disjoint classes of games. In what follows, we treat a strategy
as function from signals to (mixed) actions, and so, the dimension of our
problem grows proportionally with the number of possible signals. Thus, it
will be assumed throughout that the number of possible signals is finite.

The new formulation suggests a practical method for the computation of
equilibrium strategies in the observable GI/G/1 queue, a generalization of
Naor’s [59] model, which we use as a guiding example in this section. The
special case with Poisson arrivals, which itself is a generalization of [59],
was thoroughly studied in [47].

###  5.1 Model description

As before, we assume a renewal sequence of customer arrivals
$\\{A_{n}\\}_{n\geq 1}$, each customer chooses an action from the set
$\mathcal{A}=\\{a_{1},\dots,a_{k}\\}$. The state of the system at any moment
in time is represented by a vector $x$ of buffer contents in the state space
$\mathcal{X}\subset\mathbb{R}_{+}^{d}$, in which $0^{d}\in\mathcal{X}$ is an
element. Upon arrival, a customer is provided with an information signal which
depends on the state at which they arrive. Let $\mathcal{I}$ be a countable
set of information signals, thus, arriving at state $x\in\mathcal{X}$, a
customer receives information $\psi(x)$, where
$\psi:\mathcal{X}\to\mathcal{I}$ is a mapping from states to signals. We
therefore reintroduce the concept of a strategy, now taking the interpretation
of a _behavioral strategy_ , i.e., a mapping $\mathbf{p}:\mathcal{I}\to\Delta$
from signals to distributions over actions. We denote by $\mathcal{P}$ the set
of possible strategies. Hence, given $\mathbf{p}\in\mathcal{P}$, a customer
who arrives at state $x\in\mathcal{X}$ will choose their action according to
the distribution
$\mathbf{p}(\psi(x))=\big{(}p_{1}(\psi(x)),\dots,p_{k}(\psi(x))\big{)}\in\Delta$,
which assigns a probability $p_{j}(\psi(x))$ for each action
$a_{j}\in\mathcal{A}$. Consistent with the previous sections, we assume each
strategy $\mathbf{p}\in\mathcal{P}$ gives rise to a regenerative system-state
process $\\{X_{n}(\mathbf{p})\\}_{n\geq 1}$ embedded at arrival instants, with
$0^{d}$ being a regenerative state and $X_{1}(\mathbf{p})=0^{d}$. Thus, given
$\mathbf{p}\in\mathcal{P}$, the stationary state $X(\mathbf{p})$, the cycle
length $L(\mathbf{p})$, and the $r$’th moment of the cycle length
$\ell^{r}(\mathbf{p})$, all follow the definitions introduced in Section 2.
Consistent with Sections 2 and 3, when facing state $x\in\mathcal{X}$ (hence a
signal $\psi(x)$ is observed) and a realization $y$ of a random outcome $Y$,
the value vector is
$\mathbf{v}(x,y)=\big{(}v_{1}(x,y),\dots,v_{k}(x,y)\big{)}$, and we let
$\mathbf{\overline{v}}(X)=\mathrm{E}_{Y}\left[\mathbf{v}\big{(}X,Y\big{)}\>\big{|}\>X\right]$.

Assuming the state is drawn from the stationary distribution induced by
$\mathbf{p}$, the signal observed is given by the r.v. $\psi(X(\mathbf{p}))$
supported on $\mathcal{I}$. We denote the (point) probability function of
$\psi(X(\mathbf{p}))$ by
$\xi_{\mathbf{p}}(s)=\mathrm{P}_{\mathbf{p}}(\psi(X(\mathbf{p}))=s)$ and its
support by $\mathcal{I}(\mathbf{p})=\\{s\in\mathcal{I}\text{ s.t.
}\xi_{\mathbf{p}}(s)>0\\}$. For conciseness, we ignore any irrelevant signals,
i.e., we assume
$\mathcal{I}=\bigcup_{\mathbf{p}\in\mathcal{P}}\mathcal{I}(\mathbf{p})$. When
$s\in\mathcal{I}(\mathbf{p})$, we say that the signal $s$ is ‘recurrent’ for
$\mathbf{p}$, otherwise we call it ‘transient’.

For any recurrent signal $s\in\mathcal{I}(\mathbf{p})$ we define, analogous to
(2), the ($k$-dimensional) _conditional_ expected value vector as

$\mathbf{u}(\mathbf{p}\mid
s)=\mathrm{E}_{\mathbf{p}}\bigg{[}\mathbf{v}\big{(}X(\mathbf{p}),Y\big{)}\>\bigg{|}\>\psi(X(\mathbf{p}))=s\bigg{]}=\frac{\mathrm{E}_{\mathbf{p}}\bigg{[}\mathbf{v}\big{(}X(\mathbf{p}),Y\big{)}\cdot\mathit{1}(\psi(X(\mathbf{p}))=s)\bigg{]}}{\xi_{\mathbf{p}}(s)}.$
(13)

For completeness, when $s\notin\mathcal{I}(\mathbf{p})$ we let
$\mathbf{u}(\mathbf{p}\mid s)=\mathbf{0}$. Finally, the best response set for
a strategy $\mathbf{p}$ is given by

$\mathcal{BR}(\mathbf{p})=\operatorname*{arg\,max}_{\mathbf{q}\in\mathcal{P}}\sum_{s\in\mathcal{I}(\mathbf{p})}\xi_{\mathbf{p}}(s)\cdot\mathbf{u}(\mathbf{p}\mid
s)^{\prime}\mathbf{q}(s),$ (14)

and similarly to Definition 2, we have that $\mathbf{p}$ is a symmetric Nash
equilibrium if $\mathbf{p}\in\mathcal{BR}(\mathbf{p})$. The objective in the
optimization problem of (14) can be thought of as the ex-ante expected utility
(prior to receiving state information) of a customer who plays $\mathbf{q}$,
when all other customers play $\mathbf{p}$. It can be further seen that the
optimization problem is separable, in the sense that it can be solved
independently for every $s\in\mathcal{I}(\mathbf{p})$. In other words, to
constitute a best response, a strategy has to prescribe an optimal play,
separately, for any state information observed.

###  5.2 Simulation and stochastic approximation

Again, we construct our utility estimator by simulating a single regenerative
cycle starting with $X_{1}(\mathbf{p})=0^{d}$, with $X_{j}$, $j=1,\dots,L$
being the state at the $j$-th arrival, and $L$ being the cycle length. To
apply the SA algorithm we construct an estimator for every $s\in\mathcal{I}$
defined by

$\mathbf{G}(s)=\sum_{j=1}^{L}\mathbf{\overline{v}}(X_{j})\cdot\mathit{1}(\psi(X_{j})=s).$

Assuming $\ell^{2}(\mathbf{p})<\infty$, an immediate extension of Lemma 4
yields, for all $s\in\mathcal{I}$,

$\mathrm{E}_{\mathbf{p}}[\mathbf{G}(s)]=\ell(\mathbf{p})\cdot\xi_{\mathbf{p}}(s)\cdot\mathbf{u}(\mathbf{p}\mid
s).$ (15)

Note that if $s\notin\mathcal{I}(\mathbf{p})$, then $\xi_{\mathbf{p}}(s)=0$,
and equality in (15) trivially holds as we have that $\mathbf{G}(s)=0$ with
probability 1.

In the case with state information, our algorithm constructs at each iteration
a point estimate $\mathbf{G}(s)$ for every $s\in\mathcal{I}$. Thus, for
practical reasons we make the assumption hereafter that $\mathcal{I}$ is a
finite set. In many applications, a trivial bound on $|\mathcal{I}|$ can be
derived, as demonstrated in Example 3 below. In other applications, different
heuristics can be applied to cap the number of signals, with arguably
“minimal” effect on the best-response function, however such considerations
exceed the scope of the this paper.

The progression of the SA algorithm is similar to the one described in the
previous sections: Let the step-size sequence be $\\{\gamma_{n}\\}_{n\geq 1}$
and the initial strategy be $\mathbf{p}^{(1)}\in\mathcal{P}$. At iteration
$n$, with a corresponding strategy $\mathbf{p}^{(n)}$, we use one regeneration
cycle to construct $|\mathcal{I}|$ estimators,
$\\{\mathbf{G}^{(n)}(s)\\}_{s\in\mathcal{I}}$. The update at the next
iteration for every $s\in\mathcal{I}$ is done by

$\mathbf{p}^{(n+1)}(s)=\pi_{\Delta}\left(\mathbf{p}^{(n)}(s)+\gamma_{n}\mathbf{G}^{(n)}(s)\right),\
n\geq 0,$ (16)

The following assumptions are sufficient for the convergence of the algorithm
to an equilibrium strategy:

######  Assumption B1.

The second moment of the cycle length $L(\mathbf{p})$ is bounded for every
strategy; $\ell^{2}(\mathbf{p})<\infty$ for all $\mathbf{p}\in\mathcal{P}$.

######  Assumption B2.

The second moment of the estimator $\mathbf{G}(s)$ is uniformly bounded on
$\mathcal{P}$;
$\sup_{\mathbf{p}\in\mathcal{S}}\mathrm{E}_{\mathbf{p}}\|\mathbf{G}(s)\|^{2}<\infty$
for all $s\in\mathcal{I}$.

######  Assumption B3.

For all $s\in\mathcal{I}$,
$\xi_{\mathbf{p}}(s)\cdot\mathbf{H}\mathbf{u}(\mathbf{p}\mid s)$ is continuous
with respect to $\mathbf{p}$, and there exists a function
$u^{*}:\mathcal{P}\to\mathbb{R}$ such that for all $\mathbf{p}\in\mathcal{P}$
and $s\in\mathcal{I}$,
$-\xi_{\mathbf{p}}(s)\cdot\mathbf{H}\mathbf{u}(\mathbf{p}\mid
s)=\nabla_{\mathbf{p}(s)}u^{*}$.

Theorem 10 introduces the key result of the section, which is a generalization
of Theorem 5:

######  Theorem 10.

Suppose Assumptions B1–B3 and A4 are satisfied. As $n\to\infty$,
$\mathbf{p}^{(n)}\>\to_{\rm as}\>\mathcal{S}^{e}$, where
$\mathcal{S}^{e}\subseteq\mathcal{P}$ is a set of equilibrium strategies,
namely, $\mathbf{p}\in\mathcal{BR}(\mathbf{p})$ for every
$\mathbf{p}\in\mathcal{S}^{e}$.

The proof of Theorem 10 is similar to that of Theorem 5, after introducing a
definition of the the surrogate best-response function, analogous to the one
defined in Equation (3). This is described in further elaboration in Appendix
D. Noticeably, Theorem 10 also implies that an equilibrium strategy exists. In
fact, if one only cares to prove that an equilibrium strategy exists, then
only the continuity of
$\xi_{\mathbf{p}}(s)\cdot\mathbf{H}\mathbf{u}(\mathbf{p}\mid s)$ is needed
(see Appendix D).

As in Section 3.4, the existence of a potential function $u^{*}$, as referred
to in Assumption B3, ensures that (with probability 1) in the limit, the
iterate will not oscillate between non-equilibrium points. In the absence of
Assumption B3, even a non-stochastic variant of the algorithm can exhibit such
undesirable oscillations in the limit. Although the implementation of the
algorithm does not require $u^{*}$ to be specified, verifying Assumption B3
can be difficult in practical settings. It is nevertheless the case that if
the limit $\lim_{n\to\infty}\mathbf{p}^{(n)}$ exists then it must satisfy the
equilibrium condition, regardless of whether Assumption B3 holds or not.

######  Example 3.

Consider a first-come first-served GI/G/$1$ queue with observable queue
lengths: The state $x=(x^{[1]},x^{[2]},\dots)$ is a vector of residual service
times, where $x^{[i]}$ is the residual service of the $i$’th customer in line,
which, for $i=1$, is that of the one in service. For simplicity, we assume
here that service times are strictly positive iid random variables with mean
$1/\mu<\infty$. Customer utility when joining the queue is linear in the
waiting time, and is equal to zero when balking. Formally,
$\mathbf{v}(x,y)=(v_{1}(x,y),v_{2}(x,y))=(R-C\cdot(y+\|x\|_{1}),0)$ for
positive constants $R,C>0$, implying that for a random state $X$,
$\overline{\mathbf{v}}(X)=(R-C\cdot(1/\mu+\|X\|_{1}),0)$.

The information available to a customer arriving at state $x$ is the number of
jobs in the system, $\|x\|_{0}=|\\{i\text{ s.t. }x^{[i]}>0\\}|$. Importantly,
following the FCFS assumption, any customer not in service has an expected
residual service time of $1/\mu$. Hence, from an arriving customer’s stand
point, the (state-dependent) utility function $\overline{\mathbf{v}}$ can be
further simplified into
$\overline{\mathbf{v}}(X)=(R-C\cdot(\|X\|_{0}/\mu+X^{[1]}),0)$. It is
therefore immediate that when $\|x\|_{0}\geq R\mu/C$, an arriving customer
strictly prefers balking over joining. By defining $K=\lfloor R\mu/C\rfloor$
we can thus _assume_ that customers balk when observing $K+1$ (or more)
customers in queue, and we restrict attention to decision epochs, which are
the arrival instants of customers to states in $\mathcal{X}=\\{x\text{ s.t.
}\|x\|_{0}\leq K\\}$. This implies that at decision epochs, the state-
information function $\psi(x)=\|x\|_{0}$ defines a finite set of possible
information signals (queue lengths), $\mathcal{I}=\\{0,\dots,K\\}$. A strategy
$\mathbf{p}$ is a mapping from $\mathcal{I}$ to the two-dimensional simplex,
$\mathbf{p}(s)=(p(s),1-p(s))\in\Delta$, with $p(s)\in[0,1]$ being the
probability of joining when $s\in\mathcal{I}$ customers are observed. Hence,
$\xi_{\mathbf{p}}(s)$ is the stationary probability of observing
$s\in\mathcal{I}$ customers in queue at a decision epoch, provided all
customers play according to $\mathbf{p}$.

Given $\mathbf{p}$, the dynamics of the process imply that if customers balk
with probability $1$ at queue length ${s}$, i.e., if $p(s)=0$, then any queue
length $\tilde{s}>s$ is transient. Put formally, we have for all
$\mathbf{p}\in\mathcal{P}$, $\mathcal{I}(\mathbf{p})=\\{s\in\mathcal{I}\text{
s.t }\forall\tilde{s}<s,p(\tilde{s})>0\\}$. An equilibrium, adapted to this
example, is a strategy $\mathbf{p}\in\mathcal{P}$ such that for all
$s\in\mathcal{I}(\mathbf{p})$, if $u_{1}(\mathbf{p}\mid s)>0$ then $p(s)=1$,
and if $u_{1}(\mathbf{p}\mid s)<0$ then $p(s)=0$. We note that if $\mathbf{p}$
is an equilibrium, there could be at most one $s\in\mathcal{I}(\mathbf{p})$
satisfying $u_{1}(\mathbf{p}\mid s)<0$, which also constitutes a threshold on
the length of the queue (as long as customers play according to $\mathbf{p}$).

Figure 4 depicts a simulation in which we set $C=1$ and $R=1.7$, implying that
$K=1$ and $\mathcal{I}=\\{0,1\\}$. Thus, every customer that observes two (or
more) other customers automatically balks. In both panels, interarrival times
are exponentially distributed with mean $1$, which is equal to the mean
service time, yet the service distribution varies between the two panels. In
panel (a), services are exponentially distributed, hence it is a particular
case of Naor’s observable M/M/1 model ([59]) in which the unique equilibrium,
$\\{\mathbf{p}^{e}(s)\\}_{s\in\\{0,1\\}}=\\{(p^{e}(s),1-p^{e}(s))\\}_{s\in\\{0,1\\}}$
is easily found, with $p^{e}(0)=1$ and $p^{e}(1)=0$, i.e., join if and only if
the system is empty. In panel (b), services are uniformly distributed over
$[0,2]$. We note that the uniform distribution is of a decreasing mean
residual lifetime, and by Proposition 3.1 in Kerner ([47]) the equilibrium
exists uniquely. In both experiments we run the simulation for $N=10^{5}$
iterations of the algorithm, with $\gamma_{n}=2/n$ and initial strategy
$\mathbf{p}^{(1)}(s)=(1/2,1/2)$ for all $s\in\\{0,1\\}$. Interestingly, in the
right panel, the equilibrium to which the algorithm converges is of a
_randomized threshold_ type; when a customer observes another customer in
service, they join with probability $p^{(N)}(1)\approx 0.37$.

(a) exponential service

(b) uniform service

Figure 4: The convergence of the SA algorithm as $n\to\infty$. The joining
probability $p^{(n)}(s)$ is plotted vs. $n$ on a square-root scale, for
$s=0,1$. In both panels, the blue curve depicts ${p}^{(n)}(0)$ and the green
curve depicts ${p}^{(n)}(1)$.

######  Remark 4.

It is possible that some information signals are transient under some
strategies but recurrent under different strategies. For instance, in the
observable GI/G/1 queue, any strategy $\mathbf{p}(s)=(p(s),1-p(s))$ with
$\overline{s}\in\mathcal{I}(\mathbf{p})$ such that $p(\overline{s})=0$
(meaning customers balk when observing queue length $\overline{s}$), imply
that the probability of observing a queue longer than $\overline{s}$ is zero,
i.e., for all $s_{0}>\overline{s}$, $s_{0}\notin\mathcal{I}(\mathbf{p})$. As a
matter of fact, the equilibrium condition
$\mathbf{p}\in\mathcal{BR}(\mathbf{p})$ does not impose any restriction on
customer decisions in the (null) event they receive a zero-probability signal.
Hassin and Haviv [36] introduce an adaptation of the concept of sub-game
perfect Nash equilibria to queueing games, which is a refinement of the
equilibrium concept, adding the requirement that customers play optimally also
when arriving to transient states. In our setup, convergence to a sub-game
perfect equilibrium cannot be guaranteed: Assuming the algorithm converges to
an equilibrium $\mathbf{p}^{e}$ with $s_{0}\notin\mathcal{I}(\mathbf{p}^{e})$,
the $s_{0}$-th (vector) component of the strategy will be updated, if at all,
only finitely often; the tail of the sequence
$\\{\mathbf{p}^{(n)}(s_{0})\\}_{n\geq 0}$ is a constant vector (almost
surely); and the limiting strategy $\mathbf{p}^{e}$ is possibly suboptimal at
$s_{0}$.

##  6 Concluding remarks

This work introduces a robust stochastic-approximation algorithm that computes
symmetric Nash equilibria in a general class of queueing games. The method
involves simulating the system once and updating the strategy at regeneration
times using the Robbins-Monro algorithm. As opposed to classical stochastic-
optimization settings, a Nash equilibrium of a game with a discrete action
space is not given by a first-order gradient condition. Therefore, an
equivalent condition for a Nash equilibrium is formulated as the root of a
continuous function. This enables the construction of an iterative fixed-point
method resembling those used in stochastic approximation. Results from renewal
theory are used in order to construct an unbiased estimator of the total
utility observed during a regeneration period. Verifiable conditions for
almost-sure convergence of the algorithm are further provided. The method is
shown to be useful for various interesting applications.

We suspect that the suggested algorithm may be further used in order to find a
socially-optimal strategy, by considering the social cost of deviation instead
of the individual optimality condition (see Haviv and Oz [42]). Loosely
speaking, the concept of social cost of deviation basically allows one to
reformulate the problem of identifying socially-optimal strategies as a
problem of finding equilibrium strategies in a game, and therefore deploy the
SA algorithm to approximate these strategies.

The framework developed here can be extended and improved in several ways, in
particular, by introducing more accurate estimators for the utility function.
A direct extension is to update the strategy after multiple regeneration
cycles. Having several observations of regeneration cycles will reduce the
variance of the estimator and will improve the accuracy of the iteration
steps. Of course, this method entails an inherent trade-off because the
algorithm can potentially “waste time” on simulating the process for
strategies that are far off from the correct solution. Thus, a suitable,
perhaps dynamic choice of the simulation length at each iteration should be
implemented.

The regeneration cycle approach comes along with both theoretical and
practical advantages – its implementation is straightforward and the estimator
for $\mathbf{g}(\mathbf{p}^{(n)})$ is independent and unbiased. In this
setting, the theory of stochastic approximations provides verifiable
sufficient conditions for almost-sure convergence of the algorithm. An
interesting future research direction is to explore the applicability of the
SA framework to a more flexible iterative scheme that updates the strategy at
every arrival and not just at regeneration times. This can be seen as a
reinforcement-learning approach, although the resulting process is clearly not
Markov. Implementation of such an algorithm should be just as easy, but
theoretical analysis is much more involved. In particular, the estimation step
is biased and is not independent of previous iterations. The theory on SA
algorithms with correlated noise may be useful for constructing meaningful
convergence conditions in this setting (see [49, Ch. 6]). This modification is
appealing because it has the potential of coping with cases in which some
strategies drive the system out of the stability region. Proving that such a
scheme of simple sequential interactions among agents eventually converges to
an equilibrium will further strengthen the plausibility of equilibrium
emergence in realistic systems and highlight the importance of studying
equilibrium analysis in queues.

##  7 Acknowledgements

The authors are grateful to Refael Hassin and Moshe Haviv for their advice and
comments on the paper, as well as to Jim Dai and Binyamin Oz for several
fruitful discussions. This research was supported in part by the Shenzhen
Research Institute for Big Data International Postdoctoral Fellowship.

## References

  * [1] I. Adiri and U. Yechiali (1974).  Optimal priority-purchasing and pricing decisions in nonmonopoly and monopoly queues.  Operations Research, 22(5): 1051–1066.
  * [2] E. Altman and Z. Altman (2003).  S-modular games and power control in wireless networks.  IEEE Transactions on Automatic Control, 48(5):839–842.
  * [3] E. Altman and N. Shimkin (1998).  Individual equilibrium and learning in processor sharing systems.  Operations Research, 46(6):776–784.
  * [4] T. W. Anderson (1962).  An introduction to multivariate statistical analysis, 3rd edition. Wiley, New York.
  * [5] S. Andradóttir (1996).  Optimization of the transient and steady-state behavior of discrete event systems.  Management Science, 42(5):717–737.
  * [6] S. Asmussen (2003).  Applied Probability and Queues, 2nd edition. Springer, New York.
  * [7] S. Asmussen and P. W. Glynn (2007).  Stochastic simulation: algorithms and analysis, Springer Science & Business Media.
  * [8] O. Baron, X. Chen and Y. Li (2022).  Omnichannel services: The false premise and operational remedies.  Management Science.
  * [9] C.E. Bell and S. Jr. Stidham (1983).  Individual versus social optimization in the allocation of customers to alternative servers.  Management Science, 29(7): 831-839.
  * [10] I. Ben-Shahar, A. Orda and N. Shimkin. (2000).  Dynamic service sharing with heterogeneous preferences.  Queueing Systems, 35(1):83–103.
  * [11] D. J. Bertsimas, D. Nakazato. (1992).  Transient and busy period analysis of the GI/G/1 queue: the method of stages.  Queueing Systems, 10(3): 153–184.
  * [12] V.S. Borkar (2008).  Stochastic approximation: a dynamical systems viewpoint, Cambridge University Press.
  * [13] R. Buche and H. J. Kushner (2000).  Stochastic approximation and user adaptation in a competitive resource sharing system.  IEEE Transactions on Automatic Control, 45(5):844–853.
  * [14] R. Buche and H. J. Kushner (2002).  Rate of convergence for constrained stochastic approximation algorithms.  SIAM journal on control and optimization, 40(4):1011–1041.
  * [15] A. Burnetas and A. Economou (2007).  Equilibrium customer strategies in a single server Markovian queue with setup times.  Queueing Systems, 56(3): 213–228.
  * [16] F. Castro, P. Frazier, H. Ma, H. Nazerzadeh and C. Yan (2020).  Matching Queues, Flexibility and Incentives.  Available at SSRN, https://ssrn.com/abstract=3627920
  * [17] M. Chen, M. Hu and J. Wang (2022).  Food delivery service and restaurant: Friend or foe?  Management Science.
  * [18] X. Chen, Y. Liu and G. Hong (2020).  An online learning approach to dynamic pricing and capacity sizing in service systems.  arXiv preprint, arXiv:2009.02911.
  * [19] M. A. Crane and D.L. Iglehart (1975).  Simulating stable stochastic systems: III. Regenerative processes and discrete-event simulations.  Operations Research, 23(1):33–45.
  * [20] S. Cui, Z. Wang and L. Yang (2020).  The economics of line-sitting.  Management Science, 66(1): 227–242.
  * [21] J. G. Dai and M. Gluzman (2021).  Queueing Network Controls via Deep Reinforcement Learning.  Stochastic Systems.
  * [22] D. J. Daley and D.R. Jacobs (1969).  The total waiting time in a busy period of a stable single-server queue, II.  Journal of Applied Probability, 6(3):565–572.
  * [23] C. Daskalakis, A. Mehta and C. Papadimitriou (2009).  TA note on approximate Nash equilibria.  Theoretical Computer Science, 410(17):1581–1588.
  * [24] De Meyer, A., and J. L. Teugels (1980).  On the asymptotic behaviour of the distributions of the busy period and service time in M/G/1.  Journal of Applied Probability, 17(3): 802–813.
  * [25] L.G. Debo, L.B. Toktay and L.N. Van Wassenhove (2008).  Queuing for expert services.  Management Science, 54(8): 1497–1512.
  * [26] N. M. Edelson and D. K. Hilderbrand (1975).  Congestion tolls for poisson queuing processes.  Econometrica, 43(1), 81–92.
  * [27] M. C. Fu (1990).  Convergence of a stochastic approximation algorithm for the GI/G/1 queue using infinitesimal perturbation analysis.  Journal of Optimization Theory and Applications, 65(1):149–160.
  * [28] M. C. Fu (1994).  Optimization via simulation: A review.  Annals of operations research, 53(1):199–247.
  * [29] D. Gamarnik (2002).  On deciding stability of constrained homogeneous random walks and queueing systems.  Mathematics of Operations Research, 27(2):272–293.
  * [30] A. Gosh, A. Bassamboo and M. Larviere (2020).  The Queue Behind the Curtain: Information Disclosure in Omnichannel Services.  Available at SSRN, https://ssrn.com/abstract=3730482
  * [31] P. Guo and R. Hassin (2001).  Strategic behavior and social optimization in Markovian vacation queues.  Operations research, 59(4): 986–997.
  * [32] P. Guo and P. Zipkin (2007).  Analysis and comparison of queues with different levels of delay information.  Management Science, 53(6): 962–970.
  * [33] R. Hassin (1996).  On the advantage of being the first server.  Management Science, 42(4): 618–623.
  * [34] R. Hassin. (2016). Rational Queueing, CRC Press.
  * [35] R. Hassin and M. Haviv (1997).  Equilibrium threshold strategies: The case of queues with priorities.  Operations Research, 45(6): 966–973.
  * [36] R. Hassin and M. Haviv. (2002).  Nash equilibrium and subgame perfection in observable queues.  Annals of Operations Research, 113(1):15–26.
  * [37] R. Hassin and M. Haviv. (2003). To Queue or Not to Queue: Equilibrium Behavior in Queueing Systems, Springer.
  * [38] R. Hassin and R.I. Snitkovsky (2017).  Strategic customer behavior in a queueing system with a loss subsystem.  Queueing Systems, 86(3):361–387.
  * [39] R. Hassin and R.I. Snitkovsky (2021).  On the Supermarket Game with finitely many servers.  Working paper.
  * [40] R. Hassin and R. Roet-Green (2017).  The impact of inspection cost on equilibrium, revenue, and social welfare in a single-server queue.  Operations Research, 65(3): 804–820.
  * [41] R. Hassin and R. Roet-Green (2021).  On Queue-Length Information when Customers Travel to a Queue.  Manufacturing & Service Operations Management, 23(4): 989–1004.
  * [42] M. Haviv and B. Oz (2021).  A busy period approach to some queueing games.  Queueing Systems, 97(3):261–77.
  * [43] M. Haviv and L. Ravner (2021).  A survey of queueing systems with strategic timing of arrivals.  Queueing Systems, 99(1), 163–198.
  * [44] M. Hu, J. Wang and H. Wen (2020).  Share or solo? Individual and social choices in ride-hailing.  Available at SSRN, https://ssrn.com/abstract=3675050
  * [45] J. Jacob and R. Roet-Green (2021).  Ride solo or pool: Designing price-service menus for a ride-sharing platform.  European Journal of Operational Research, 295(3): 1008–1024.
  * [46] K. Jagannathan, I. Menache, E. Modiano and G. Zussman (2012).  Non-cooperative spectrum access—The dedicated vs. free spectrum choice.  IEEE Journal on Selected Areas in Communications, 30(11): 2251–2261.
  * [47] Y. Kerner (2011).  Equilibrium joining probabilities for an M/G/1 queue.  Games and Economic Behavior, 71(2):521–526.
  * [48] J. Kiefer and J. Wolfowitz (1952).  Stochastic estimation of the maximum of a regression function.  The Annals of Mathematical Statistics, 23(3):462–466.
  * [49] H. Kushner and G. G. Yin (2003). Stochastic approximation and recursive algorithms and applications. Springer Science & Business Media.
  * [50] P. L’Ecuyer and P. W. Glynn (1994).  Stochastic optimization by simulation: Convergence proofs for the GI/G/1 queue in steady-state.  Management Science, 40(11):1562–1578.
  * [51] P. L’Ecuyer, N. Giroux and P. W. Glynn (1994).  Stochastic optimization by simulation: Numerical experiments with the M/M/1 queue in steady-state.  Management Science, 40(10):1245–1261.
  * [52] B. Liu, Q. Xie and E. Modiano (2019).  Reinforcement learning for optimal control of queueing systems.  2019 57th Annual Allerton Conference on Communication, Control, and Computing (Allerton), 663–670. IEEE
  * [53] M. Mandjes, B. Patch, and N. S. Walton. (2019).  Detecting Markov chain instability: a monte carlo approach.  Stochastic Systems, 7(2):289–314.
  * [54] C. Maglaras and A. Zeevi (2003).  Pricing and capacity sizing for systems with shared resources: Approximate solutions and scaling relations.  Management Science, 49(8): 1018–1038.
  * [55] A. Manou, P.G. Canbolat and F. Karaesmen (2017).  Pricing in a transportation station with strategic customers.  Production and Operations Management, 26(9): 1632–1645.
  * [56] A. Manou, A. Economou and F. Karaesmen (2014).  Strategic customers in a transportation station: when is it optimal to wait?.  Operations Research, 62(4): 910–925.
  * [57] H. Mendelson (1985).  Pricing computer services: Queueing effects.  Communications of the ACM, 28(3): 312–321.
  * [58] M. Maschler, E. Solan and S. Zamir. (2013). Game Theory, Cambridge University Press.
  * [59] P. Naor. (1969).  The regulation of queue size by levying tolls.  Econometrica, 37(1):15–24.
  * [60] A. Nemirovski, A. Juditsky, G. Lan and A. Shapiro (2009).  Robust stochastic approximation approach to stochastic programming.  SIAM Journal on optimization, 19(4):1574–1609. 37(1):15–24.
  * [61] A. K. Parlaktürk and S. Kumar (2004).  Self-interested routing in queueing networks.  Management Science, 50(7): 949–966.
  * [62] H. Robbins and S. Monro (1951).  A stochastic approximation method.  The annals of mathematical statistics, Sep: 400–407.
  * [63] R. Roet-Green and Y. Yuan (2020).  Information visibility in omnichannel queues.  Available at SSRN, https://ssrn.com/abstract=3485810
  * [64] R. Y. Rubinstein and R. Marcus (1985).  Efficiency of multivariate control variates in Monte Carlo simulation.  Operations Research, 33(3): 661–677.
  * [65] Y. Sakuma, H. Masuyama and E. Fukuda (2020).  A discrete-time single-server Poisson queueing game: Equilibria simulated by an agent-based model.  European Journal of Operational Research, 283(1): 253–264.
  * [66] O. Sebbouh,R. M. Gower and A. Defazio (2021).  Almost sure convergence rates for stochastic gradient descent and stochastic heavy ball.  2021 Conference on Learning Theory, 3935–3971, PMLR.
  * [67] K. Sun, Y. Liu and L. Yang (2020).  Order Ahead for Pickup: Promise or Peril?.  Available at SSRN, https://ssrn.com/abstract=3673617
  * [68] R. Suri and M. A. Zazanis (1988).  Perturbation analysis gives strongly consistent sensitivity estimates for the M/G/1 queue.  Management Science, 34(1): 39–64.
  * [69] R. Suri and Y. T. Leung (1989).  Single run optimization of discrete event simulations—An empirical study using the M/M/l queue.  IIE transactions, 21(1): 35–49.
  * [70] T.A. Taylor (2018).  On-demand service platforms.  Manufacturing & Service Operations Management, 20(4): 704–720.
  * [71] H. Thorisson. (1985).  The queue GI/G/1: Finite moments of the cycle variables and uniform rates of convergence.  Stochastic processes and their applications, 19(1):85–99.
  * [72] J. Wang, S. Cui, and Z. Wang (2019).  Equilibrium strategies in M/M/1 priority queues with balking.  Production and Operations Management, 28(1): 43–62.
  * [73] W. Whitt (1974).  The continuity of queues.  Advances in Applied Probability, 6:175–183.
  * [74] W. Whitt (1989).  Planning queueing simulations.  Management Science, 35(11):1341–1366.
  * [75] J. Xu and B. Hajek (2013).  The supermarket game.  Stochastic Systems, 3(2):405–441.
  * [76] W. Yang and B. L. Nelson (1992).  Multivariate batch means and control variates.  Management Science, 38(10): 1415–1431.
  * [77] D. D. Yao (1995).  S-modular games, with queueing applications.  Queueing Systems, 21(3):449–475.

##  Appendix A Proofs

###  A.1 Proof of Lemma 2

###### Proof.

We show first that every symmetric Nash equilibrium $\mathbf{p}$ satisfies
$\mathbf{p}=\mathbf{f}(\mathbf{p})$. Assume
$\mathbf{p}\in\operatorname*{arg\,max}_{\mathbf{q}\in\Delta}\mathbf{u(p)^{\prime}q}$.
Then any $\mathbf{q}\in\Delta$ satisfies
$\mathbf{u(p)^{\prime}q}\leq\mathbf{u(p)^{\prime}p}$, or equivalently,
$0\leq\mathbf{u(p)^{\prime}(p-q)}$. Assuming $\mathbf{q}\neq\mathbf{p}$, we
have that

$\displaystyle\|\mathbf{p}+\mathbf{u}(\mathbf{p})-\mathbf{q}\|^{2}$
$\displaystyle=\|\mathbf{p}-\mathbf{q}\|^{2}+2\mathbf{u}(\mathbf{p})^{\prime}(\mathbf{p}-\mathbf{q})+\|\mathbf{u}(\mathbf{p})\|^{2}$
$\displaystyle>\|\mathbf{u}(\mathbf{p})\|^{2}=\|\mathbf{p}+\mathbf{u}(\mathbf{p})-\mathbf{p}\|^{2},$

meaning that
$\mathbf{p}=\pi_{\Delta}\big{(}\mathbf{p}+\mathbf{u}(\mathbf{p})\big{)}=\mathbf{f}(\mathbf{p})$.

We proceed to showing that if $\mathbf{p}$ is _not_ a symmetric equilibrium
strategy then $\mathbf{p}\neq\mathbf{f}(\mathbf{p})$. Assume now that
$\mathbf{p}\notin\operatorname*{arg\,max}_{\mathbf{q}\in\Delta}\mathbf{u(p)^{\prime}q}$
so that there exists some strategy $\mathbf{q}\in\Delta$,
$\mathbf{q}\neq\mathbf{p}$, such that $\mathbf{u(p)^{\prime}(q-p)}>0$. Hence,
we can choose $\theta\in(0,1]$ such that

$\theta<\frac{2\mathbf{u(p)^{\prime}(q-p)}}{\|\mathbf{q-p}\|^{2}}.$

Define $\tilde{\mathbf{q}}=(1-\theta)\mathbf{p}+\theta\mathbf{q}$, so that
$\mathbf{p}-\tilde{\mathbf{q}}=\theta(\mathbf{p-q})$. Note that $\Delta$ is a
convex set, thus $\tilde{\mathbf{q}}\in\Delta$. Then,

$\displaystyle\|\mathbf{p}+\mathbf{u(p)}-\tilde{\mathbf{q}}\|^{2}$
$\displaystyle=\|\mathbf{p}-\tilde{\mathbf{q}}\|^{2}+2\mathbf{u(p)}^{\prime}(\mathbf{p}-\tilde{\mathbf{q}})+\|\mathbf{u(p)}\|^{2}$
$\displaystyle=\theta\left(\theta\|\mathbf{q}-\mathbf{p}\|^{2}-2\mathbf{u(p)}^{\prime}(\mathbf{q}-\mathbf{p})\right)+\|\mathbf{u(p)}\|^{2}$
$\displaystyle<\|\mathbf{u(p)}\|^{2}=\|\mathbf{p}+\mathbf{u(p)}-\mathbf{p}\|^{2},$

therefore
$\mathbf{p}\neq\pi_{\Delta}\big{(}\mathbf{p}+\mathbf{u}(\mathbf{p})\big{)}$,
meaning that $\mathbf{p}\neq\mathbf{f}(\mathbf{p})$. ∎

###  A.2 Proof of Theorem 5 and auxiliary results

The proof of Theorem 5 essentially relies on casting the algorithm in (5) into
a specific form, in which the conditions of Theorems 2.4 in [49, Ch. 5.2] can
be verified. Roughly speaking, this theorem relates the asymptotic behavior of
the iterate $\mathbf{p}^{(n)}$ in (5) to solution paths of an appropriate ODE
(Equation (17) below), in which the strategy ${\mathbf{p}}(\cdot)$ is
understood as a continuous function of time. The proof will make use of the
following result:

######  Lemma 11.

Let $\mathbf{p}\in\Delta$ and $\mathbf{x}\in\mathbb{R}^{k}$, then

$\pi_{\Delta}(\mathbf{p+x})=\pi_{\Delta}(\mathbf{p+Hx}),$