# Queues with Delayed Information: Analyzing the Impact of the Choice Model
Function

Philip Doldo
Center for Applied Mathematics
Cornell University
657 Rhodes Hall, Ithaca, NY 14853
<EMAIL_ADDRESS>
Jamol Pender
School of Operations Research and Information Engineering
Center for Applied Mathematics
Cornell University
228 Rhodes Hall, Ithaca, NY 14853
<EMAIL_ADDRESS>

###### Abstract

In this paper, we study queueing systems with delayed information that use a
generalization of the multinomial logit choice model as its arrival process.
Previous literature assumes that the functional form of the multinomial logit
model is exponential. However, in this work we generalize this to different
functional forms. In particular, we compute the critical delay and analyze how
it depends on the choice of the functional form. We highlight how the
functional form of the model can be interpreted as an exponential model where
the exponential rate parameter is uncertain. Furthermore, the rate parameter
distribution is given by the inverse Laplace-Stieltjes transform of the
functional form when it exists. We perform numerous numerical experiments to
confirm our theoretical insights.

##  1 Introduction

The field of queueing theory has been dedicated to understanding the dynamics
of queueing systems and their impact on society. In an age where information
is readily available, it is an important problem to understand how information
affects the decision making of customers who want to join queueing systems. In
some applications, such as amusement parks, transportation systems, and even
telecommunication systems, it has been observed that the information given to
the customers can potentially be delayed Doldo and Pender [5], Novitzky and
Pender [14], Pender et al. [19, 18], Nirenberg et al. [12]. This delay in
information can be caused by factors such as the fact that it can take time to
process and send information about queue lengths to the customers and thus the
information that customers actually receive describes the queue lengths from
some amount of time in the past. Queueing systems where customers are provided
with delayed information have been studied extensively, see for example
Novitzky et al. [16], Doldo and Pender [3, 2], Doldo et al. [6], Novitzky and
Pender [13], Novitzky et al. [15], Atar and Lipshutz [1], Whitt [23].

There are many different types of delayed information that customers can be
provided with and different choices can lead to different choice models and
thus different queueing system dynamics. One of the most common choices used
in the literature is where the customer is provided with the lengths of the
queues from some amount of time in the past. This delay in information is
often modeled with a constant delay, as is the case in Mitzenmacher [11],
Raina and Wischik [20], Kamath et al. [8], Raina et al. [21], Lipshutz and
Williams [10], Lipshutz [9], Atar and Lipshutz [1], Nirenberg et al. [12],
Novitzky et al. [16], Doldo et al. [6], Novitzky et al. [15], Pender et al.
[17, 18], Whitt [23].

Different types of queue length information with a constant delay have been
used, see Doldo and Pender [3] for example, which considers so-called mean-
field queue length information with a constant delay. Other types of
information that have been considered include delayed velocity information and
moving average information, see Dong et al. [7], Novitzky et al. [16], Pender
et al. [17], Novitzky et al. [15]. Some of the literature has examined
queueing systems where the delay is not constant. For example, a queueing
model that uses delayed updating information is examined in Doldo and Pender
[5], Novitzky and Pender [14]. Additionally, queueing models with distributed
delays have been studied by considering distributed delay equation models in
Doldo and Pender [2], Novitzky and Pender [13], though it is worth noting that
distributed delay models exhibit different dynamics in general than simply
replacing the constant delay with a random variable in standard models, as
discussed in Doldo and Pender [4].

In this paper, we will assume a similar model in spirit of Pender et al. [17].
However, our model uses a generalized version of the multinomial logit choice
model which typically assumes an exponential form. The exponential form is
related to a linear utility function, however our functional form will
correspond to a nonlinear utility function. One of our main goals of this work
is to understand how this nonlinearity of the utility function affects the
critical delay which induces a Hopf bifurcation. As discussed in more detail
in Section 3, we will focus our attention on utilities that incorporate the
complementary cumulative distribution function of some probability
distribution. We will examine how the choice of probability distribution used
to induce the choice model can lead to different system dynamics.

###  1.1 Main Contributions of Paper

The contributions of this work can be summarized as follows:

  * •

We consider queueing system choice models informed by more general utilities
that are built from the complementary cumulative distribution functions of
some probability distribution.

  * •

We solve for an equilibrium solution and compute the critical delay at which
the queueing system exhibits a change in stability in terms of the hazard
function corresponding to the given probability distribution.

  * •

We numerically examine how the queueing dynamics change when different
probability distributions induce the choice model, including how stability
properties differ between various common probability distributions.

###  1.2 Organization of Paper

The remainder of this paper is organized as follows. In Section 2 we introduce
the class of queueing models that we will be analyzing. In Section 3, we find
an equilibrium solution and compute the critical delay of our queueing system.
In Section 4 we consider various probability distributions and how their
corresponding choice models impact the dynamics of the queueing system. In
Section 5 we give concluding remarks and discuss possible extensions.

##  2 Generalized Multinomial Logit Queueing Model

In this section we describe the DDE system that we are interested in
analyzing. The previous literature examines delay differential equation (DDE)
fluid models that are based on the following model of an $N$-queue system

$\displaystyle\overset{\bullet}{q}_{i}(t)$ $\displaystyle=$
$\displaystyle\lambda\frac{\exp(-\theta
q_{i}(t-\Delta))}{\sum_{j=1}^{N}\exp(-\theta q_{j}(t-\Delta))}-\mu
q_{i}(t),\hskip 14.22636pti=1,...,N$ (2.1)

where $q_{i}(t)$ represents the length of the $i^{\text{th}}$ queue at time
$t$ for $i=1,...,N$, $\lambda>0$ is the arrival rate, and $\mu>0$ is the
service rate of the infinite-server queue, $\theta>0$ is a choice model
parameter, and $\Delta>0$ is a time delay. The arrival rate $\lambda$ is
multiplied by the probability of joining the $i^{\text{th}}$ queue which is
determined by the choice model which, in this case, is based on a multinomial
logit model. The choice model takes the delayed queue lengths as input to
model providing customers with delayed information about the queue lengths.
More specifically, the logit choice model used in Equation 2.1 assumes that
the utility of an agent joining the $i^{\text{th}}$ queue is

$\displaystyle U_{i}=-\theta q_{i}(t-\Delta)+\epsilon_{i}$ (2.2)

where $\epsilon_{i}$ is distributed according to a standard Gumbel
distribution. See [22] for the details regarding how to construct such a
choice model. We note that the deterministic portion of the utility is
negative so that maximizing the utility corresponds to joining a queue with a
shorter length (neglecting the random component of the utility as some
realizations of the utility could be larger despite having a larger queue
length if the realization of the corresponding random component of the utility
happens to be small enough).

In this paper, we will focus on analyzing the more general DDE system

$\displaystyle\overset{\bullet}{q}_{i}(t)$ $\displaystyle=$
$\displaystyle\lambda\frac{\bar{G}(q_{i}(t-\Delta))}{\sum_{j=1}^{N}\bar{G}(q_{j}(t-\Delta))}-\mu
q_{i}(t),\hskip 14.22636pti=1,...,N$ (2.3)

which uses a more general functional form of the multinomial logit model which
is based on a more general utility in the form

$\displaystyle U_{i}=\log\left(\bar{G}(q_{i}(t-\Delta))\right)+\epsilon_{i}$
(2.4)

where $\bar{G}:\mathbb{R}\to[0,\infty)$ is a decreasing function (so that
longer queue lengths correspond to smaller utilities, neglecting randomness)
and $\epsilon_{i}$ is again distributed according to a standard Gumbel
distribution. This allows the utility to depend nonlinearly on the delayed
queue length. We note that the choice model used in Equation 2.1 corresponds
to when

$\displaystyle\bar{G}(x)=\exp(-\theta x)$ (2.5)

which happens to be the complementary cumulative distribution function
corresponding to an exponential distribution with parameter $\theta$. This
motivated us to focus on the case where $\bar{G}$ is a complementary
cumulative distribution function corresponding to some probability
distribution as this may be a particularly interesting class of decreasing
functions to consider. Under such a choice model, the probability of joining
the $i^{\text{th}}$ queue is

$\displaystyle\frac{\bar{G}(q_{i}(t-\Delta))}{\sum_{j=1}^{N}\bar{G}(q_{j}(t-\Delta))}.$
(2.6)

Additionally, some complementary cumulative distribution functions can be
written as the moment-generating function of some random variable $Y<0$
(negativity of $Y$ is important so that the moment-generating function is a
decreasing function) so that

$\displaystyle\bar{G}(x)=M_{Y}(x):=\mathbb{E}[e^{Yx}]$ (2.7)

where the expectation is taken over the random variable $Y$. An example of
such a case is when the complementary cumulative distribution function
corresponds to a hyperexponential distribution so that

$\displaystyle\bar{G}(x)=\sum_{k=1}^{m}p_{k}e^{-\theta_{k}x}=\mathbb{E}[e^{Yx}]$
(2.8)

where $Y$ is a discrete random variable where $Y=-\theta_{k}<0$ with
probability $p_{k}$ for $k=1,...,m$. In such a case, we note that the
deterministic part of the utility can be viewed as the cumulant-generating
function of the random variable $Y$, that is

$\displaystyle U_{i}=\log(\mathbb{E}[e^{Yq_{i}(t-\Delta)}]),\hskip
14.22636pti=1,...,N.$ (2.9)

We can interpret such a utility as one corresponding to an exponential
distribution where we have some uncertainty about the parameter (modeled by
the random variable $-Y>0$) of the distribution and account for this
uncertainty by averaging over the distribution that we assume the parameter
follows. If we view the utility as a nonlinear transformation of a linear
quantity (namely, $-Yq_{i}(t-\Delta)$), then we can view this as averaging
over the possible values for the slope of this linear quantity. Of course,
such choices of $\bar{G}$ can be interpreted as Laplace transforms. For ease
of notation, let $Z:=-Y>0$ so that

$\displaystyle\bar{G}(x)=M_{Y}(x)$ $\displaystyle=\mathbb{E}[e^{Yx}]$ (2.10)
$\displaystyle=\mathbb{E}[e^{-Zx}]$ (2.11)
$\displaystyle=\int_{0}^{\infty}e^{-zx}f_{Z}(z)dz$ (2.12)
$\displaystyle=\mathcal{L}[f_{Z}](x)$ (2.13)

where $f_{Z}$ is the probability density function of $Z$ and $\mathcal{L}$
denotes the Laplace transform. It follows that we can directly recover the
density $f_{Z}$ from $\bar{G}(x)$ by applying the inverse Laplace transform so
that

$\displaystyle f_{Z}=\mathcal{L}^{-1}[\bar{G}].$ (2.14)

This relation gives us a way of interpreting the choice models corresponding
to various distributions as being induced by an exponential distribution with
uncertainty in its parameter. Ultimately, we will concern ourselves with
better understanding the dynamics of queueing systems with choice models
induced by different probability distributions.

##  3 Queueing System Stability Analysis

In this section, we will perform a stability analysis of the queueing system
with the generalized functional form of the multinomial logit model. We will
find an equilibrium solution of the system 2.3 and use this equilibrium
solution to perform a linearization which will allow us to compute the
critical delay of this system.

It is well-known that when treating the delay $\Delta$ as a bifurcation
parameter, the DDE system 2.1 undergoes a Hopf bifurcation. In particular,
there exists a critical delay value $\Delta_{\text{cr}}$ such that if
$\Delta>\Delta_{\text{cr}}$, then the system is unstable an oscillates
periodically. Computing this critical delay is an important part of
understanding the dynamics of the queueing system. In this work, we will focus
our attention on the choice model and consider how different choice models can
impact the system dynamics. In particular, we will think of a choice model as
being induced by a probability distribution in the sense that the probability
of joining the $i^{\text{th}}$ queue can be written in the form

$\displaystyle\frac{\bar{G}(q_{i}(t-\Delta))}{\sum_{j=1}^{N}\bar{G}(q_{j}(t-\Delta))}$
(3.15)

where $\bar{G}$ is the complementary cumulative distribution function
corresponding to the given probability distribution. As we noted above, the
choice model in system 2.1 is induced by an exponential distribution with
parameter $\theta$. We will examine how the dynamics of the queueing system
change when the choice model is induced by different probability
distributions.

Suppose we are given some random variable $X:\Omega\to\mathbb{R}$ on the
probability space $(\Omega,\mathcal{F},\mathbb{P})$ and we assume that $X$ has
a cumulative distribution function $G(x)=\mathbb{P}(X\leq x)$ and a
corresponding probability density function $g$ exists. We define the
complementary cumulative distribution function to be
$\bar{G}(x)=1-G(x)=\mathbb{P}(X>x)$. We also define the so-called hazard
function corresponding to this distribution to be
$h(x)=\frac{g(x)}{\bar{G}(x)}$. We consider the DDE system

$\displaystyle\overset{\bullet}{q}(t)$ $\displaystyle=$
$\displaystyle\lambda\frac{\bar{G}(q_{i}(t-\Delta))}{\sum_{j=1}^{N}\bar{G}(q_{j}(t-\Delta))}-\mu
q_{i}(t),\hskip 14.22636pti=1,...,N.$ (3.16)

We are interested in better understanding how using complementary cumulative
distribution functions of different probability distributions to build the
choice model can impact the dynamics of the queueing system. We note that the
system 2.3 is a generalization of the system 2.1, the latter of which had its
choice model induced by an exponential distribution with parameter $\theta>0$.
Understanding how the critical delay of the system changes as different choice
models are used is very important to understanding how the system dynamics
change as the critical delay tells us when the system is stable or unstable.
In order to find the critical delay of the system 2.3, we will perform a
linearization analysis for which an equilibrium solution is needed. Thus, we
introduce an equilibrium solution with Theorem 3.1.

######  Theorem 3.1.

The system 2.3 has an equilibrium solution at $q_{1}=\cdots=q_{N}=q^{*}$ where

$q^{*}=\frac{\lambda}{N\mu}.$

###### Proof.

Letting $q_{1}=\cdots=q_{N}=q^{*}$ for some constant $q^{*}\in\mathbb{R}$ to
be determined, we get the system of $N$ identical equations in the form

$\displaystyle 0$
$\displaystyle=\lambda\frac{\bar{G}(q^{*})}{\sum_{j=1}^{N}\bar{G}(q^{*})}-\mu
q^{*}$ (3.17) $\displaystyle=\lambda\frac{\bar{G}(q^{*})}{N\bar{G}(q^{*})}-\mu
q^{*}$ (3.18) $\displaystyle=\frac{\lambda}{N}-\mu q^{*}$ (3.19)

and so it follows that $q^{*}=\frac{\lambda}{N\mu}$ is a solution. ∎

A nice property of system 2.3 is that the equilibrium solution from Theorem
3.1 always exists regardless of the probability distribution that is used to
define the choice model. Now that we have an equilibrium solution, we can
perform the linearization analysis used to get Theorem 3.2. Before proceeding,
we want to note that the hazard function $h(x)$ corresponding to the given
probability distribution is an important quantity that will show up in our
analysis and thus we make a brief comment on hazard functions. We can rewrite
the hazard function corresponding to the distribution of a continuous random
variable $X$ as follows

$\displaystyle h(x)=\frac{g(x)}{\bar{G}(x)}=\frac{g(x)}{\mathbb{P}(X>x)}$
$\displaystyle=\lim_{\delta\to
0}\frac{G(x+\delta)-G(x)}{\delta}\cdot\frac{1}{\mathbb{P}(X>x)}$ (3.20)
$\displaystyle=\lim_{\delta\to 0}\frac{\mathbb{P}(x<X\leq
x+\delta)}{\delta}\cdot\frac{1}{\mathbb{P}(X>x)}$ (3.21)
$\displaystyle=\lim_{\delta\to 0}\frac{1}{\delta}\frac{\mathbb{P}(\\{X\leq
x+\delta\\}\cap\\{X>x\\})}{\mathbb{P}(X>x)}$ (3.22)
$\displaystyle=\lim_{\delta\to 0}\frac{\mathbb{P}(X\leq
x+\delta|X>x)}{\delta}.$ (3.23)

######  Theorem 3.2.

The system 2.3 exhibits a change in stability at the critical value
$\Delta_{\text{cr}}$ of the delay parameter $\Delta$ where

$\displaystyle\Delta_{\text{cr}}$ $\displaystyle=$
$\displaystyle\frac{\text{arccos}\left(\frac{\mu}{C}\right)}{\sqrt{C^{2}-\mu^{2}}}$
(3.24)

where

$\displaystyle
C=-\frac{\lambda}{N}\frac{g\left(\frac{\lambda}{N\mu}\right)}{\bar{G}\left(\frac{\lambda}{N\mu}\right)}=-\frac{\lambda}{N}h\left(\frac{\lambda}{N\mu}\right)$
(3.25)

and $\Delta_{\text{cr}}$ is valid when $C<-|\mu|$.

###### Proof.

We begin by linearizing the system 2.3 about the equilibrium point
$q_{1}=\cdots=q_{N}=q^{*}$ where $q^{*}=\frac{\lambda}{N\mu}$ by introducing
the variables

$q_{i}(t)=q^{*}+u_{i}(t),\hskip 14.22636pti=1,...,N.$

Letting $u(t)=(u_{1}(t),...,u_{N}(t))^{T}$, we obtain the linearized system

$\displaystyle\overset{\bullet}{u}(t)$ $\displaystyle=$ $\displaystyle
Au(t-\Delta)-\mu u(t)$ (3.26)

where

$A=C\begin{bmatrix}\frac{N-1}{N}&-\frac{1}{N}&-\frac{1}{N}\dots-\frac{1}{N}\\\
-\frac{1}{N}&\frac{N-1}{N}&-\frac{1}{N}\dots-\frac{1}{N}\\\
-\frac{1}{N}&-\frac{1}{N}&\frac{N-1}{N}\dots-\frac{1}{N}\\\
\vdots&&\ddots\vdots\\\
-\frac{1}{N}&-\frac{1}{N}&-\frac{1}{N}\dots\frac{N-1}{N}\end{bmatrix}$

where

$C=-\frac{\lambda}{N}\frac{g\left(\frac{\lambda}{N\mu}\right)}{\bar{G}\left(\frac{\lambda}{N\mu}\right)}=-\frac{\lambda}{N}h\left(\frac{\lambda}{N\mu}\right).$

The matrix $A$ has eigenvalues $0$ with multiplicity $1$ and $C$ with
multiplicity $N-1$ and we can diagonalize the system according to the change
of variables $u(t)=Ev(t)$ where $v(t)=(v_{1}(t),...,v_{N}(t))^{T}$ and $E$ is
a matrix whose columns are eigenvectors of $A$ so that $AE=DE$ where
$D=\text{diag}(0,C,...,C)$ is a diagonal matrix with the eigenvalues of $A$ on
its diagonal. Applying this change of variables yields

$\displaystyle\overset{\bullet}{u}(t)$ $\displaystyle=Au(t-\Delta)-\mu u(t)$
(3.27) $\displaystyle\iff$ $\displaystyle E\overset{\bullet}{v}(t)$
$\displaystyle=AEv(t-\Delta)-\mu Ev(t)$ (3.28) $\displaystyle\iff$
$\displaystyle E\overset{\bullet}{v}(t)$ $\displaystyle=EDv(t-\Delta)-\mu
Ev(t)$ (3.29) $\displaystyle\iff$ $\displaystyle\overset{\bullet}{v}(t)$
$\displaystyle=Dv(t-\Delta)-\mu v(t)$ (3.30)

which can be written as

$\displaystyle\overset{\bullet}{v}_{1}(t)$ $\displaystyle=-\mu v_{1}(t)$
(3.31) $\displaystyle\overset{\bullet}{v}_{2}(t)$
$\displaystyle=Cv_{2}(t-\Delta)-\mu v_{2}(t)$ (3.32) $\displaystyle\vdots$
(3.33) $\displaystyle\overset{\bullet}{v}_{N}(t)$
$\displaystyle=Cv_{N}(t-\Delta)-\mu v_{N}(t)$ (3.34)

The first of these equations, Equation 3.31, is an ordinary differential
equation with a solution that decays with time since $\mu>0$. We thus turn our
attention to the remaining $N-1$ equations with are all DDEs in the same form.
By assuming a solution in the form $v_{j}(t)=e^{rt}$ for $j\in{2,...,N}$ and
$r\in\mathbb{C}$, we get the following characteristic equation.

$\displaystyle r-Ce^{-r\Delta}+\mu$ $\displaystyle=$ $\displaystyle 0$ (3.35)

Stability changes when the real part of $r$ changes signs and thus we consider
when $r$ is on the imaginary axis so that $r=i\omega$ for some
$\omega\in\mathbb{R}$. Separating into real and imaginary parts yields

$\displaystyle\mu$ $\displaystyle=C\cos(\omega\Delta)$ (3.36)
$\displaystyle\omega$ $\displaystyle=-C\sin(\omega\Delta)$ (3.37)

so that

$\mu^{2}+\omega^{2}=C^{2}$

and we take

$\omega=\sqrt{C^{2}-\mu^{2}}.$

Looking back at the characteristic equation 3.35, we have that

$\displaystyle r-Ce^{-r\Delta}+\mu$ $\displaystyle=0$ (3.38)
$\displaystyle\iff$ $\displaystyle\frac{C}{\mu+i\omega}$
$\displaystyle=e^{i\omega\Delta}$ (3.39) $\displaystyle\iff$
$\displaystyle\frac{C(\mu-i\omega)}{\mu^{2}+\omega^{2}}$
$\displaystyle=e^{i\omega\Delta}$ (3.40) $\displaystyle\iff$
$\displaystyle\frac{\mu}{C}-i\frac{\omega}{C}$
$\displaystyle=e^{i\omega\Delta}$ (3.41) $\displaystyle\iff$
$\displaystyle\log\left(\frac{\mu}{C}-i\frac{\omega}{C}\right)$
$\displaystyle=i\omega\Delta$ (3.42) $\displaystyle\iff$
$\displaystyle\ln\left(\frac{\mu^{2}+\omega^{2}}{C^{2}}\right)+i\text{arg}\left(\frac{\mu}{C}-i\frac{\omega}{C}\right)$
$\displaystyle=i\omega\Delta$ (3.43) $\displaystyle\iff$ $\displaystyle\Delta$
$\displaystyle=\frac{1}{\omega}\text{arg}\left(\frac{\mu}{C}-i\frac{\omega}{C}\right)$
(3.44)

where we note that $\frac{\mu}{C}-i\frac{\omega}{C}$ is of unit modulus and
thus we can alternatively write

$\displaystyle\Delta_{\text{cr}}=\frac{1}{\omega}\text{arccos}\left(\frac{\mu}{C}\right)$
(3.45)

provided that $\text{arccos}\left(\frac{\mu}{C}\right)\in\mathbb{R}$ because
we require that
$\text{arg}\left(\frac{\mu}{C}-i\frac{\omega}{C}\right)\in\mathbb{R}$ by
definition. This requirement forces $\left|\frac{\mu}{C}\right|\leq 1$ or
equivalently that $C^{2}\geq\mu^{2}$. Additionally, we require that
$\omega\in\mathbb{R}$ so that this inequality becomes strict so that
$C^{2}>\mu^{2}$. One can show that the system is always unstable when
$C>|\mu|$ and we require that $C<-|\mu|$. These conditions ensure that a real
value of $\Delta_{\text{cr}}$ exists. To account for multivaluedness, we take
the smallest value of $\Delta_{\text{cr}}>0$.

For completeness, we can write the critical delay as

$\displaystyle\Delta_{\text{cr}}=\frac{\text{arccos}\left(\frac{\mu}{C}\right)}{\sqrt{C^{2}-\mu^{2}}}=\frac{\text{arccos}\left(-\frac{N\mu\bar{G}\left(\frac{\lambda}{N\mu}\right)}{\lambda
g\left(\frac{\lambda}{N\mu}\right)}\right)}{\sqrt{\left(\frac{\lambda
g\left(\frac{\lambda}{N\mu}\right)}{N\bar{G}\left(\frac{\lambda}{N\mu}\right)}\right)^{2}-\mu^{2}}}.$
(3.46)

∎

As seen in the proof of Theorem 3.2, we can express the critical delay in
terms of the repeated eigenvalue $C$ of the linearized system. In particular,
we see that $C$ directly depends on the hazard function of the given
probability distribution. Using this information, we can explore how using
different probability distributions to create the choice model can give the
system different stability properties.

##  4 Dynamics for Different Distributions

In this section, we numerically examine the dynamics of the queueing system
2.3 for complementary cumulative distribution functions $\bar{G}$
corresponding to various different probability distributions. In analyzing the
dynamics of the queueing system, it is crucial to understand the critical
delay as this provides important information about the qualitative behavior of
the system. As seen in Section 3, the critical delay of the queueing system
depends directly on the eigenvalue $C$ of the linearized system according to
the relation

$\displaystyle\Delta_{\text{cr}}=\frac{\text{arccos}\left(\frac{\mu}{C}\right)}{\sqrt{C^{2}-\mu^{2}}}$
(4.47)

where the eigenvalue $C$ is directly related to the hazard rate of the chosen
probability distribution by

$\displaystyle C=-\frac{\lambda}{N}h\left(\frac{\lambda}{N\mu}\right).$ (4.48)

Therefore, this eigenvalue is a useful quantity to examine to help understand
the qualitative system dynamics. In the remainder of this section, we consider
when the given probability distribution is exponential, normal, log-normal,
Weibull, gamma, or of phase-type. For each of these distributions, we will
provide useful quantities including the probability density function, denoted
$g$, the complementary cumulative distribution function, denoted $\bar{G}$,
the mean, the variance, and the eigenvalue $C$. Additionally, we will
numerically examine how the critical delay depends on the mean and variance of
the given distribution.

###  4.1 Exponential Distribution

In this section, we assume that the complementary cumulative distribution
function $\bar{G}$ that characterizes the choice model is given by an
exponential distribution. Below we provide some useful quantities relating to
the exponential distribution.

$\displaystyle X$ $\displaystyle\sim\text{Exp}(\theta),\hskip
14.22636pt\theta>0$ (4.49) $\displaystyle g(x)$ $\displaystyle=\theta
e^{-\theta x},\hskip 14.22636ptx\geq 0$ (4.50) $\displaystyle\bar{G}(x)$
$\displaystyle=e^{-\theta x}$ (4.51) $\displaystyle\mathbb{E}[X]$
$\displaystyle=\frac{1}{\theta}$ (4.52) $\displaystyle\text{Var}(X)$
$\displaystyle=\frac{1}{\theta^{2}}$ (4.53) $\displaystyle C$
$\displaystyle=-\frac{\lambda\theta}{N}$ (4.54)

The exponential distribution is a well-known nonnegative continuous
probability distribution. The DDE system 2.3 corresponding to an exponential
distribution is equivalent to a DDE queueing system with a multinomial logit
choice model, which has already been extensively studied in the existing
literature. Such a choice model uses $-\theta q_{i}(t-\Delta)$ for the
deterministic part of the utility of joining the $i^{\text{th}}$ queue so that
the utility depends linearly on the delayed queue length. We note that the
exponential distribution is determined by a single parameter $\theta>0$ as we
see that its mean and variance are both fully determined by this parameter.
Additionally, the exponential distribution has a constant hazard function
given by $h(x)=g(x)/\bar{G}(x)=\theta$.

The exponential complementary cumulative distribution function can be viewed
as the moment generating function of a random variable $Y$ that takes on the
value $-\theta$ with probability $1$ and thus, in accordance with Equation
2.14, we have

$\displaystyle f_{-Y}(x)=\mathcal{L}^{-1}[e^{-\theta x}]=\delta(x-\theta)$
(4.55)

where $\delta(x)$ denotes a Dirac delta centered at the origin. Obviously,
this is a trivial case where we have complete certainty of the parameter of
our exponential distribution.

In Figure 1 we show queue length plots and phase diagrams on each side of the
critical delay. In particular, we observe that the queue lengths approach an
equilibrium when the delay is below the critical delay and the queue lengths
oscillate with some limiting amplitude when the delay is larger than the
critical delay. It is worth noting that the mean and variance of the
exponential distribution are determined by a single parameter and thus the
exponential distribution has a fixed variance for a fixed mean. In Figure 2 we
see how the value of the critical delay varies as the mean of the distribution
is varied. We see that there is a critical value of the mean that, when
approached, can make the critical delay arbitrarily large. This is useful
information because instability can lead to inefficiencies in the queueing
system as the resulting queue length oscillations can cause some servers to be
overworked and others to be underworked. Thus, being able to make the critical
delay large would be beneficial for the overall productivity of the queueing
system as large delays could still result in a stable system.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 1: Before and after the change in stability using the choice model
induced by an exponential distribution with $\theta=1$ with constant history
function on $[-\Delta,0]$ with $q_{1}=4.99$ and $q_{2}=5.01$,
$N=2,\lambda=10$, $\mu=1$. The top two plots are queue length versus time with
$\Delta=.3$ (a) and $\Delta=.7$ (b). The bottom two plots are phase plots of
the queue length derivative with respect to time against queue length for
$\Delta=.3$ (c) and $\Delta=.7$ (d). The critical delay is
$\Delta_{\text{cr}}=.3617$.

Figure 2: Left: The critical delay plotted against the mean of the exponential
distribution that induces the choice model. Right: A plot of the probability
density function used for some selected values of the mean.

###  4.2 Normal Distribution

In this section, we assume that the complementary cumulative distribution
function $\bar{G}$ that characterizes the choice model is given by a normal
distribution. Below we provide some useful quantities relating to the normal
distribution.

$\displaystyle X$ $\displaystyle\sim\text{Normal}(\alpha,\sigma^{2})$ (4.56)
$\displaystyle g(x)$
$\displaystyle=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{x-\alpha}{\sigma}\right)^{2}\right)$
(4.57) $\displaystyle\bar{G}(x)$
$\displaystyle=\frac{1}{2}\left[1-\text{erf}\left(\frac{x-\alpha}{\sigma\sqrt{2}}\right)\right]$
(4.58) $\displaystyle\mathbb{E}[X]$ $\displaystyle=\alpha$ (4.59)
$\displaystyle\text{Var}(X)$ $\displaystyle=\sigma^{2}$ (4.60) $\displaystyle
C$
$\displaystyle=-\sqrt{\frac{2}{\pi}}\frac{\lambda}{N\sigma}\frac{\exp\left(-\frac{1}{2}\left(\frac{\frac{\lambda}{N\mu}-\alpha}{\sigma}\right)^{2}\right)}{\left[1-\text{erf}\left(\frac{\frac{\lambda}{N\mu}-\alpha}{\sigma\sqrt{2}}\right)\right]}$
(4.61)

An important difference between the normal distribution and the exponential
distribution is that the normal distribution is characterized by two
parameters instead of just one parameter, its mean parameter $\alpha$ and
standard deviation parameter $\sigma$. This is interesting to consider because
we can vary the variance of the normal distribution for a fixed mean and vice
versa and thus we have the freedom to adjust the value of the critical delay
even if one of these parameters remains fixed.

The normal distribution has an interesting connection to hazard functions. In
particular, if one considers the mean of a so-called truncated normal
distribution, one can show that for $X\sim\text{Normal}(\alpha,\sigma^{2})$

$\displaystyle\mathbb{E}[X|X<a]$ $\displaystyle=\alpha-\sigma
h\left(\frac{a-\alpha}{\sigma}\right)$ (4.62) $\displaystyle\mathbb{E}[X|X>a]$
$\displaystyle=\alpha+\sigma h\left(\frac{a-\alpha}{\sigma}\right)$ (4.63)

where $h(x)$ is the hazard function of the standard normal distribution. This
is true largely because the probability density function $\phi$ of the
standard normal distribution has the nice property that
$\frac{d}{dx}\phi(x)=-x\phi(x)$.

In Figure 3 we show queue length plots and phase diagrams on each side of the
critical delay. In Figure 4 we see how the value of the critical delay varies
as the mean is varied with a fixed variance and in Figure 5 we see how it
changes as the variance is varied with a fixed mean. We see that in the cases
considered, there is a critical value at which the critical delay can be made
arbitrarily large when approached.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 3: Before and after the change in stability using the choice model
induced by a normal distribution with $\alpha=1$ and $\sigma=1$ with constant
history function on $[-\Delta,0]$ with $q_{1}=4.99$ and $q_{2}=5.01$,
$N=2,\lambda=10$, $\mu=1$. The top two plots are queue length versus time with
$\Delta=.05$ (a) and $\Delta=.1$ (b). The bottom two plots are phase plots of
the queue length derivative with respect to time against queue length for
$\Delta=.05$ (c) and $\Delta=.1$ (d). The critical delay is
$\Delta_{\text{cr}}=.0767$.

Figure 4: Left: The critical delay plotted against the mean of the normal
distribution that induces the choice model with a fixed variance of 1. Right:
A plot of the probability density function used for some selected values of
the mean.

Figure 5: Left: The critical delay plotted against the variance of the normal
distribution that induces the choice model with a fixed mean of 1. Right: A
plot of the probability density function used for some selected values of the
variance.

###  4.3 Log-Normal Distribution

In this section, we assume that the complementary cumulative distribution
function $\bar{G}$ that characterizes the choice model is given by a log-
normal distribution. Below we provide some useful quantities relating to the
log-normal distribution.

$\displaystyle X$ $\displaystyle\sim\text{Log-Normal}(\alpha,\sigma^{2})$
(4.64) $\displaystyle g(x)$
$\displaystyle=\frac{1}{x\sigma\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{\ln(x)-\alpha}{\sigma}\right)^{2}\right)$
(4.65) $\displaystyle\bar{G}(x)$
$\displaystyle=\frac{1}{2}\left[1-\text{erf}\left(\frac{\ln(x)-\alpha}{\sigma\sqrt{2}}\right)\right]$
(4.66) $\displaystyle\mathbb{E}[X]$
$\displaystyle=\exp\left(\alpha+\frac{\sigma^{2}}{2}\right)$ (4.67)
$\displaystyle\text{Var}(X)$
$\displaystyle=(\exp(\sigma^{2})-1)\exp(2\alpha+\sigma^{2})$ (4.68)
$\displaystyle C$
$\displaystyle=-\sqrt{\frac{2}{\pi}}\frac{\alpha}{\sigma}\frac{\exp\left(-\frac{1}{2}\left(\frac{\ln\left(\frac{\lambda}{N\mu}\right)-\alpha}{\sigma}\right)^{2}\right)}{\left[1-\text{erf}\left(\frac{\ln\left(\frac{\lambda}{N\mu}\right)-\alpha}{\sigma\sqrt{2}}\right)\right]}$
(4.69)

The log-normal distribution is a continuous probability distribution such that
the natural logarithm applied to a log-normal random variable with parameters
$\alpha$ and $\sigma$ results in a random variable that has a normal
distribution with mean $\alpha$ and standard deviation $\sigma$. In Figure 6
we show queue length plots and phase diagrams on each side of the critical
delay. In Figure 7 we see how the value of the critical delay varies as the
mean is varied with a fixed variance and in Figure 8 we see how it changes as
the variance is varied with a fixed mean. As the mean varies, there is a
critical value of the mean such that the critical delay can be made
arbitrarily large when approached. This is similar to the case shown for the
normal distribution, but it differs in that we see that in the case of the
log-normal distribution, the critical delay gets larger as the mean approaches
zero as well. Additionally, the log-normal distribution appears to have a
larger critical variance than the normal distribution does for unit mean, but
the critical delay is larger in the log-normal case for small variance than it
is in the normal case.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 6: Before and after the change in stability using the choice model
induced by a log-normal distribution with $\alpha=-\frac{1}{2}\log(2)$ and
$\sigma=\sqrt{\log(2)}$ (which results in mean 1 and variance 1) with constant
history function on $[-\Delta,0]$ with $q_{1}=4.99$ and $q_{2}=5.01$,
$N=2,\lambda=10$, $\mu=1$. The left two plots are queue length versus time
with $\Delta=.5$ (Left) and $\Delta=1$ (Right). The right two plots are phase
plots of the queue length derivative with respect to time against queue length
for $\Delta=.5$ (Left) and $\Delta=1$ (Right). The critical delay is
$\Delta_{\text{cr}}=.6148$.

Figure 7: Left: The critical delay plotted against the mean of the log-normal
distribution that induces the choice model with a fixed variance of 1. Right:
A plot of the probability density function used for some selected values of
the mean.

Figure 8: Left: The critical delay plotted against the variance of the log-
normal distribution that induces the choice model with a fixed mean of 1.
Right: A plot of the probability density function used for some selected
values of the variance.

###  4.4 Weibull Distribution

In this section, we assume that the complementary cumulative distribution
function $\bar{G}$ that characterizes the choice model is given by a Weibull
distribution. Below we provide some useful quantities relating to the Weibull
distribution.

$\displaystyle X$ $\displaystyle\sim\text{Weibull}(\alpha,\beta)$ (4.70)
$\displaystyle g(x)$ $\displaystyle=\beta\alpha x^{\alpha-1}e^{-\beta
x^{\alpha}},\hskip 14.22636ptx,\alpha,\beta>0$ (4.71)
$\displaystyle\bar{G}(x)$ $\displaystyle=e^{-\beta x^{\alpha}}$ (4.72)
$\displaystyle\mathbb{E}[X]$
$\displaystyle=\beta^{-1/\alpha}\Gamma(1+1/\alpha)$ (4.73)
$\displaystyle\text{Var}(X)$
$\displaystyle=\frac{1}{\beta^{2\alpha}}\left[\Gamma\left(1+\frac{2}{\alpha}\right)-\left(\Gamma\left(1+\frac{1}{\alpha}\right)\right)^{2}\right]$
(4.74) $\displaystyle C$
$\displaystyle=-\frac{\lambda\beta\alpha\left(\frac{\lambda}{\mu
N}\right)^{\alpha-1}}{N}$ (4.75)

The Weibull distribution is a continuous probability distribution with two
positive parameters $\alpha$ and $\beta$. Some of the above quantities are
defined in terms of the gamma function

$\displaystyle\Gamma(x):=\int_{0}^{\infty}z^{x-1}e^{-z}dz,\hskip
14.22636pt\text{Re}(x)>0$ (4.76)

which has the property that $\Gamma(n)=(n-1)!$ when $n$ is a positive integer
and can be viewed as a smooth interpolation of the factorial function. The
Weibull distribution has a polynomial hazard function

$\displaystyle h(x)=\beta\alpha x^{\alpha-1}$ (4.77)

so that when $\alpha<1$ it is a decreasing function and when $\alpha>1$ it is
an increasing function. We note that when $\alpha=1$, the hazard rate is a
constant $\beta$, just like that of the exponential distribution. Furthermore,
the Weibull distribution reduces to an exponential distribution with parameter
$\beta>0$ when $\alpha=1$.

In Figure 9 we show queue length plots and phase diagrams on each side of the
critical delay. In Figures 10 and 12 we see how the value of the critical
delay varies as the mean is varied with a fixed variance and in Figures 11 and
13 we see how it changes as the variance is varied with a fixed mean. We see
that in the cases considered, there is a critical mean value and a critical
variance value at which the critical delay is unbounded. The shape of the
critical delay against variance plot appears similar to the one in the log-
normal case, noting the initial concave down section for small variance.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 9: Before and after the change in stability using the choice model
induced by a Weibull distribution with $\alpha=2$ and
$\beta=\left(\Gamma(1+\frac{1}{\alpha})\right)^{\alpha}$ with constant history
function on $[-\Delta,0]$ with $q_{1}=4.99$ and $q_{2}=5.01$,
$N=2,\lambda=10$, $\mu=1$. The top two plots are queue length versus time with
$\Delta=.03$ (a) and $\Delta=.06$ (b). The bottom two plots are phase plots of
the queue length derivative with respect to time against queue length for
$\Delta=.03$ (c) and $\Delta=.06$ (d). The critical delay is
$\Delta_{\text{cr}}=.0407$.

Figure 10: Left: The critical delay plotted against the mean of the Weibull
distribution that induces the choice model with fixed $\alpha=2$. Right: A
plot of the probability density function used for some selected values of the
mean.

Figure 11: Left: The critical delay plotted against the variance of the
Weibull distribution that induces the choice model with fixed $\alpha=2$.
Right: A plot of the probability density function used for some selected
values of the variance.

Figure 12: Left: The critical delay plotted against the mean of the Weibull
distribution that induces the choice model with fixed $\alpha=\frac{1}{2}$.
Right: A plot of the probability density function used for some selected
values of the mean.

Figure 13: Left: The critical delay plotted against the variance of the
Weibull distribution that induces the choice model with fixed
$\alpha=\frac{1}{2}$. Right: A plot of the probability density function used
for some selected values of the variance.

###  4.5 Gamma Distribution

In this section, we assume that the complementary cumulative distribution
function $\bar{G}$ that characterizes the choice model is given by a gamma
distribution. Below we provide some useful quantities relating to the gamma
distribution.

$\displaystyle X$ $\displaystyle\sim\text{Gamma}(\alpha,\beta),\hskip
14.22636pt\theta>0$ (4.78) $\displaystyle g(x)$
$\displaystyle=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}$
(4.79) $\displaystyle\bar{G}(x)$
$\displaystyle=1-\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta x)$ (4.80)
$\displaystyle\mathbb{E}[X]$ $\displaystyle=\frac{\alpha}{\beta}$ (4.81)
$\displaystyle\text{Var}(X)$ $\displaystyle=\frac{\alpha}{\beta^{2}}$ (4.82)
$\displaystyle C$
$\displaystyle=-\frac{\lambda\frac{\beta^{\alpha}}{\Gamma(\alpha)}\left(\frac{\lambda}{N\mu}\right)^{\alpha-1}e^{-\beta\frac{\lambda}{N\mu}}}{N\left(1-\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta\frac{\lambda}{N\mu})\right)}$
(4.83)

We consider the gamma distribution which can be viewed as a generalization of
the exponential distribution in some sense. In particular, if the parameter
$\alpha$ is a positive integer, then the gamma distribution reduces to an
Erlang distribution which reduces to an exponential distribution if
$\alpha=1$. Some of the above quantities are defined in terms of the gamma
function defined in Equation 4.76, but we note that some quantities also
depend on the lower incomplete gamma function

$\displaystyle\gamma(x,t):=\int_{0}^{t}z^{x-1}e^{-z}dz,\hskip
14.22636pt\text{Re}(x)>0.$ (4.84)

Alternatively, one could choose to express quantities involving the lower
incomplete gamma function in terms of the upper incomplete gamma function

$\displaystyle\Gamma(x,t):=\int_{t}^{\infty}z^{x-1}e^{-z}dz,\hskip
14.22636pt\text{Re}(x)>0$ (4.85)

as we have the apparent relation

$\displaystyle\Gamma(x)=\gamma(x,t)+\Gamma(x,t)$ (4.86)

so that one could instead write

$\displaystyle\bar{G}(x)=\frac{1}{\Gamma(\alpha)}\Gamma(\alpha,\beta x).$
(4.87)

With this in mind, one could alternatively express the eigenvalue as

$\displaystyle
C=-\frac{\lambda}{N}\frac{\beta^{\alpha}\left(\frac{\lambda}{N\mu}\right)^{\alpha-1}e^{-\beta\frac{\lambda}{N\mu}}}{\Gamma(\alpha,\beta\frac{\lambda}{N\mu})}.$
(4.88)

In Figure 14 we show queue length plots and phase diagrams on each side of the
critical delay. In Figure 15 we see how the value of the critical delay varies
as the mean is varied with a fixed variance and in Figure 16 we see how it
changes as the variance is varied with a fixed mean. For unit variance, we see
that there appear to be two critical mean values at which the critical delay
becomes unbounded. For unit mean, there appears to be a single critical value
of the variance.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 14: Before and after the change in stability using the choice model
induced by a Gamma distribution with $\alpha=\beta=1$ (which reduces to an
exponential distribution) with constant history function on $[-\Delta,0]$ with
$q_{1}=4.99$ and $q_{2}=5.01$, $N=2,\lambda=10$, $\mu=1$. The top two plots
are queue length versus time with $\Delta=.3$ (a) and $\Delta=.7$ (b). The
bottom two plots are phase plots of the queue length derivative with respect
to time against queue length for $\Delta=.3$ (c) and $\Delta=.7$ (d). The
critical delay is $\Delta_{\text{cr}}=.3617$.

Figure 15: Left: The critical delay plotted against the mean of the gamma
distribution that induces the choice model with a fixed variance of 1. Right:
A plot of the probability density function used for some selected values of
the mean.

Figure 16: Left: The critical delay plotted against the variance of the gamma
distribution that induces the choice model with a fixed mean of 1. Right: A
plot of the probability density function used for some selected values of the
variance.

###  4.6 Phase Type Distribution

In this section, we assume that the complementary cumulative distribution
function $\bar{G}$ that characterizes the choice model is given by a phase-
type distribution. Below we provide some useful quantities relating to phase-
type distributions.

$\displaystyle X$ $\displaystyle\sim\text{Ph}(\bm{\alpha},\mathbf{S^{0}},S)$
(4.89) $\displaystyle g(x)$
$\displaystyle=\bm{\alpha}\exp({S}x)\mathbf{S^{0}}$ (4.90)
$\displaystyle\bar{G}(x)$ $\displaystyle=\bm{\alpha}\exp({S}x)\mathbf{1}$
(4.91) $\displaystyle\mathbb{E}[X^{n}]$
$\displaystyle=(-1)^{n}n!\bm{\alpha}{S}^{-n}\mathbf{1}$ (4.92)
$\displaystyle\text{Var}(X)$
$\displaystyle=2\bm{\alpha}{S}^{-2}\mathbf{1}-(\bm{\alpha}{S}^{-1}\mathbf{1})^{2}$
(4.93) $\displaystyle C$
$\displaystyle=-\frac{\lambda\bm{\alpha}\exp\left({S}\frac{\lambda}{N\mu}\right)\mathbf{S^{0}}}{N\bm{\alpha}\exp\left({S}\frac{\lambda}{N\mu}\right)\mathbf{1}}$
(4.94)

Phase-type distributions describe the distribution of time spent to reach an
absorbing state in a continuous-time Markov chain with a finite state space
with one absorbing state and the remaining states being transient. Consider a
continuous-time Markov chain with $p+1$ states, $p\in\mathbb{Z}^{+}$ of which
are transient states and the remaining state is an absorbing state. Let the
state space be $\\{0,1,...,p\\}$ with state $0$ the absorbing state and let
$\bm{\alpha}\in\mathbb{R}^{1\times p}$ be a vector of probabilities where the
$i^{\text{th}}$ entry of $\bm{\alpha}$ corresponds to the probability of the
continuous-time Markov chain starting at the $i^{\text{th}}$ state. The
continuous-time Markov chain has transition-rate matrix

$\displaystyle Q=\begin{bmatrix}0&\bm{0}\\\ \bm{S}^{0}&S\end{bmatrix}$ (4.95)

where $\bm{0}\in\mathbb{R}^{1\times p}$ is a vector with each entry $0$,
$S\in\mathbb{R}^{p\times p}$, and $\bm{S}^{0}=-S\bm{1}$ where
$\bm{1}\in\mathbb{R}^{p\times 1}$ is a vector where each entry is $1$. Special
cases of phase-type distributions include the Erlang distribution, which is
the distribution of a sum of exponential random variables, and the
hyperexponential distribution, which is a mixture of exponential
distributions.

For example, a hyperexponential distribution has probability density function
and complementary cumulative distribution function

$\displaystyle g(x)$
$\displaystyle=\sum_{k=1}^{m}p_{k}\theta_{k}e^{-\theta_{k}x}$ (4.96)
$\displaystyle\bar{G}(x)$ $\displaystyle=\sum_{k=1}^{m}p_{k}e^{-\theta_{k}x}$
(4.97)

with $\sum_{k=1}^{m}p_{k}=1$ and $p_{k}\in[0,1]$ for $k=1,...,m$. It follows
that, under the corresponding choice model, the probability of joining the
$i^{\text{th}}$ queue is

$\displaystyle\frac{\bar{G}(q_{i}(t-\Delta))}{\sum_{j=1}^{N}\bar{G}(q_{j}(t-\Delta))}=\frac{\sum_{k=1}^{m}p_{k}e^{-\theta_{k}q_{i}(t-\Delta)}}{\sum_{j=1}^{N}\sum_{k=1}^{m}p_{k}e^{-\theta_{k}q_{j}(t-\Delta)}}$
(4.98)

for $i=1,...,N$. It turns out that we can view this choice of $\bar{G}$ as the
moment-generating function of a discrete random variable $Y<0$ that takes on
the value $Y=\theta_{k}$ with probability $p_{k}$ for $k=1,...,m$.

In Figure 17 we show queue length plots and phase diagrams on each side of the
critical delay when considering the special case of a hyperexponential
distribution. In Figure 18 we see how the value of the critical delay varies
as the mean is varied with a fixed variance and in Figure 19 we see how it
changes as the variance is varied with a fixed mean, both for the case of a
hyperexponential distribution. For the cases considered, there appears to be a
critical value of the variance as well as a few critical values for the mean.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 17: Before and after the change in stability using the choice model
induced by a phase-type distribution with $\bm{\alpha}=(.3,.7)$ and
$S=\text{diag}(-1.8367,-.8367)$ (which results in a hyperexponential
distribution with mean 1 and variance 1) with constant history function on
$[-\Delta,0]$ with $q_{1}=4.99$ and $q_{2}=5.01$, $N=2,\lambda=10$, $\mu=1$.
The top two plots are queue length versus time with $\Delta=.3$ (a) and
$\Delta=.5$ (b). The bottom two plots are phase plots of the queue length
derivative with respect to time against queue length for $\Delta=.3$ (c) and
$\Delta=.5$ (d). The critical delay is $\Delta_{\text{cr}}=.4443$.

Figure 18: Left: The critical delay plotted against the mean of a
hyperexponential distribution (in this case, composed of two evenly-weighted
exponential distributions) that induces the choice model with a fixed variance
of 1. Right: A plot of the probability density function used for some selected
values of the mean.

Figure 19: Left: The critical delay plotted against the variance of the
hyperexponential distribution (in this case, composed of two evenly-weighted
exponential distributions) that induces the choice model with a fixed mean of
1. Right: A plot of the probability density function used for some selected
values of the variance.

##  5 Conclusion

In this paper, we examine choice models informed by utilities that are
functions of complementary cumulative distribution functions for some
probability distribution and consider an infinite-server fluid queueing system
where customers are informed by delayed queue length information. We determine
the critical delay of this queueing system in terms of the hazard function of
the given probability distribution. We consider how using choice models with
functional forms based on various different probability distributions can
impact the dynamics of the queueing system. In particular, we see that it is
often possible to choose a probability distribution with specific mean and
variance to make the critical delay arbitrarily large. This information is
useful because such probability distributions can result in a queueing system
that is robust to large delays in information. Naturally, there is room for
extending this work by considering other information than standard delayed
queue length information, such as delayed queue velocity information or
updating queue length information. Additionally, other probability
distributions could be considered and for larger ranges of parameters. One
could also consider other classes of decreasing functions to base the
functional form of the choice model on. It could be interesting to focus more
attention on choice models induced by an exponential distribution where the
distribution parameter has uncertainty and could potentially be viewed as a
random variable according to various distributions.

##  6 Acknowledgements

We would like to thank the Center for Applied Mathematics at Cornell
University for sponsoring Philip Doldo’s research. Finally, we acknowledge the
gracious support of the National Science Foundation (NSF) for Jamol Pender’s
Career Award CMMI # 1751975.

## References

  * Atar and Lipshutz [2021] Rami Atar and David Lipshutz.  Heavy traffic limits for join-the-shortest-estimated-queue policy using delayed information.  _Mathematics of Operations Research_ , 46(1):268–300, 2021.
  * Doldo and Pender [2021a] Philip Doldo and Jamol Pender.  Multi-Delay Differential Equations: A Taylor Expansion Approach.  _International Journal of Bifurcation and Chaos to appear_ , 2021a.
  * Doldo and Pender [2021b] Philip Doldo and Jamol Pender.  Mean field queues with delayed information.  _arXiv preprint arXiv:2112.05899_ , 2021b.
  * Doldo and Pender [2021c] Philip Doldo and Jamol Pender.  A note on the interpretation of distributed delay equations.  _arXiv preprint arXiv:2106.11413_ , 2021c.
  * Doldo and Pender [2021d] Philip Doldo and Jamol Pender.  Queues with updating information: Finding the amplitude of oscillations.  _International Journal of Bifurcation and Chaos to appear_ , 2021d.
  * Doldo et al. [2021] Philip Doldo, Jamol Pender, and Richard Rand.  Breaking the symmetry in queues with delayed information.  _International Journal of Bifurcation and Chaos_ , 31(09):2130027, 2021.
  * Dong et al. [2019] Jing Dong, Elad Yom-Tov, and Galit B Yom-Tov.  The impact of delay announcements on hospital network coordination and waiting times.  _Management Science_ , 65(5):1969–1994, 2019\.
  * Kamath et al. [2015] Gopal Krishna Kamath, Krishna Jagannathan, and Gaurav Raina.  Car-following models with delayed feedback: local stability and hopf bifurcation.  In _2015 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , pages 538–545. IEEE, 2015.
  * Lipshutz [2017] David Lipshutz.  Exit time asymptotics for small noise stochastic delay differential equations.  _arXiv preprint arXiv:1710.09771_ , 2017.
  * Lipshutz and Williams [2015] David Lipshutz and Ruth J Williams.  Existence, uniqueness, and stability of slowly oscillating periodic solutions for delay differential equations with nonnegativity constraints.  _SIAM Journal on Mathematical Analysis_ , 47(6):4467–4535, 2015.
  * Mitzenmacher [2000] Michael Mitzenmacher.  How useful is old information?  _IEEE Transactions on Parallel and Distributed Systems_ , 11(1):6–20, 2000.
  * Nirenberg et al. [2018] Samantha Nirenberg, Andrew Daw, and Jamol Pender.  The impact of queue length rounding and delayed app information on disney world queues.  In _Simulation Conference (WSC), 2018 Winter_. IEEE, 2018.
  * Novitzky and Pender [2020a] Sophia Novitzky and Jamol Pender.  Queues with delayed information: a probabilistic perspective.  _Cornell University, Ithaca NY_ , 14853, 2020a.
  * Novitzky and Pender [2020b] Sophia Novitzky and Jamol Pender.  To update or not to update: queues with information updates.  _Cornell University, Ithaca NY_ , 14853, 2020b.
  * Novitzky et al. [2019] Sophia Novitzky, Jamol Pender, Richard H Rand, and Elizabeth Wesson.  Nonlinear dynamics in queueing theory: Determining the size of oscillations in queues with delay.  _SIAM Journal on Applied Dynamical Systems_ , 18(1):279–311, 2019.
  * Novitzky et al. [2020] Sophia Novitzky, Jamol Pender, Richard H Rand, and Elizabeth Wesson.  Limiting the oscillations in queues with delayed information through a novel type of delay announcement.  _Queueing Systems_ , 95(3):281–330, 2020.
  * Pender et al. [2017] Jamol Pender, Richard H Rand, and Elizabeth Wesson.  Queues with choice via delay differential equations.  _International Journal of Bifurcation and Chaos_ , 27(04):1730016, 2017.
  * Pender et al. [2018] Jamol Pender, Richard H Rand, and Elizabeth Wesson.  An analysis of queues with delayed information and time-varying arrival rates.  _Nonlinear Dynamics_ , 91(4):2411–2427, 2018\.
  * Pender et al. [2020] Jamol Pender, Richard Rand, and Elizabeth Wesson.  A stochastic analysis of queues with customer choice and delayed information.  _Mathematics of Operations Research_ , 45(3):1104–1126, 2020.
  * Raina and Wischik [2005] Gaurav Raina and Damon Wischik.  Buffer sizes for large multiplexers: Tcp queueing theory and instability analysis.  In _Next Generation Internet Networks, 2005_ , pages 173–180. IEEE, 2005.
  * Raina et al. [2015] Gaurav Raina, Sreelakshmi Manjunath, Sai Prasad, and Krishnamurthy Giridhar.  Stability and performance analysis of compound tcp with rem and drop-tail queue management.  _IEEE/ACM Transactions on Networking_ , 24(4):1961–1974, 2015.
  * Train [2009] Kenneth E Train.  _Discrete choice methods with simulation_.  Cambridge university press, 2009.
  * Whitt [2021] Ward Whitt.  On the many-server fluid limit for a service system with routing based on delayed information.  _Operations Research Letters_ , 49(3):316–319, 2021.


for the duality gap or the distance to the optimal solutions [29, 28], while
in this manuscript we have focused on the primal optimality gap of (1.1).
Whether the logarithmic factors reported in this manuscript can be removed
seems to remain open.

##  6 Concluding Remarks

This manuscript presents a catalyst scheme applicable to both smooth convex
optimization and convex-concave minimax problems. We establish the optimal
iteration and sample complexities of the catalyst scheme when accelerating a
simple variant of stochastic gradient method for convex optimization. We
further develop a novel variant of the stochastic extragradient method for
solving strongly-monotone variational inequalities with optimal complexities.
By catalyzing the proposed extragradient variant, the catalyst scheme further
obtains optimal iteration and sample complexities for strong-convex-strongly-
concave minimax problems, up to a potential logarithmic factor. It remains
highly rewarding to further investigate whether such a logarithmic factor is
indeed unavoidable for minimax problems or could be removed with refined
methods.

## Acknowledgement

The authors express their sincere appreciation to Tianyi Lin for valuable
discussions that inspire the development of this manuscript.

## References

  * [1] Yuanhao Wang and Jian Li.  Improved algorithms for convex-concave minimax optimization.  Advances in Neural Information Processing Systems, 33:4800–4810, 2020.
  * [2] Yujia Jin, Aaron Sidford, and Kevin Tian.  Sharper rates for separable minimax and finite sum optimization via primal-dual extragradient methods.  In Conference on Learning Theory, pages 4362–4415. PMLR, 2022.
  * [3] Renbo Zhao.  Accelerated stochastic algorithms for convex-concave saddle-point problems.  Mathematics of Operations Research, 47(2):1443–1473, 2022.
  * [4] Xuan Zhang, Necdet Serhat Aybat, and Mert Gürbüzbalaban.  Robust accelerated primal-dual methods for computing saddle points.  arXiv preprint arXiv:2111.12743, 2021.
  * [5] Arkadi Nemirovski.  Prox-method with rate of convergence o (1/t) for variational inequalities with lipschitz continuous monotone operators and smooth convex-concave saddle point problems.  SIAM Journal on Optimization, 15(1):229–251, 2004.
  * [6] Galina M Korpelevich.  The extragradient method for finding saddle points and other problems.  Matecon, 12:747–756, 1976.
  * [7] Anatoli Juditsky, Arkadi Nemirovski, and Claire Tauvel.  Solving variational inequalities with stochastic mirror-prox algorithm.  Stochastic Systems, 1(1):17–58, 2011.
  * [8] Yunmei Chen, Guanghui Lan, and Yuyuan Ouyang.  Accelerated schemes for a class of variational inequalities.  Mathematical Programming, 165:113–149, 2017.
  * [9] Yurii Nesterov.  Dual extrapolation and its applications to solving variational inequalities and related problems.  Mathematical Programming, 109(2-3):319–344, 2007.
  * [10] Yurii Nesterov and Laura Scrimali.  Solving strongly monotone variational and quasi-variational inequalities.  2006\.
  * [11] Arkadi Nemirovski, Anatoli Juditsky, Guanghui Lan, and Alexander Shapiro.  Robust stochastic approximation approach to stochastic programming.  SIAM Journal on optimization, 19(4):1574–1609, 2009.
  * [12] Guanghui Lan.  First-order and stochastic optimization methods for machine learning, volume 1.  Springer, 2020.
  * [13] Farzad Yousefian, Angelia Nedić, and Uday V Shanbhag.  On smoothing, regularization, and averaging in stochastic approximation methods for stochastic variational inequality problems.  Mathematical Programming, 165:391–431, 2017.
  * [14] Alfredo N Iusem, Alejandro Jofré, Roberto Imbuzeiro Oliveira, and Philip Thompson.  Extragradient method with variance reduction for stochastic variational inequalities.  SIAM Journal on Optimization, 27(2):686–724, 2017.
  * [15] Shisheng Cui and Uday V Shanbhag.  On the analysis of variance-reduced and randomized projection variants of single projection schemes for monotone stochastic variational inequality problems.  Set-Valued and Variational Analysis, 29:453–499, 2021.
  * [16] Renato DC Monteiro and Benar Fux Svaiter.  On the complexity of the hybrid proximal extragradient method for the iterates and the ergodic mean.  SIAM Journal on Optimization, 20(6):2755–2787, 2010.
  * [17] Cong D Dang and Guanghui Lan.  On the convergence properties of non-euclidean extragradient methods for variational inequalities with generalized monotone operators.  Computational Optimization and applications, 60:277–310, 2015.
  * [18] Yu Malitsky.  Projected reflected gradient methods for monotone variational inequalities.  SIAM Journal on Optimization, 25(1):502–520, 2015.
  * [19] Aleksandr Beznosikov, Valentin Samokhin, and Alexander Gasnikov.  Distributed saddle-point problems: Lower bounds, optimal algorithms and federated gans.  arXiv preprint arXiv:2010.13112, 2021.
  * [20] Aryan Mokhtari, Asuman Ozdaglar, and Sarath Pattathil.  A unified analysis of extra-gradient and optimistic gradient methods for saddle point problems: Proximal point approach.  In International Conference on Artificial Intelligence and Statistics, pages 1497–1507. PMLR, 2020.
  * [21] Leonid Denisovich Popov.  A modification of the arrow-hurwicz method for search of saddle points.  Mathematical notes of the Academy of Sciences of the USSR, 28:845–848, 1980.
  * [22] Paul Tseng.  On linear convergence of iterative methods for the variational inequality problem.  Journal of Computational and Applied Mathematics, 60(1-2):237–252, 1995.
  * [23] Georgios Kotsalis, Guanghui Lan, and Tianjiao Li.  Simple and optimal methods for stochastic variational inequalities, i: operator extrapolation.  SIAM Journal on Optimization, 32(3):2041–2073, 2022.
  * [24] Yu Nesterov.  Smooth minimization of non-smooth functions.  Mathematical programming, 103:127–152, 2005.
  * [25] Antonin Chambolle and Thomas Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of mathematical imaging and vision, 40:120–145, 2011.
  * [26] Antonin Chambolle and Thomas Pock.  On the ergodic convergence rates of a first-order primal–dual algorithm.  Mathematical Programming, 159(1-2):253–287, 2016.
  * [27] Yunmei Chen, Guanghui Lan, and Yuyuan Ouyang.  Optimal primal-dual methods for a class of saddle point problems.  SIAM Journal on Optimization, 24(4):1779–1814, 2014.
  * [28] Yuyuan Ouyang and Yangyang Xu.  Lower complexity bounds of first-order methods for convex-concave bilinear saddle-point problems.  Mathematical Programming, 185(1-2):1–35, 2021.
  * [29] Junyu Zhang, Mingyi Hong, and Shuzhong Zhang.  On lower iteration complexity bounds for the saddle point problems.  arXiv preprint arXiv:1912.07481, 2019.
  * [30] Adam Ibrahim, Waïss Azizian, Gauthier Gidel, and Ioannis Mitliagkas.  Lower bounds and conditioning of differentiable games.  arXiv preprint arXiv:1906.07300, 2019.
  * [31] Tianyi Lin, Chi Jin, and Michael I Jordan.  Near-optimal algorithms for minimax optimization.  In Conference on Learning Theory, pages 2738–2779. PMLR, 2020.
  * [32] Osman Güler.  New proximal point algorithms for convex minimization.  SIAM Journal on Optimization, 2(4):649–664, 1992.
  * [33] Saverio Salzo and Silvia Villa.  Inexact and accelerated proximal point algorithms.  Journal of Convex analysis, 19(4):1167–1192, 2012.
  * [34] Renato DC Monteiro and Benar Fux Svaiter.  An accelerated hybrid proximal extragradient method for convex optimization and its implications to second-order methods.  SIAM Journal on Optimization, 23(2):1092–1125, 2013.
  * [35] Hongzhou Lin, Julien Mairal, and Zaid Harchaoui.  Catalyst acceleration for first-order convex optimization: from theory to practice.  Journal of Machine Learning Research, 18(1):7854–7907, 2018.
  * [36] Guanghui Lan.  An optimal method for stochastic composite optimization.  Mathematical Programming, 133(1-2):365–397, 2012.

##  Appendix A Supplementary Proofs

###### Proof of Lemma 2.2.

Denoting $\delta_{t}=\nabla\phi(u_{t})-g_{t}$, one clearly has
$\mathbb{E}_{\xi_{t}}\left[\delta_{t}\right]=0,~{}\mathbb{E}_{\xi_{t}}\left[\lVert\delta_{t}\rVert_{*}^{2}\right]\leq\sigma^{2}.$
In addition,

$\displaystyle\phi(u_{t})$
$\displaystyle\overset{(a)}{\leq}\phi(u_{t-1})+\langle\nabla\phi(u_{t-1}),u_{t}-u_{t-1}\rangle+\tfrac{L_{\phi}}{2}\lVert
u_{t}-u_{t-1}\rVert^{2}$ $\displaystyle=\phi(u_{t-1})+\langle
g_{t-1},u_{t}-u_{t-1}\rangle+\tfrac{1}{2\eta_{t}}\lVert
u_{t}-u_{t-1}\rVert^{2}-\left(\tfrac{1}{2\eta_{t}}-\tfrac{L_{\phi}}{2}\right)\lVert
u_{t}-u_{t-1}\rVert^{2}+\langle\delta_{t-1},u_{t}-u_{t-1}\rangle$
$\displaystyle\overset{(b)}{\leq}\phi(u_{t-1})+\langle
g_{t-1},u-u_{t-1}\rangle+\tfrac{1}{2\eta_{t}}\lVert
u-u_{t-1}\rVert^{2}-\tfrac{1}{2\eta_{t}}\lVert
u-u_{t}\rVert^{2}+\tfrac{\lVert\delta_{t-1}\rVert_{*}^{2}}{2(1/\eta_{t}-L_{\phi})}$
$\displaystyle=\phi(u_{t-1})+\langle\nabla\phi(u_{t-1}),u-u_{t-1}\rangle+\tfrac{1}{2\eta_{t}}\lVert
u-u_{t-1}\rVert^{2}-\tfrac{1}{2\eta_{t}}\lVert
u-u_{t}\rVert^{2}+\tfrac{\lVert\delta_{t-1}\rVert_{*}^{2}}{2(1/\eta_{t}-L_{\phi})}-\langle\delta_{t-1},u-u_{t-1}\rangle$
$\displaystyle\overset{(c)}{\leq}\phi(u)+\left(\tfrac{1}{2\eta_{t}}-\tfrac{\mu_{\phi}}{2}\right)\lVert
u-u_{t-1}\rVert^{2}-\tfrac{1}{2\eta_{t}}\lVert
u-u_{t}\rVert^{2}+\tfrac{\lVert\delta_{t-1}\rVert_{*}^{2}}{2(1/\eta_{t}-L_{\phi})}-\langle\delta_{t-1},u-u_{t-1}\rangle,$

where $(a)$ follows from $\phi$ being $L_{\phi}$ smooth, $(b)$ follows from
the Young’s inequality together with $\eta_{t}<1/L_{\phi}$, and $(c)$ follows
from $\phi$ being $\mu_{\phi}$ strongly-convex. Define
$\Delta_{t}\coloneqq\mathbb{E}\left[\phi(u_{t})-\phi(u)\right]$ and
$D_{t}^{X}\coloneqq\mathbb{E}\lVert u-u_{t}\rVert^{2}$. Further taking
expectation of the above inequality and using
$\mathbb{E}\left[\langle\delta_{t-1},u-u_{t-1}\rangle\right]=0$ yield

$\displaystyle\Delta_{t}\leq\left(\tfrac{1}{2\eta_{t}}-\tfrac{\mu_{\phi}}{2}\right)D_{t-1}^{X}-\tfrac{1}{2\eta_{t}}D_{t}^{X}+\tfrac{\sigma^{2}}{2(1/\eta_{t}-L_{\phi})}.$

Now multiply both sides of the above relation by
$\tfrac{2\eta_{t}\mu_{\phi}}{\Lambda_{t}}$, we obtain

$\displaystyle\tfrac{2\eta_{t}\mu_{\phi}}{\Lambda_{t}}\Delta_{t}+\tfrac{\mu_{\phi}}{\Lambda_{t}}D^{X}_{t}$
$\displaystyle\leq\tfrac{(1-\eta_{t}\mu_{\phi})\mu_{\phi}}{\Lambda_{t}}D^{X}_{t-1}+\tfrac{\mu_{\phi}\eta_{t}\sigma^{2}}{\Lambda_{t}(1/\eta_{t}-L_{\phi})}.$

Summing up of the above inequality from $t=1$ to $T$, and making use of the
definition of $\Lambda_{t}$, we obtain

$\displaystyle\textstyle\sum\nolimits_{t=1}^{T}\tfrac{2\eta_{t}\mu_{\phi}}{\Lambda_{t}}\Delta_{t}+\tfrac{\mu_{\phi}}{\Lambda_{T}}D_{T}^{X}\leq\mu_{\phi}D_{0}^{X}+\textstyle\sum\nolimits_{t=1}^{T}\tfrac{\mu_{\phi}\eta_{t}\sigma^{2}}{\Lambda_{t}(1/\eta_{t}-L_{\phi})}.$
(A.1)

In addition, it is clear that

$\displaystyle\textstyle\sum\nolimits_{t=1}^{T}\tfrac{2\eta_{t}\mu_{\phi}}{\Lambda_{t}}=\textstyle\sum\nolimits_{t=1}^{T}2\left(\tfrac{1}{\Lambda_{t}}-\tfrac{1}{\Lambda_{t-1}}\right)=\tfrac{2(1-\Lambda_{T})}{\Lambda_{T}}.$

Combining the above observation with (A.1) and the definition of
$\overline{u}_{T}$, we obtain

$\displaystyle\tfrac{2(1-\Lambda_{T})}{\Lambda_{T}}\mathbb{E}\left[\phi(\overline{u}_{T})-\phi(s^{*})\right]+\tfrac{\mu_{\phi}}{\Lambda_{T}}D_{T}^{X}\leq\mu_{\phi}D_{0}^{X}+\textstyle\sum\nolimits_{t=1}^{T}\tfrac{\mu_{\phi}\eta_{t}\sigma^{2}}{\Lambda_{t}(1/\eta_{t}-L_{\phi})}.$

The desired claim follows immediately after simple rearrangements. ∎

###### Proof of Proposition 2.4.

Since $\eta_{t}\leq\tfrac{1}{2L_{\phi}}$, Lemma 2.2 applies and one can
immediately obtain (2.21). Moreover, direct computation yields

$\displaystyle\Lambda_{t}=\Lambda_{t-1}\left(1-\eta_{t}\mu_{\phi}\right)=\tfrac{(t_{0}+1)(t_{0}+2)}{(t+t_{0}+1)(t+t_{0}+2)}.$
(A.2)

and

$\displaystyle\tfrac{\eta_{t}^{2}}{\Lambda_{t}(1-L_{\phi}\eta_{t})}\overset{(a)}{\leq}\tfrac{2\eta_{t}^{2}}{\Lambda_{t}}\overset{(b)}{\leq}\tfrac{32}{\mu^{2}t_{0}^{2}},$
(A.3)

where $(a)$ follows from $\eta_{t}\leq\tfrac{1}{2L_{\phi}}$, and $(b)$ follows
from (A.2) and the definition of $\eta_{t}$. From (A.2) and the definition of
$t_{0}$, it is also clear that when $T\geq t_{0}$, $\Lambda_{T}\leq 1/2$.
Combining this observation with (A.2) and (A.3), we obtain

$\displaystyle\varepsilon$
$\displaystyle=\tfrac{\Lambda_{T}}{1-\Lambda_{T}}\leq 1;$
$\displaystyle\delta$
$\displaystyle=\tfrac{\Lambda_{T}}{2(1-\Lambda_{T})}\textstyle\sum\nolimits_{t=1}^{T}\tfrac{\mu_{\phi}\eta_{t}^{2}\sigma^{2}}{\Lambda_{t}(1-L_{\phi}\eta_{t})}\leq\tfrac{t_{0}^{2}}{T^{2}}\cdot\tfrac{32\sigma^{2}T}{\mu_{\phi}t_{0}^{2}}\leq\tfrac{32\sigma^{2}}{\mu_{\phi}T}.$

The above relations immediately imply (2.23) and thus conclude the proof. ∎