Limiting the Oscillations in Queues with Delayed Information Through a Novel Type of Delay Announcement

TL;DR

This paper models queue dynamics with delayed information and velocity signals, showing how optimal velocity information can reduce oscillations and improve system stability, with practical bounds for implementation.

Contribution

It introduces a neutral delay differential equation model incorporating velocity information and derives bounds for optimal velocity information use.

Findings

Velocity information can eliminate or reduce oscillations.

Too much velocity information can induce oscillations.

Derived bounds guide optimal velocity information usage.

Abstract

Many service systems use technology to notify customers about their expected waiting times or queue lengths via delay announcements. However, in many cases, either the information might be delayed or customers might require time to travel to the queue of their choice, thus causing a lag in information. In this paper, we construct a neutral delay differential equation (NDDE) model for the queue length process and explore the use of velocity information in our delay announcement. Our results illustrate that using velocity information can have either a beneficial or detrimental impact on the system. Thus, it is important to understand how much velocity information a manager should use. In some parameter settings, we show that velocity information can eliminate oscillations created by delays in information. We derive a fixed point equation for determining the optimal amount of velocity information that should be used and find closed form upper and lower bounds on its value. When the oscillations cannot be eliminated altogether, we identify the amount of velocity information that minimizes the amplitude of the oscillations. However, we also find that using too much velocity information can create oscillations in the queue lengths that would otherwise be stable.

Figures

Equations

$$\stackrel{\bullet }{q}(t)=\lambda p_i(q_1,\dots ,q_N,\stackrel{\bullet }{q},\dots ,\stackrel{\bullet }{q},\Delta )-\mu q_i(t),\forall i\in \{1,...,N\},$$

$${\text{information about i}}^{th}\text{ queue}=q_i(t-\Delta )+\delta \stackrel{\bullet }{q}(t-\Delta ).$$

$$\displaystyle p_{i}(q_{1},\dots,q_{N},\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{q}}$}_{1},\dots,\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{q}}$}_{N},\Delta)=\frac{\exp\Big{(}-\theta\big{(}q_{i}(t-\Delta)+\delta\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{q}}$}_{i}(t-\Delta)\big{)}\Big{)}}{\sum_{j=1}^{N}\exp\Big{(}-\theta\big{(}q_{j}(t-\Delta)+\delta\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{q}}$}_{j}(t-\Delta)\big{)}\Big{)}},$$

$$\displaystyle p_{1}(q_{1},q_{2})=\frac{\exp(-\theta q_{1})}{\exp(-\theta q_{1})+\exp(-\theta q_{2})}=\frac{1}{1+\exp(-\theta\big{(}q_{2}-q_{1})\big{)}}.$$

$$\displaystyle\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{q}}$}_{i}(t)=\lambda\cdot\frac{\exp\Big{(}-\theta\big{(}q_{i}(t-\Delta)+\delta\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{q}}$}_{i}(t-\Delta)\big{)}\Big{)}}{\sum_{j=1}^{N}\exp\Big{(}-\theta\big{(}q_{j}(t-\Delta)+\delta\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{q}}$}_{i}(t-\Delta)\big{)}\Big{)}}-\mu q_{i}(t),\quad\forall i\in\{1,...,N\}$$

$$q_i(t)=f_i(t),\stackrel{\bullet }{q}(t)=\stackrel{\bullet }{f}(t),t\in [-\Delta ,0].$$

$$\displaystyle\lambda\cdot\frac{\exp\Big{(}-\theta\big{(}q_{i}(t-\Delta)+\delta\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{q}}$}_{i}(t-\Delta)\big{)}\Big{)}}{\sum_{j=1}^{N}\exp\Big{(}-\theta\big{(}q_{j}(t-\Delta)+\delta\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{q}}$}_{i}(t-\Delta)\big{)}\Big{)}}-\mu q_{i}(t),\quad\forall i\in\{1,...,N\}$$

$$q_i(t)=q_i^{\ast }=\frac{\lambda }{N\mu },1\le i\le N.$$

$$\Phi (R,\Delta )$$

$$u_i(t)=q_i(t)-q_i^{\ast }=q_i(t)-\frac{\lambda }{N\mu }.$$

$$\displaystyle\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{u}}$}_{i}(t)\approx-\frac{\lambda\theta}{N}\Big{(}u_{i}(t-\Delta)+\delta u^{\prime}_{i}(t-\Delta)\Big{)}+\frac{\lambda\theta}{N^{2}}\sum_{j=1}^{N}\Big{(}u_{j}(t-\Delta)+\delta u^{\prime}_{j}(t-\Delta)\Big{)}-\mu u_{i}(t)$$

$$\displaystyle\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{{\bf u}}}$}(t)=-\frac{\lambda\theta}{N}\cdot\big{(}{\bf u}(t-\Delta)+\delta{\bf u^{\prime}}(t-\Delta)\big{)}+\frac{\lambda\theta}{N^{2}}A\big{(}{\bf u}(t-\Delta)+\delta{\bf u^{\prime}}(t-\Delta)\big{)}-\mu{\bf u}(t),$$

$$A=VDM,\text{where}V,D,M\in {\mathbb{R}}^{N\times N},$$

$$VM=MV=I,D_{ij}=0\text{if}i\ne j.$$

$$\mathbf{u}(t)=V\mathbf{w}(t).$$

$$\displaystyle V\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{{\bf w}}}$}(t)=-\frac{\lambda\theta}{N}V\big{(}{\bf w}(t-\Delta)+\delta{\bf w^{\prime}}(t-\Delta)\big{)}$$

$$\displaystyle+\frac{\lambda\theta}{N^{2}}VDMV\big{(}{\bf w}(t-\Delta)+\delta{\bf w^{\prime}}(t-\Delta)\big{)}-\mu V{\bf w}(t).$$

$$\displaystyle\raisebox{-0.4pt}{$\stackrel{{\scriptstyle\bullet}}{{{\bf w}}}$}(t)=-\frac{\lambda\theta}{N}\big{(}{\bf w}(t-\Delta)+\delta{\bf w^{\prime}}(t-\Delta)\big{)}$$

$$\displaystyle+\frac{\lambda\theta}{N^{2}}D\big{(}{\bf w}(t-\Delta)+\delta{\bf w^{\prime}}(t-\Delta)\big{)}-\mu{\bf w}(t).$$

$$\stackrel{\bullet }{w}(t)$$

$$\displaystyle-(a+\mu)N=e^{-a\Delta}\lambda\theta\Big{(}(1+a\delta)\cos(b\Delta)+b\delta\sin(b\Delta)\Big{)}$$

$$\displaystyle bN=e^{-a\Delta}\lambda\theta\Big{(}-b\delta\cos(b\Delta)+(1+a\delta)\sin(b\Delta)\Big{)}.$$

$$\cos (b\Delta )$$

$$\sin (b\Delta )$$

$$b=\sqrt{\frac{{\lambda }^2{\theta }^2(a\delta +1{)}^2-e^{2a\Delta }{N}^2(a+\mu {)}^2}{e^{2a\Delta }{N}^2-{\delta }^2{\lambda }^2{\theta }^2}}.$$

$$\frac{\delta \theta \lambda }{N}>e^{a\Delta }>\frac{\theta \lambda (1+a\delta )}{N(a+\mu )}.$$

$$\frac{\delta \theta \lambda }{N}<e^{a\Delta }<\frac{\theta \lambda (1+a\delta )}{N(a+\mu )}.$$

$$\cos (b\Delta )=-\frac{e^{a\Delta }N}{\lambda \theta \delta },\sin (b\Delta )=0,$$

$$\displaystyle-(a+\mu)N=e^{-a\Delta}\lambda\theta\Big{(}(1+a\delta)\cos(b\Delta)+b\delta\sin(b\Delta)\Big{)}$$

$$\displaystyle bN=e^{-a\Delta}\lambda\theta\Big{(}-b\delta\cos(b\Delta)+(1+a\delta)\sin(b\Delta)\Big{)}.$$

$$\cos (b\Delta )$$

$$\sin (b\Delta )$$

$$b=\sqrt{\frac{{\lambda }^2{\theta }^2(a\delta +1{)}^2-e^{2a\Delta }{N}^2(a+\mu {)}^2}{e^{2a\Delta }{N}^2-{\delta }^2{\lambda }^2{\theta }^2}}.$$

Full text

Keywords: neutral delay-differential equation, Hopf bifurcation, perturbations method, operations research, queueing theory, fluid limits, delay announcement, velocity

AMS subject classifications: 34K40, 34K18, 41A10, 37G15, 34K27

1 Introduction

Many corporations and services eagerly adopt new technologies that allow the service managers to interact with their customers. One highly important aspect of the communication is the delay announcement, which informs the customers of their estimated waiting time or queue length. Delay announcements are used in variety of service systems. For example, some hospital emergency rooms display their expected waiting times online. Telephone call centers warn the customers that are placed on hold about extensive waiting times. Amusement parks update the waiting times for different rides, and some transportation networks warn about heavy traffic or delays via road signs.

The popularity of delay announcements among service managers extends beyond customer satisfaction. Informing customers about their waiting times allows managers to influence customer decisions, and the overall system dynamics, which are crucial to a company’s productivity and underlying revenue. As a result, many important questions arise when a service decides to implement a delay announcement for its customers. What type of information, if any, should the service provider give to customers? Are there circumstances when the delay announcement hurts the service provider? How long does it take for the service provider to calculate the delay announcement and disseminate it to customers?

The existing literature explores different ways to give a delay announcement, as well as different response behaviors of the customers. For example, Ibrahim et al. [1] study a service system with applications to telephone call centers. Upon calling and receiving the delay announcement, customers have the option to join the queue, balk (leave immediately), or abandon the queue after spending some time waiting. The authors develop methods for determining the accuracy of the last-to-enter-service (LES) delay announcement, which estimates the waiting time for an incoming customer as the waiting time of the most recent customer who entered service. Under the same options for customer behavior, Jouini et al. [2] consider providing different percentiles of the waiting time distribution as information to their customers. They determine the amount of information that maximizes the number of customers who end up receiving service. Armony and Maglaras [3] model a queue where upon arrival the customers are told the steady state expected waiting time, and in addition are given an option to request a call back. The authors propose a staffing rule that picks the minimum number of service agents that satisfies a set of operational constraints on the performance of the system. Guo and Zipkin [4, 5] allow customers to either join the queue or balk, when the customers are presented either with no information, partial information, or full information about the queues are disclosed. The authors discuss for what situation the extra information is beneficial, and when the addition information can hurt the customers or the service provider.

This paper also explores the impact of the delay announcement on the dynamics of the queueing process. However, the current literature focuses only on services that give the delay announcements to their customers in real-time, while we consider the scenarios when the information itself is delayed. Lags in information are common in services that inform their customers about the waiting times prior to customers’ arrival to the service. Such services are prominent in the context of hospital emergency rooms, highway transportation, amusement park rides, and internet buffer sizing [6, 7, 8, 9]. One specific example is the Citibike bike-sharing network in New York City [10, 11]. Riders can search the availability of bikes on a smartphone app, as shown on Figure (1). However, in the time that it takes for the riders to leave their home and get to a station, all of the bikes could have been taken from that station. Thus, the information they used is delayed and is somewhat unreliable by the time they arrive to the station.

In this paper, we present a deterministic fluid-like model of queues. This may seem as a counter-intuitive choice given that queues usually comprise of discrete units such as the number of people, jobs, or automobiles. However, a queueing system with a heavy traffic flow can be well approximated by a fluid model. These approximations are common in queueing literature [1, 12], where the queues are modeled as stochastic processes and then shown to converge in some limiting regime to deterministic equations. Fluid models are especially useful in settings like manufacturing systems, traffic networks, cloud-computing jobs, busy call centers, and crowds at Disneyland park where the demand for service is large [13, 14, 12, 15, 7].

1.1 Contributions of Paper

We present a fluid model of $N$ queues where customers choose which queue to join, giving preference to the shorter queue based on delayed information. Similar models were previously considered by the authors in [16, 17, 18]. The size of delay in information determines whether the queues approach a stable equilibrium or Hopf bifurcations occur and the queues oscillate indefinitely. The threshold at which the queues become unstable can be affected by the type of information that is revealed to customers. This paper analyzes what kind of information the service managers should provide to customers in order to distribute the workload evenly among the queues. This benefits both the customers who will avoid excessive waits at the longer queues and the servers who will avoid getting overworked or underworked. In many settings, the operator knows not only the current queue lengths, but also the rate at which the queues are changing, namely the queue velocity.

• •

We develop a new queueing model, where customers are told a weighted sum of the queue length and the queue velocity.

• •

We show that the queueing system can undergo a Hopf bifurcation if the delay due to the customer travel time is large. We derive the exact point where the Hopf bifurcation occurs.

• •

We specify how the weight coefficient of velocity information should be chosen so that queues can maintain their stability under greater lags in time. Specifically, we prove that there exists an optimal weight that maximizes the delay where bifurcation occurs.

• •

We derive a fixed point equation for the optimal weight, as well as closed-form expressions for upper and lower bounds on that weight. We also provide upper and lower bounds on the maximum delay where the bifurcation occurs.

• •

When the oscillations in queues cannot be prevented, we use the second-order approximation of amplitude via Lindstedt’s method to determine the weight of velocity information that minimizes the amplitude of oscillations.

• •

When the weights are chosen inadequately, the velocity information can be harmful to the system. We specify the threshold for the weight coefficient where the adverse effects take place.

1.2 Organization of Paper

The remainder of the paper is organized as follows. Section (2) presents a mathematical model for $N$ queues and describes the qualitative behavior of the queueing system. In particular, we prove the existence and uniqueness of the equilibrium and give conditions under which the equilibrium is locally stable. We show that for certain values of parameters, infinitely many Hopf bifurcations may occur.

For some parameters, the queues converge to an equilibrium for sufficiently small delay in information, but as the delay exceeds a certain threshold $\Delta_{cr}$, the equilibrium becomes unstable. Section (3) discusses how the velocity information affects $\Delta_{cr}$, and since the queues are stable only when the delay is less than $\Delta_{cr}$, it becomes our objective to maximize the threshold delay to provide.

Section (4) considers a queueing system with two queues, which is a special case of our $N$-queue model. We prove that all Hopf bifurcations are supercritical. We use a perturbations technique to develop a highly accurate approximation of the amplitude near the bifurcation point, and show that the amplitude of oscillations in queues can be decreased with the right choice of the velocity information weight parameter.

2 The Queueing Model

Customers arrive at a rate $\lambda>0$ to a system of $N$ queues, where they are given information about the waiting times at each queue based on the current queue length and the rate at which the queue is changing. Each customer chooses one of $N$ queues to join, giving probabilistic preference to the shorter queue. Then, customers travel for $\Delta>0$ time units to reach the queue of their choice. The model assumes infinite-server queues with service rate $\mu>0$, which is customary in the operations research literature [19, 20, 21]. This assumption implies that the departure rate for a queue is the service rate $\mu$ multiplied by the total number of customers in that queue. Figure (2) shows the queueing system for $N$ queues. The queue length for each queue $i$ is given by

[TABLE]

where the function $p_{i}$ represents the probability that a customer chooses the $i^{th}$ queue.

Customers are told a weighted sum of the queue length and the queue velocity,

[TABLE]

This information helps customers decide which queue to join. The probability of a customer choosing the $i$-th queue is given by the Multinomial Logit Model (MNL), commonly used to model customer choice in fields of operations research, economics, and applied psychology [22, 23, 24, 25]

[TABLE]

where $\theta>0$ is a standard coefficient of the MNL, $\delta\geq 0$ is the weight of the information about queue’s velocity, and $\Delta>0$ is the delay in time due to customers travelling to the service.

The parameter $\theta$ determines how strong the customer preference is for the shortest queue. For intuition, we will illustrate the MNL model on the simplest model where there are two queues, and the parameters $\delta,\Delta$ are set to [math]. Figure (3) shows the probability of a customer joining the 1st queue, as a function of $\theta$ and the difference in queue lengths $q_{2}-q_{1}$,

[TABLE]

When $\theta\to 0$, customers choose queues arbitrarily giving no preference based on the queue length. Figure (3) indicates the $\theta=0$ by the yellow line, and $p_{1}=p_{2}=0.5$ for any difference in queue lengths. When $\theta\to\infty$, customers always choose the shortest queue, even when the difference in lengths is marginal. This is marked by the black line in Figure (3). However, for simplicity one can set $\theta=1$, which is denoted by the red curve. In this case, when the queues are roughly of equal length, they will be joined with roughly the same probabilities, but ones the difference in the queue lengths increases, the shorter queue will become more preferable.

Complete model.

The incorporation of the probabilities $p_{i}$ into the queueing system provides a system of neutral delay differential equations (NDDE’s) for the queue lengths

[TABLE]

with the initial conditions specified by nonnegative continuous functions $f_{i}$

[TABLE]

2.1 Conditions for Stability and Hopf Bifurcations

In this Section, we describe the behavior of the queues from Equation (2.5). We begin by establishing the existence and uniqueness of solution to the initial value problem (2.5) - (2.6). We note that there exists an extensive analysis of functional differential equations, see for example [26, 27, 28]. The existence and uniqueness of solution for our specific model directly follows from Driver [29], as stated in the result below.

Theorem 2.1**.**

Let $f_{i}(t)$ from Equation (2.6) be absolutely continuous on $t\in[-\Delta,0]$, and \raisebox{-0.4pt}{\stackrel{{\scriptstyle\bullet}}{{f}}}_{i}(t) be bounded for almost all $t\in[-\Delta,0]$ for every $1\leq i\leq N$. Then there exists a solution $q_{i},\dots,q_{N}$ for all $t>0$ that satisfies Equations (2.5) - (2.6). Further, the solution is unique.

Proof.

The existence of the solution is given by Theorem 1 of Driver [29]. The uniqueness of the solution follows from Theorem 2 of Driver [29], but we first need to ensure that the conditions of Theorem 2 are fulfilled. The theorem requires that the function

[TABLE]

satisfies the local Lipschitz condition with respect to $(q_{1}(t),\dots,q_{N}(t))$ with Lipschitz constant $L$ where $L$ is a continuous function of (\raisebox{-0.4pt}{\stackrel{{\scriptstyle\bullet}}{{q}}}_{1}(t-\Delta),\dots,\raisebox{-0.4pt}{\stackrel{{\scriptstyle\bullet}}{{q}}}_{N}(t-\Delta)). Here $q_{i}(t)$ and $q_{i}(t-\Delta)$ are treated as different variables, so the local Lipschitz condition with respect to $(q_{1}(t),\dots,q_{N}(t))$ is satisfied trivially with $L=2\mu$. Therefore, the solution to the system (2.5) - (2.6) is guaranteed to be unique. ∎

Theorem 2.2**.**

The unique equilibrium of $q_{i}(t)$ from Equation (2.5) is given by

[TABLE]

Proof.

See the Appendix for the proof. ∎

The stability of the equilibrium can be determined by the stability of the linearized system of equations [26, 30]. Hence, we proceed by linearizing $q_{i}$ about the equilibrium, and finding the characteristic equation.

Proposition 2.3**.**

The characteristic equation of (2.5) is given by

[TABLE]

Proof.

We introduce the functions $u_{i}(t)$ that represent the deviation of $q_{i}(t)$ from the equilibrium:

[TABLE]

Once the NDDE’s are linearized, \raisebox{-0.4pt}{\stackrel{{\scriptstyle\bullet}}{{u}}}_{i}(t) are approximated as

[TABLE]

In the vector form, we have

[TABLE]

where $A\in\mathbb{R}^{N\times N}$, and $A_{ij}=1$ for $1\leq i,j\leq N$. The matrix $A$ can be diagonalized:

[TABLE]

Since all rows of $A$ are identical, $A$ has only one eigenvalue. This implies diagonal matrix $D$ has only one nonzero element $D_{11}=N$. This property can be exploited with the introduction of a vector ${\bf w}(t)$:

[TABLE]

That is an acceptable form of definition because $V$ is invertible. Equation (2.12) becomes

[TABLE]

Pre-multiplying this equation by $M$ yields the following simplification,

[TABLE]

Writing out $D$ explicitly reduces the system of $N$ equations down to just two equations:

[TABLE]

Equation (2.20) has a solution of the form $w_{1}(t)=ae^{-\mu t}$, so $w_{1}\to 0$ over time. By assuming a solution of the form $w_{i}(t)=e^{Rt}$, the characteristic equation then follows from (2.21). ∎

The equilibrium is stable when all eigenvalues $R$ of the characteristic equation have negative real parts. It is evident that any real root $R$ must be negative. However, there are also infinitely many complex roots, and they depend on the delay $\Delta$. When $\delta>\frac{N}{\lambda\theta}$, the equilibrium cannot be stable for any $\Delta>0$ because there are infinitely many eigenvalues with positive real parts. We demonstrate this in the result below.

Proposition 2.4**.**

Suppose $\delta>\frac{N}{\lambda\theta}$. Then for any $\Delta>0$, there are infinitely many eigenvalues of the characteristic equation that have positive real parts.

Proof.

Suppose $\delta>\frac{N}{\lambda\theta}$. We assume $R=a+ib$ with $a,b\in\mathbb{R}$. We can assume $b\geq 0$ without loss of generality. Plugging in $R$ and separating the real and imaginary parts:

[TABLE]

We find the expressions for sine and cosine to be

[TABLE]

The identity $\sin^{2}(b\Delta)+\cos^{2}(b\Delta)=1$ gives an expression for $b$,

[TABLE]

We will now show that there are infinitely many eigenvalues $R$, where Re$[R]=a>0$, by separately considering the cases when $\delta\mu>1$, $\delta\mu<1$, and $\delta\mu=1$.

Case 1: $\delta\mu>1$.

We will construct an interval $(a_{1},a_{2})$ with $0<a_{1}<a_{2}$, which contains infinitely many values Re$[R]=a$ that together with $b$ from Equation (2.26) satisfy the characteristic equation. We will choose $a$ to be such that both the numerator and the denominator of $b$ are negative, therefore guaranteeing $b$ to be be real. This yields two inequalities

[TABLE]

Since $\frac{\delta\theta\lambda}{N}>1$ by the assumption that $\delta>\frac{N}{\theta\lambda}$, then the inequality $\frac{\delta\theta\lambda}{N}>e^{a\Delta}$ holds for $a\in[0,a_{2})$, where $a_{2}=\frac{1}{\Delta}\ln(\frac{\delta\theta\lambda}{N})>0$. Further, as $a$ increases, the exponent $e^{a\Delta}$ must inevitably outgrow $\frac{\theta\lambda(1+a\delta)}{N(a+\mu)}$, so there exists $a_{1}\geq 0$ such the second part of inequality from Equation (2.27) holds for all $a\geq a_{1}$. Lastly, note that the condition $\frac{\delta\theta\lambda}{N}>\frac{\theta\lambda(1+a\delta)}{N(a+\mu)}$ holds for all $a\geq 0$ because $\delta\mu>1$, so we can choose $a_{1}$ to be less than $a_{2}$, i.e. $a_{1}\in(0,a_{2})$. This shows that there exists an interval $(a_{1},a_{2})$ with $0<a_{1}<a_{2}$ where the inequalities from Equation (2.27) hold, so by Equation (2.26) we have $0\neq b\in R$ for all $a\in(a_{1},a_{2})$.

If $b\in\mathbb{R}$ satisfies Equation (2.25) for some value of $a$, then $R$ is an eigenvalue of the characteristic equation. To show that there are infinitely many eigenvalues with real parts in $(a_{1},a_{2})$, we consider the limit $a\to a_{2}^{-}$, when the denominator of $b$ approaches zero and $b\to\infty$. In this limit, the right hand side of Equation (2.25) will oscillate between $-1$ and $1$ infinite number of times, while the left hand side of Equation (2.25) will converge to [math]. Hence, there are infinitely many solutions to Equation (2.25) with $a\in(a_{1},a_{2})$, and so there are infinitely many eigenvalues with positive real parts.

Case 2: $\delta\mu<1$.

The argument here is analogous to Case 1, except to guarantee that $b$ from Equation (2.26) is real-valued, we will determine an interval in the range of $a$ where the numerator and the denominator of $b$ are positive. We get the condition

[TABLE]

At $a=0$, $\frac{\theta\lambda(1+a\delta)}{N(a+\mu)}=\frac{\theta\lambda\delta}{N\mu\delta}>\frac{\delta\theta\lambda}{N}>1$, so there is an interval $[0,a_{2})$ for $a$ where $e^{a\Delta}<\frac{\theta\lambda(1+a\delta)}{N(a+\mu)}$ holds. Further, $\frac{\delta\theta\lambda}{N}<\frac{\theta\lambda(1+a\delta)}{N(a+\mu)}$ holds for all $a\geq 0$ because $\delta\mu<1$, therefore $a_{1}=\frac{1}{\Delta}\ln(\frac{\delta\theta\lambda}{N})>0$ must be smaller than $a_{2}$. Therefore for all $a\in(a_{1},a_{2})$, with $0<a_{1}<a_{2}$, $b\in\mathbb{R}$.

Just as in Case 1, when $a\to a_{1}^{+}$, $b\to\infty$ so the right hand side of Equation (2.25) will oscillate between 1 and -1 infinitely many times, while the left hand side will converge to [math]. Thus, there will be infinitely many eigenvalues that satisfy the characteristic equation (2.9).

Case 3: $\delta\mu=1$.

In this case, the expressions for sine and cosine simplify to

[TABLE]

so $b=(2k-1)\pi/\Delta$ for $k=1,2,\dots$, and $a=\frac{1}{\Delta}\ln(\frac{\lambda\theta\delta}{N})$. Since $\delta>\frac{N}{\theta\lambda}$, then $a>0$, and the characteristic equation(2.9) has infinitely many eigenvalues with positive real parts. ∎

However, when the weight coefficient $\delta$ is sufficiently small, i.e. $\delta<\frac{N}{\lambda\theta}$, then given a sufficiently small delay the queues converge to a stable equilibrium. As the next result shows, the stability is due to all complex eigenvalues having negative real parts.

Proposition 2.5**.**

Suppose $\delta<\frac{N}{\lambda\theta}$. When $\Delta$ is sufficiently small, all eigenvalues of the characteristic equation have negative real parts.

Proof.

To reach contradiction, let us assume that for any $\Delta_{0}>0$ there exists some $\Delta\in(0,\Delta_{0})$ and an eigenvalue $R=a+ib$ with $a\geq 0$ that satisfy the characteristic equation (2.9). We can assume $b\geq 0$ without loss of generality. Plugging in $R$ and separating the real and imaginary parts:

[TABLE]

Solving for sine and cosine, we find

[TABLE]

The identity $\sin^{2}(b\Delta)+\cos^{2}(b\Delta)=1$ gives an expression for $b$,

[TABLE]

Since $e^{2a\Delta}\geq 1$ and $N>\delta\lambda\theta$ by assumption, the denominator of $b$ is positive, so the numerator of $b$ must be non-negative. Therefore we get inequalities

[TABLE]

From the first inequality, we obtain an upper bound on $a$, $a\leq\frac{\lambda\theta-N\mu}{N-\lambda\theta\delta}$. If $\lambda\theta<N\mu$, then $a<0$ so we reached a contradiction. If $\lambda\theta\geq N\mu$, then we use (2.34) and (2.35) to find an upper bound on $b$:

[TABLE]

However, we note from the cosine equation (2.32) that $\cos(b\Delta)<0$. Since $b$ is non-negative then $b\Delta>\frac{\pi}{2}$ so $b>\frac{\pi}{2\Delta}$ for any $\Delta$. Choose $\Delta_{0}=\frac{\pi}{4B}$. Then for any $\Delta<\Delta_{0}$ we get a contradiction with (2.36):

[TABLE]

Hence when $\Delta_{0}$ is sufficiently small, then for any $\Delta<\Delta_{0}$ the real part of any eigenvalue is negative. ∎

An interesting edge case, however, is when $\delta=\frac{N}{\lambda\theta}$. If the equality holds, three different behaviors may be observed. The equilibrium will be stable, regardless of the size of the delay, if $\delta\mu>1$. However, the equilibrium will be unstable if $\delta\mu<1$. Further, if $\delta\mu=1$ then the behavior of the queues cannot be determined from the characteristic equation (2.9), as all eigenvalues will be purely imaginary. We justify these findings in the result below.

Proposition 2.6**.**

Suppose $\delta=\frac{N}{\lambda\theta}$. If $\delta\mu<1$, then for any $\Delta$ there exists at least one eigenvalue with positive real part. If $\delta\mu>1$, then all eigenvalues have negative real parts. Further, if $\delta\mu=1$ then all eigenvalues are purely imaginary.

Proof.

As in Propositions (2.4) - (2.5), we express the eigenvalue as $R=a+ib$, and then separate the real and imaginary parts of the characteristic equation. The assumption $\delta=\frac{N}{\lambda\theta}$ simplifies the expressions to be

[TABLE]

We will address the three cases separately.

Case 1: $\delta\mu>1$.

To reach contradiction, suppose there exists an eigenvalue with a nonnegative real part, $a\geq 0$. The expression for $b$ is given by

[TABLE]

where the denominator is positive, so the numerator must be nonnegative for $b$ to be real. Therefore $aN+\theta\lambda-Ne^{a\Delta}(a+\mu)>0$. However, the assumption $\delta\mu>1$ is equivalent to $\lambda\theta<N\mu$, so we can show that

[TABLE]

and we reached a contradiction. Thus, if $\delta\mu>1$ then any eigenvalue must have a negative real part.

Case 2: $\delta\mu<1$.

This condition is equivalent to $\lambda\theta<N\mu$. Again, $b$ satisfies Equation (2.40). As $a\to 0^{+}$, $b\to\infty$ so $\sin(b\Delta)$ oscillates between 1 and $-1$ infinitely quickly. Further, as $a\to 0^{+}$ the right hand side of Equation (2.38) goes to zero. Therefore, Equation (2.38) will have infinitely many roots, while Equation (2.39) will be satisfied at each root automatically since $b$ is given by Equation (2.40). Therefore $\delta\mu<1$ implies that the characteristic equation will have infinitely many eigenvalues with positive real parts.

Case 3: $\delta\mu=1$.

This case is equivalent to the condition $\lambda\theta=N\mu$, which simplifies the Equations (2.38) - (2.39) to be

[TABLE]

Hence, $b=(2k+1)\pi/\Delta$, where $k=0,1,2,\dots$, and $1=e^{{a\Delta}}$ so $a=0$. Therefore the roots of the characteristic equation (2.9) are purely imaginary. ∎

Hence, the equilibrium is stable when $\delta=\frac{N}{\lambda\theta}$ and $\delta\mu>1$, or when $\delta<\frac{N}{\lambda\theta}$ and the delay $\Delta$ is sufficiently small. Further, the only way for the equilibrium to become unstable given that $\delta<\frac{N}{\lambda\theta}$, is if a pair of complex eigenvalues crosses from the negative real side of the complex plane into the positive real side. We will determine the threshold value of delay where the stability of the equilibrium may change by finding where the eigenvalues (if any) on the complex plane reach the imaginary axis.

Proposition 2.7**.**

The characteristic equation (2.9) has a pair of purely imaginary solutions $R=\pm i\omega_{cr}$ with $\omega_{cr}$ being real and positive, at each root $\Delta_{cr}$, given that

[TABLE]

and $\Delta_{cr}$ satisfies the transcendental equation

[TABLE]

Proof.

Assume that $R$ from the characteristic equation (2.9) is purely imaginary, $R=\pm i\omega_{cr}$. Plugging in $R$, the real and imaginary parts produce two equations

[TABLE]

We can solve for the values of the sine and cosine functions i.e.

[TABLE]

and by the trigonometric identity $\sin^{2}(\omega_{cr}\Delta_{cr})+\cos^{2}(\omega_{cr}\Delta_{cr})=1$, $\omega_{cr}$ is found. The cosine equation from (2.47) then gives the equation for $\Delta_{cr}$. ∎

Proposition (2.7) provides the infinitely many critical delays $\Delta_{cr}$ as well as the necessary conditions on the other parameters ($\omega_{cr}\in\mathbb{R}$, $\omega_{cr}\neq 0$) for when Hopf bifurcations may occur. This information allows us to prove that a Hopf bifurcation occurs at every $\Delta_{cr}$.

Theorem 2.8**.**

Suppose $\omega_{cr}$ from Equation (2.43) is real and nonzero. Then a Hopf bifurcation occurs at $\Delta=\Delta_{cr}$, where $\Delta_{cr}$ is any positive root of

[TABLE]

Proof.

By Proposition (2.7), at each $\Delta_{cr}$ there is a pair of purely imaginary eigenvalues $R=i\omega_{cr}$, $\bar{R}=-i\omega_{cr}$. A Hopf bifurcation can only occur if $\frac{d}{d\Delta}$Re$[R(\Delta_{cr})]\neq 0$. To verify this, we assume that $R(\Delta)=\alpha(\Delta)+i\omega(\Delta)$. The characteristic equation (2.9) is differentiated with respect to delay, and we find that at $\Delta_{cr}$ where $\alpha=0$ and $\omega=\omega_{cr}$, $\frac{d}{d\Delta}$Re$[R]$ is given by

[TABLE]

The assumption $\omega_{cr}>0$ guarantees the numerator of $\frac{d\alpha}{d\Delta}(\Delta_{cr})$ to be nonzero. To show that the denominator $D$ is nonzero as well, note that it is a quadratic function of $\Delta$, with an absolute minimum at $\Delta^{*}$ such that

[TABLE]

Once $\Delta^{*}$ and $\omega=\omega_{cr}$ from Equations (2.50) and (2.43) are substituted into the denominator $D(\Delta)$ from Equation (2.49), we find that the minimum of $D$ with respect to $\Delta$ is positive:

[TABLE]

Hence the denominator of $\frac{d\alpha}{d\Delta}(\Delta_{cr})$ is positive for any delay $\Delta$, so

[TABLE]

In fact, if $\delta<\frac{N}{\lambda\theta}$ then $\frac{d\alpha}{d\Delta}(\Delta_{cr})>0$ so the eigenvalues always cross from left to right on the complex plane. If $\delta>\frac{N}{\lambda\theta}$ then $\frac{d\alpha}{d\Delta}(\Delta_{cr})<0$ so the eigenvalues always cross from right to left. At each root of $\Delta_{cr}$ there is one purely imaginary pair of eigenvalues, but all other eigenvalues necessarily have a nonzero real part. Hence all roots $\Lambda_{j}\neq R,\bar{R}$ satisfy $\Lambda_{j}\neq mR,m\bar{R}$ for any integer $m$. Therefore all conditions of the infinite-dimensional version of the Hopf Theorem from Hale and Lunel [26] are satisfied, so a Hopf bifurcation occurs at every root $\Delta_{cr}$. ∎

In the proof of Theorem (2.8), for $\delta<\frac{N}{\lambda\theta}$ it is shown that any pair of complex eigenvalues, which crosses the imaginary axis on the complex plane, necessarily crosses from left to right. The implication here is that once the real part of an eigenvalue becomes positive, it remains positive as the delay increases. This allows us to state the conditions for the local stability of the equilibrium.

Theorem 2.9**.**

When $\lambda\theta>N\mu$ and $\delta<\frac{N}{\lambda\theta}$, the equilibrium is locally stable for sufficiently small delay $\Delta$. When either $\lambda\theta\leq N\mu$ and $\delta<\frac{N}{\lambda\theta}$ or $\lambda\theta<N\mu$ and $\delta=\frac{N}{\lambda\theta}$, the equilibrium is locally stable for all $\Delta$.

Proof.

If $\delta=\frac{N}{\lambda\theta}$ and $\lambda\theta<N\mu$ then by Proposition (2.6), for any delay all eigenvalues of the characteristic equation have negative real parts, therefore the equilibrium is locally stable.

If $\delta<\frac{N}{\lambda\theta}$, then by Proposition (2.5) there exists a sufficiently small $\Delta$ such that all eigenvalues of the characteristic equation have negative real parts. The only way for the equilibrium to become unstable is for an eigenvalue to reach the imaginary axis for some $\Delta$. For that to happen, $\omega_{cr}=\sqrt{\frac{\lambda^{2}\theta^{2}-N^{2}\mu^{2}}{N^{2}-\delta^{2}\lambda^{2}\theta^{2}}}\in\mathbb{R}$, $\omega_{cr}\neq 0$ must hold. In case when $\lambda\theta\leq N\mu$ and $\delta<\frac{N}{\lambda\theta}$, then either $\omega_{cr}\notin\mathbb{R}$ or $\omega_{cr}=0$ so the eigenvalues have negative real parts for all $\Delta$. Therefore, again the eigenvalues have negative real parts for all (finite) $\Delta$. Finally, ∎

To summarize, the behavior of the queues from Equation (2.5) can be categorized by two cases, when $\lambda\theta<N\mu$ and $\lambda\theta>N\mu$. In each case, two different types of behavior can be observed, depending on the size of the parameter $\delta$. Hence, there can be four qualitatively different scenarios, as shown in Figure (4). In the following discussion of the two cases, we will refer to this diagram and will explain it in detail.

Case 1: $\lambda\theta<N\mu$.

This case is represented by the regions A and B that are to the left of the vertical line $\lambda\theta=N\mu$ from Figure (4). When $\delta\leq\frac{N}{\lambda\theta}$, or region B, the queues approach a stable equilibrium for any delay $\Delta$. Here all eigenvalues stay on the negative (real) side of the complex plane. As $\Delta$ increases, the complex eigenvalues approach the imaginary axis, but never reach it, as shown in Figure (6). However, when $\delta>\frac{N}{\lambda\theta}$, which is region A of Figure (4), the queues will never be stable, and will undergo infinitely many Hopf bifurcations as the delay increases. For sufficiently small delay $\Delta$, the complex eigenvalues will be on the positive (real) side of the complex plane, and as $\Delta$ will increase, the complex pairs will cross the imaginary axis from right to left, causing Hopf bifurcations to occur as shown in Figure (6). Note, however, that queues will never gain stability because for any delay $\Delta$ there will be eigenvalues with positive real parts.

Case 2: $\lambda\theta>N\mu$.

This case is represented by the regions C and D in Figure (4). When $\delta<\frac{N}{\lambda\theta}$, or region D, the queues will approach a stable equilibrium for a sufficiently small delay $\Delta$. All the eigenvalues will be on the negative (real) side of the complex plane. As the delay $\Delta$ increases, the complex pairs of eigenvalues will move towards the imaginary axis, crossing the axis eventually one by one from left to right as indicated in Figure (8). During the crossing of each pair, a Hopf bifurcation occurs. When $\delta\geq\frac{N}{\lambda\theta}$, which is region C of the diagram (4), the complex eigenvalues cannot reach the imaginary axis, and they all stay to the right side of the imaginary axis on the complex plane as show in Figure (8), so there will never be a stable equilibrium.

Another aspect to point out is the dependence on the MNL parameter $\theta$. When customers join the queues at random, or $\theta\to 0$, the parameters inevitably end up in region B of Figure (4), so the queues will stable for any delay. Alternatively, if customers always join the shortest queue, or $\theta\to\infty$, then for any $\delta>0$ we inevitably end up in region C of Figure (4), so the queues will always be unstable.

3 Achieving Maximum Stability

In physical settings, it is often important to preserve the stability of the queues. Stability evens out the individual waiting times of the customers, minimizing the negative experience. It is therefore useful to know when providing extra information helps to postpone the point of the bifurcation, and when the extra information makes the bifurcation to happen sooner. For example, consider a numerical example from Figures (10) - (12), with two queues and fixed parameters $\lambda,\theta,\mu,$ and $\Delta$. In Figures (10) and (12), $\Delta<\Delta_{cr}$ so the queues converge to an equilibrium over time. However, in Figures (10) and (12) we have $\Delta>\Delta_{cr}$, so the queues oscillate indefinitely. Although the delay $\Delta$ is the same, the change of behavior results from tweaking the parameter $\delta$, which consequently regulates the bifurcation threshold $\Delta_{cr}$.

In this Section, we will consider the scenario $\lambda\theta>N\mu$, where the equilibrium of the queues can become unstable. We will study how the bifurcation threshold $\Delta_{cr}$ changes depending on the weight of the velocity information $\delta$. Our next result shows that the threshold $\Delta_{cr}$ is a concave function of $\delta$.

Proposition 3.1**.**

Suppose $\lambda\theta>N\mu$. Then the function $\Delta_{cr}(\delta)$ is concave for all $\delta\in[0,\frac{N}{\lambda\theta})$.

Proof.

The critical delay $\Delta_{cr}$ is given by Equation (2.48). It is clear that the second derivative $\frac{d^{2}\Delta_{cr}}{d\delta^{2}}$ is negative for all $\delta\in[0,\frac{N}{\lambda\theta})$:

[TABLE]

∎

Proposition (3.1) allows us to show that there exists a specific size of the weight $\delta$ that makes the queueing system optimally stable. We call this size of the weight $\delta_{max}$, and it is such that $\delta=\delta_{max}$ maximizes the threshold $\Delta_{cr}$. In Proposition (3.2), we give an equation that determines $\delta_{max}$, and provide closed-form expressions for an upper and a lower bound of $\delta_{max}$.

Proposition 3.2**.**

Suppose $\lambda\theta>N\mu$. There exists a unique $\delta_{max}\geq 0$ that maximizes a given root $\Delta_{cr}$ for fixed parameters $\lambda,\mu,N,\theta$. It is given by solution of

[TABLE]

Furthermore, $\delta_{max}$ is bounded by $\delta_{1}<\delta_{max}<\delta_{2}$, where

[TABLE]

Proof.

We can treat $\Delta_{cr}$ as a function of $\delta$. The implicit differentiation of (2.44) gives the rate with which $\Delta_{cr}$ changes:

[TABLE]

By Proposition (3.1), $\Delta_{cr}(\delta)$ is concave on the interval $[0,\frac{N}{\lambda\theta})$. Further, it can be shown that $\frac{d}{d\delta}\Delta_{cr}(0)=1>0$ and $\lim_{\delta\to\frac{N}{\lambda\theta}}\frac{d}{d\delta}\Delta_{cr}(\delta)=-\infty<0$, so there is a point $\delta_{max}$ where $\frac{d}{d\delta}\Delta_{cr}(\delta_{max})=0$. Therefore $\Delta_{cr}(\delta)$ reaches its absolute maximum at $\delta_{max}\in(0,\frac{N}{\lambda\theta})$. For intuition, we plot $\frac{d}{d\delta}\Delta_{cr}(\delta)$ in Figure (13).

The value $\delta_{max}$ can be found numerically by solving $\frac{d}{d\delta}\Delta_{cr}(\delta_{max})=0$ from Equation (3.63), alternatively written as (3.58). It is left to find closed-form expressions for the bounds on $\delta_{max}$. By Equation (3.62), we can express $\delta_{max}$ as

[TABLE]

where $\Delta_{0}=\Delta_{cr}(0)<\Delta_{cr}(\delta_{max})$. When solved for $\delta_{max}$, the inequality (3.65) produces an upper bound condition $\delta_{max}<\delta_{2}$ given by Equation (3.60).

To find the lower bound, we note that $\frac{d}{d\delta}\Delta_{cr}(\delta)$ is monotonically decreasing. Thus, $\frac{d}{d\delta}\Delta_{cr}(\delta)<\frac{d}{d\delta}\Delta_{cr}(0)=1$ for all $\delta\in(0,\frac{N}{\lambda\theta})$, and $\Delta_{cr}(\delta)\leq\delta+\Delta_{cr}(0)<\frac{N}{\lambda\theta}+\Delta_{0}$. Therefore by Equation (3.64), we get

[TABLE]

which produces the bound $\delta_{max}>\delta_{1}$ from Equation (3.59) when solved for $\delta_{max}$. ∎

Figures (15) - (15) show $\Delta_{cr}$ as a function of $\lambda$ and $\delta$. For each arrival rate $\lambda$, the maximum $\Delta_{cr}$ is attained for some $\delta$ between the two curves $\delta_{1}$ and $\delta_{2}$. Similarly, Figures (17) and (17) show $\Delta_{cr}$ as a function of $\mu$ and $\delta$ with the two curves $\delta_{1}$ and $\delta_{2}$. As seen in the Figures (15) - (17), the bounds on $\delta_{max}$ are tight.

Besides knowing at which value $\delta$ the maximal bifurcation threshold $\Delta_{cr}$ may occur, it is also important to know how large that threshold actually is. In the next result, we develop bounds for the maximum $\Delta_{cr}$ that can be attained for fixed parameters $\lambda,N,\mu,$ and $\theta$.

Proposition 3.3**.**

The maximum value of a root $\Delta_{cr}$ for fixed parameters $\lambda,\mu$, and $N$, is attained at $\delta_{max}$ and is bounded by $\Delta_{1}<\Delta_{cr}(\delta_{max})<\Delta_{2}$, where

[TABLE]

Proof.

By Proposition (3.2), $\Delta_{cr}(\delta)$ attains its maximum at $\delta=\delta_{max}$. Hence the lower bound $\Delta_{1}<\Delta_{cr}(\delta_{max})$ trivially follows, since $\delta_{max}\neq\delta_{1},\delta_{2}$. To find an upper bound, note that $\frac{d}{d\delta}\Delta_{cr}(\delta)$ is a monotonically decreasing function, so $\frac{d}{d\delta}\Delta_{cr}(\delta_{1})>\frac{d}{d\delta}\Delta_{cr}(\delta)$ for all $\delta>\delta_{1}$, and also that $\frac{d}{d\delta}\Delta_{cr}(\delta_{1})>0$ since $\Delta_{cr}(\delta)$ increases while $\delta<\delta_{max}$. Hence

[TABLE]

In addition, it is known that $\frac{d}{d\delta}\Delta_{cr}(\delta)<0$ when $\delta>\delta_{max}$, so

[TABLE]

Therefore $\Delta_{cr}(\delta_{max})<\min[\Delta_{2a},\Delta_{2b}]=\Delta_{2}$, as desired. ∎

Figure (18) illustrates $\Delta_{cr}(\delta)-\Delta_{cr}(0)$ as a function of $\delta$, with the maximum attained at $\delta_{max}$ and the bounds on the maximum given by $\Delta_{1}$ and $\Delta_{2}$. Further, it is evident from Figure (18) that there is a threshold value, which we call $\delta_{cap}$, that places a cap on the potential utility of the velocity information. When $\delta$ is less than $\delta_{cap}$, the queueing system becomes more stable from the velocity information because $\Delta_{cr}(\delta)>\Delta_{cr}(0)$. However, when $\delta$ exceeds $\delta_{cap}$, the queues become more unstable in the sense that $\Delta_{cr}(\delta)<\Delta_{cr}(0)$. The result below provides an equation for $\delta_{cap}$.

Proposition 3.4**.**

Suppose $\lambda\theta>N\mu$. There exists a unique $\delta_{cap}>0$ such that $\Delta_{cr}(\delta)>\Delta_{cr}(0)$ for all $\delta<\delta_{cap}$, and $\Delta_{cr}(\delta)<\Delta_{cr}(0)$ for all $\delta>\delta_{cap}$. It is given by the solution to

[TABLE]

Proof.

As previously shown, $\Delta_{cr}(\delta)$ is monotonically increasing on $\delta\in[0,\delta_{max})$ and monotonically decreasing on $\delta\in(\delta_{max},\frac{N}{\lambda\theta})$. Further, $\lim_{\delta\to\frac{N}{\lambda\theta}}\Delta_{cr}(\delta)=0<\Delta_{cr}(0)$ since $\Delta_{cr}(0)>0$ by assumption, so there exists exactly one point $\delta_{cap}$ on the interval $(\delta_{max},\frac{N}{\lambda\theta})$ where $\Delta_{cr}(\delta_{cap})=\Delta_{cr}(0)$, and it also follows that $\Delta_{cr}(\delta_{cap})>\Delta_{cr}(0)$ for all $\delta<\delta_{cap}$ and $\Delta_{cr}(\delta_{cap})<\Delta_{cr}(0)$ for all $\delta>\delta_{cap}$. By substituting the expression for $\Delta_{cr}$ from (2.48) into $\Delta_{cr}(0)-\Delta_{cr}(\delta_{cap})=0$ we get Equation (3.73). ∎

To summarize, when $\lambda\theta>N\mu$, the queues are stable when the delay is less than $\Delta_{cr}$. We can therefore provide the most stability for the queues by choosing $\delta$ that maximizes $\Delta_{cr}$, i.e. $\delta_{max}$. Proposition (3.2) proves the existence of $\delta_{max}$, gives an equation describing $\delta_{max}$, and provides closed-form expressions for bounds $\delta_{1}$ and $\delta_{2}$ such that $\delta_{1}<\delta_{max}<\delta_{2}$. Proposition (3.3) also provides bounds $\Delta_{1}$ and $\Delta_{2}$ for the maximum value that $\Delta_{cr}$ can take as a function of $\delta$, so $\Delta_{1}<\Delta_{cr}(\delta_{max})<\Delta_{2}$. Lastly, we show that even if $\delta\neq\delta_{max}$, it is still beneficial to include the velocity information as long as $\delta<\delta_{cap}$. When $\delta$ exceeds $\delta_{cap}$, however, $\Delta_{cr}(\delta)$ becomes less than $\Delta_{cr}(0)$, so the queues are less likely to be stable than if the velocity information was omitted altogether. Proposition (3.4) proves the existence of $\delta_{cap}$ and provides the equation for it.

4 Impact of Velocity Information on the Amplitude

Now that we have a good understanding of how the velocity information impacts the critical delay, we address a more practical question. What is the impact of the velocity on the amplitude of the oscillations? This question is important because it reveals how much the queues will oscillate when they are not in equilibrium. Moreover, it can provide an estimate of how much throughput is lost because of the oscillations (amusement park capacity) or even provide valuable estimates of how much fuel or energy is lost in transportation settings.

Although our previous analysis holds for an arbitrary number of queues, in the sequel, we will demonstrate how $\delta$ affects the amplitude dynamics of the queues in the case of a two queue network. One reason for this restriction is that we must move beyond linearization techniques. In fact, we must use third order Taylor expansions to obtain information about the amplitude, see for example [18]. Thus, many of the matrix techniques we exploited for linearizing the NDDE in Section (2), cannot be used in the context of tensors for the third order Taylor expansion. Thus, for the case of two dimensions, we have the following system of equations

[TABLE]

where as usual $\Delta,\lambda,\mu,\theta>0$ and $\delta\geq 0$. Similar to Section (3), we will consider the scenario with $\lambda\theta>2\mu$, where even for a small $\delta$ the equilibrium of queues loses stability for sufficiently large delay. Our first result shows that the Hopf bifurcations that occur at each root $\Delta_{cr}$ are supercritical.

Theorem 4.1**.**

Suppose $\omega_{cr}\in\mathbb{R}$ and $\omega_{cr}\neq 0$. The NDDE system (4.74) - (4.75) undergoes a supercritical Hopf bifurcation at each root $\Delta_{cr}$. If $\delta\mu<1$ then the limit cycle is born when $\Delta\leq\Delta_{cr}$. If $\delta\mu>1$ then the limit cycle is born when $\Delta\geq\Delta_{cr}$.

Proof.

We will use the method of slow flow, or the Method of Multiple Scales, to determine the stability of the Hopf bifurcations given by Theorem (2.8). This method is often applied to systems of delay differential equations (DDEs) [31, 32, 33]. We note, however, that the stability of the limit cycles can also be determined by showing that the floquet exponent has negative real part, as outlined in Hassard et al. [34].

The first step in the method of slow flow is to consider the perturbation of $q_{1}$ and $q_{2}$ from the equilibrium point $q_{1}^{*}=q_{2}^{*}=\frac{\lambda}{2\mu}$, and to approximate the resulting derivatives by third order Taylor expansion. The two resulting DDEs can be uncoupled when their sum and their difference are taken

[TABLE]

The function $w_{1}(t)=Ce^{-\mu t}$ decays to [math], while the function $w_{2}(t)$ has a Hopf bifurcation at $\Delta_{cr}$ where the periodic solutions are born.

We set $w_{2}(t)=\sqrt{\epsilon}x(t)$ in order to prepare the NDDE for perturbation treatment:

[TABLE]

We replace the independent variable t by two new time variables $\xi=\omega t$ (stretched time) and $\eta=\epsilon t$ (slow time). Then we expand $\Delta$ and $\omega$ about the critical Hopf values:

[TABLE]

The time derivative $\stackrel{{\scriptstyle\bullet}}{{x}}$ becomes

[TABLE]

The expression for $x(t-\Delta)$ may be simplified by Taylor expansion for small $\epsilon$:

[TABLE]

where $x(\xi-\omega_{cr}\Delta_{cr},\eta)=\tilde{x}$. The function $x$ is represented as $x=x_{0}+\epsilon x_{1}+\dots$, and we get

[TABLE]

After these substitutions are made into (4.80), the resulting equation can be separated by the powers of $\epsilon$ into two equations. For the $\epsilon^{0}$ terms, we get an equation for $x_{0}$ without any terms involving $x_{1}$, namely $L(x_{0})=0$, where

[TABLE]

which is satisfied with a solution of the form

[TABLE]

The equation resulting from $\epsilon^{1}$ terms is $L(x_{1})+M(x_{0})=0$. Since $L(x_{1})=0$ is satisfied by a solution of the form (4.88), then the terms from $M(x_{0})$ involving $\cos(\xi)$ and $\sin(\xi)$ are resonant. To eliminate the resonant terms, their coefficients must be [math], which gives two equations for $A(\eta)$ and $B(\eta)$. Switching into polar coordinates, we define $R=\sqrt{A^{2}+B^{2}}$, and find

[TABLE]

In order to find the equilibrium points of $R$ and to discuss their stability, we need to show that the coefficients $c_{1},c_{2}$, and $c_{3}$ are positive. Notice that $c_{3}$ is a quadratic function of $\Delta_{cr}$ with the minimum located at $\Delta_{cr}^{*}$ such that $\frac{d}{d\Delta_{cr}}c_{3}(\Delta_{cr}^{*})=0$, hence

[TABLE]

therefore the denominator of $\frac{dR}{d\eta}$, $c_{3}$, is always positive. Also, we can show $c_{1}$ to be positive. We first note that at the Hopf, Equation (2.47) must be satisfied so $\cos(\omega_{cr}\Delta_{cr})<0$, which implies that $\omega_{cr}\Delta_{cr}>\frac{\pi}{2}$ and

[TABLE]

Next, we note that $c_{1}$ is an increasing linear function of $\Delta_{cr}$, so $c_{1}$ must be positive for any $\Delta_{cr}>\Delta_{cr}^{*}$ where $c_{1}(\Delta_{cr}^{*})=0$. This $\Delta_{cr}^{*}$ is found to be

[TABLE]

Using the inequality in (4.97), we can show by contradiction that $\Delta_{cr}$ is always greater than $\Delta_{cr}^{*}$. Suppose that for some parameters, we have $\Delta_{cr}^{*}>\frac{\pi}{2\omega_{cr}}$. From the equation (2.43), this implies that

[TABLE]

Set $\bar{\delta}=\delta^{2}$. The inequality can be written as

[TABLE]

Notice that the coefficient of $\bar{\delta}^{2}$ is always positive. It can be shown by finding $\mu^{2}$ that minimizes the coefficient, $\mu^{2}=\frac{\lambda^{2}\theta^{2}}{2\pi^{2}}$, and then finding the minimum value of that coefficient, which is \lambda^{4}\theta^{4}\Big{(}1-\frac{1}{\pi^{2}}\Big{)} so it’s clearly positive. This means that $f(\bar{\delta})$ is a convex function, with a minimum at $\bar{\delta}^{*}$

[TABLE]

where the denominator is the same as the coefficient of $\bar{\delta}^{2}$ from Equation (4.104), so it must be positive. The inequalities (4.104) and (4.106) contradict each other, and so $\Delta_{cr}^{*}\leq\frac{\pi}{2\omega_{cr}}$ for all parameters. Hence by Equation (4.97), $\Delta_{cr}>\Delta_{cr}^{*}$, which implies that $c_{1}$ must be positive.

Since $c_{1}$ is positive, the only way for $R$ from (4.89) to have a nonzero equilibrium point is for $c_{2}$ to be also positive. This produces the conditions on the direction of the Hopf

[TABLE]

Recall that $\alpha$ represents the perturbation from $\Delta_{cr}$. So when $\delta\mu<1$, the limit cycle is born when $\Delta$ exceeds $\Delta_{cr}$. If $\delta\mu>1$, then the limit cycle is born when $\Delta$ becomes less than $\Delta_{cr}$. In either case, the equilibrium points of $R(\eta)$ are given by

[TABLE]

Since $c_{1},c_{2},c_{3}>0$, $R_{0}$ is unstable and $R_{1}$ is stable. In its explicit form,

[TABLE]

and it represents the amplitude of the limit cycle near the Hopf. Since $R_{1}$ is stable, then the Hopf bifurcation is supercritical. ∎

Theorem (4.1) establishes that as $\Delta$ increases, the equilibrium becomes unstable and a stable limit cycle is born.

4.1 First-Order Approximation of Amplitude

We would like to choose the weight coefficient $\delta$ in a way that minimizes the amplitude of the oscillation in queues. To do this, we first need to know what the amplitude of the oscillations is as a function of the system parameters. In the following result, we use a perturbation method to approximate the amplitude of oscillations around the bifurcation point.

Proposition 4.2**.**

The amplitude of the oscillations of the queues near the first Hopf can be approximated by $\frac{R_{1}}{2}$, where $R_{1}$ is given by Equation (4.110).

Proof.

The radius of the limit cycle from (4.110) approximates the amplitude of the oscillations of $w_{2}(t)$ from (4.78). By the change of variables given in Equation (4.76), as $t\to\infty$, the behavior of the queues up to a phase shift is

[TABLE]

Thus, the amplitude of oscillations of queues is $\frac{R_{1}}{2}$. ∎

Therefore, when $\Delta$ exceeds $\Delta_{cr}$, the amplitude of oscillations can be approximated to first order by

[TABLE]

when $\Delta-\Delta_{cr}$ is small.

The approximation is accurate when $\delta$ is substantially smaller than the ratio $\frac{2}{\lambda\theta}$. For example, in Figure (20) the queues oscillate throughout time, and the two horizontal lines provide a good approximation of the amplitude of oscillations based on Equation (4.113). However, the approximation becomes inaccurate when $\delta$ approaches $\frac{2}{\lambda\theta}$. As demonstrated in Figure (20), when $\delta=0.195$ and $\frac{2}{\lambda\theta}=0.2$, the approximated amplitude is only about a half of what the actual amplitude is. The discrepancy is observed in Figures (22) - (24) as well. The surface plot in Figure (22) shows the true amplitude based on numerical integration as a function of the delay $\Delta$ and the coefficient $\delta$, while the surface plot in Figure (22) shows the amplitude’s first-order approximation. Furthermore, the surface plot in Figure (24) shows the error of first-order approximation, where the error increases with $\delta$. Finally, Figure (24) provides intuition for why the approximation fails as $\delta$ approaches $\frac{2}{\lambda\theta}$. Figure (24) presents a plot comparing the amplitude and its approximation as functions of delay while $\delta=0.19$ is close to the threshold $\frac{2}{\lambda\theta}=0.2$. The approximation is proportional to $\sqrt{\Delta-\Delta_{cr}}$, while the true amplitude appears to be a linear function of $(\Delta-\Delta_{cr})$ (even though it is not exactly linear).

Since we are interested in using the analytical expression of the amplitude approximation to determine the coefficient $\delta$ that minimizes the amplitude for a given delay, it is important for the approximation to be accurate. As seen from Figure (22), for a fixed delay, say $\Delta=0.5$, the point of the approximated minimum amplitude (at $\delta\approx 0.2$) does not agree with the true minimum amplitude (at $\delta\approx 0.11$). Hence, the first-order approximation of amplitude is insufficient for our purposes, and we must derive the second order of the approximation.

4.2 Second-Order Approximation of Amplitude

The first-order approximation, as seen from Equation (4.113), is of the form

[TABLE]

where $c_{0}$ is a factor determined by the system parameters and is independent of delay. The second-order approximation takes the form

[TABLE]

where $c_{1}$ is also independent of the delay. The full expression for $c_{1}$ is long and messy, so we omit it from this Section. However, the reader can refer to the Appendix (6.2) for the expressions as well as a discussion on how the second-order approximation is obtained. As shown in Figures (26) and (26), the second order approximation performs just as well as the first order approximation when $\delta$ is significantly smaller than $\frac{2}{\lambda\theta}$, but is much more accurate when $\delta$ approaches $\frac{2}{\lambda\theta}$. Figures (28) and (28) confirm that this trend holds throughout the parameter space in $\delta$ and delay $\Delta$. Figure (30) compares the true amplitude with the two approximations when $\delta=0.1$. The next plot in Figure (30) draws the same comparison but when $\delta=0.19$ is closer to its upper limit $\frac{2}{\lambda\theta}=0.2$. It is evident from the two plots that the second-order approximation is significantly more accurate than the first-order approximation, especially as $\delta\to\frac{2}{\lambda\theta}$. Figures (32) and (32) illustrate the same point more systematically, by comparing the errors of first and second order approximations. These surface plots reveal that the higher order approximation decreases the maximum error by a factor of 10.

4.3 Minimizing the Amplitude of Oscillations

Since the second-order approximation is sufficiently accurate, we proceed by using the analytical formula of the second-order approximation to determine the coefficient $\delta$ that minimizes the amplitude of oscillations. Figure (33) shows the numerically computed amplitude, together with its minimum for each delay according to the second-order approximation. The minimum of the amplitude as a function of $\delta$ is found numerically in MATLAB. It is evident that the approximated minimum closely corresponds to where the true minimum is.

Figure (33) shows that the velocity information indeed affects the amplitude of oscillations, and the amplitude can be reduced with a proper choice of the coefficient $\delta$. Figure (33) also reveals an important finding. The value $\delta_{max}$ for the coefficient $\delta$ that maximizes $\Delta_{cr}$ is not the same as $\delta_{amp}$ that minimizes the amplitude of oscillations. Specifically, $\delta_{max}$ is independent of the delay $\Delta$, while $\delta_{amp}$ is a function of the delay. The one point where the two values are guaranteed to be equal each other, $\delta_{max}=\delta_{amp}$, is when $\delta_{amp}$ is computed for the delay eqal to the maximum possible $\Delta_{cr}$, i.e. $\Delta=\Delta_{cr}(\delta_{max})$. Therefore, one should use $\delta_{max}$ as the weight coefficient as long as the delay is less than the bifurcation threshold $\Delta_{cr}$ evaluated at $\delta_{max}$, but when the delay exceeds $\Delta_{cr}$ one should use $\delta_{amp}$ for the weight coefficient instead.

5 Conclusion

This paper answers important questions with regards to businesses incorporating the queue length velocity into the information that is provided to customers via delay announcements. We consider the information passed to the customers about each of $N$ queues to be a linear combination of the current queue length and the rate at which that queue is moving, or the queue velocity,

[TABLE]

with the delay $\Delta$ being the time of customers travelling to the selected queue.

The most evident finding is that the coefficient $\delta$ that weighs the queue velocity information should always be less than the ratio $\frac{N}{\lambda\theta}$. Maintaining this limit guarantees that, at best, the queues will be locally stable for any delay in information. At worst, the queues will be stable when the delay $\Delta$ is sufficiently small, eventually undergoing a Hopf bifurcation at $\Delta=\Delta_{cr}$ and becoming unstable. Alternatively, if $\delta>\frac{N}{\lambda\theta}$, then the queues will never be stable even when the delay in information is infinitesimally small. The reader can refer to Figure (4) for more details.

Even when the condition $\delta<\frac{N}{\lambda\theta}$ is met, significant improvements can still be made by choosing $\delta$ optimally. In the case when queues become unstable as the delay in information increases (so when $\lambda\theta>N\mu$), the weight $\delta$ can shift the delay threshold $\Delta_{cr}$ at which the queues become unstable. In fact, there is exists a ”cap” on the weight $\delta_{cap}$, such that it is safe and beneficial to include the queue velocity whenever $\delta\leq\delta_{cap}$, meaning that the queues will remain stable under greater delay than if the velocity information was omitted. Further, if the threshold $\delta_{cap}$ is exceeded, then queue velocity information will be harmful to the system. In this case, the queues will lose their stability for a smaller delay $\Delta$ than if the queue velocity was omitted altogether. An edge case that exemplifies the usefulness of this discovery is as follows. If we take $\delta\to\frac{N}{\lambda\theta}$, at which point is is clear that $\delta>\delta_{cap}$, then the queues bifurcate almost immediately because $\Delta_{cr}\to 0$ even though the same queues would have remained stable under a much larger delay if $\delta$ was set to 0. Hence, it is important to keep $\delta$ smaller than $\delta_{cap}$.

We also showed that there exists an optimal value for $\delta$ called $\delta_{max}$ that gives the most stability to the queues. For $\delta=\delta_{max}$, queues will be stable for greater delay than is possible given any other choice of $\delta$. We provide an equation from which $\delta_{max}$ can be found numerically, as well as closed form expressions for upper and lower bounds on $\delta_{max}$. Choosing $\delta$ within those bounds is a safe choice for the service managers.

This leads to a natural assessment of the limitations of providing the queue velocity information. The threshold $\Delta_{cr}$ where the queues lose stability can be arbitrarily close to 0 when $\delta$ is chosen poorly, but even the best choice of $\delta$ can only help so much. We provide a formula from which the maximum attainable $\Delta_{cr}$ can be computed. Further, we give expressions on the bounds for that optimal $\Delta_{cr}$ because they don’t rely on $\delta_{max}$ and hence they may be easier to evaluate. This means that while including $\delta$ can always improve the queue dynamics to some degree, there is a limit on how much impact $\delta$ may have.

The presence of the queue velocity information can also affect the amplitude with which the queues oscillate after losing stability. From numerical integration of the queues such as in Figure (22), it is clear that incorporating the queue velocity information can decrease the amplitude of the oscillations, which is beneficial from the managerial perspective. Using a perturbation technique, we derive an analytic expression that approximates the amplitude of oscillations very accurately. Based on the analytic expression, for any delay we can determine the coefficient $\delta$ that will minimize the amplitude of the oscillations. We note that this coefficient as a function of delay, and is not necessarily equal to the coefficient $\delta_{max}$ that maximizes the delay threshold.

In the future, it would be interesting to extend this model to include terms with higher order derivatives. Under the assumption that service managers can measure the information about the queues ($q_{i}^{(2)},q_{i}^{(3)},\dots$), it would be natural to incorporate this data into the information that is provided to the customers:

[TABLE]

The equations describing such a queueing system will no longer be neutral, and may be more complicated. However, such queueing system may answer new questions. One question of significance is to determine the minimum sufficient number of higher-order derivatives ($K$) that should be included in order to guarantee that the queues will be stable for a given delay.

6 Appendix

6.1 Uniqueness and existence of the equilibrium

Proof of Theorem (2.2).

To check that $q_{i}(t)=\frac{\lambda}{N\mu}$ is an equilibrium, plug into Equation (2.5) to get

[TABLE]

To show uniqueness, we will argue by contradiction. Suppose there is another equilibrium given by $\bar{q}_{i}$, $1\leq i\leq N$, and for some $i$ we have $q^{*}_{i}\neq\bar{q}_{i}$. The following condition must hold

[TABLE]

Hence, the mean of $\bar{q}_{i}$ is $\frac{\lambda}{N\mu}$ and since $\bar{q}_{i}$ cannot all be equal to each other, there must exist some $\bar{q}_{s}$ that is smaller than the mean, and some $\bar{q}_{g}$ that is greater than the mean

[TABLE]

This leads to a contradiction:

[TABLE]

[TABLE]

[TABLE]

Since \raisebox{-1.4pt}{\stackrel{{\scriptstyle\bullet}}{{q}}}_{g}(t)\neq 0, then $\bar{q}_{i}(t)$ is not an equilibrium, and the equilibrium (2.8) is unique. ∎

6.2 Approximation to Amplitude of Oscillations in Queues

To see how the velocity information affects the behavior of the queues after a Hopf bifurcation occurs, we need to develop approximations for the amplitude of oscillations. In Section (4), we find a first-order approximation to the amplitude but observe that it is not sufficiently accurate. Hence, we require a second-order approximation. The steps to determing the second-order approximation are outlined below.

This process is very closely related to the staps taken in Theorem (4.1). We begin with Equation (4.80), and expand the time $\tau=\omega t$. Then expand our functions of interest in $\epsilon$ to the second order:

[TABLE]

where $\Delta_{0}$ and $\omega_{0}$ are the delay and frequency at bifurcation, so $\Delta_{0}=\Delta_{cr}$ and $\omega_{cr}$. By collecting all the terms with the like powers of $\epsilon$ into separate equations, we get equations from which we can solve for $x_{0}$ and $x_{1}$. From the equation for $\epsilon^{0}$ we find that $x_{0}(\tau)=A\cos(\tau)$ is a solution. Next, we use the equation for $\epsilon^{1}$ terms to solve for $A$, which has the expression given by Equation (4.110). We can now find $x_{1}$ that has a solution of the form $x_{1}(\tau)=a_{1}\sin(\tau)+a_{2}\cos(\tau)+a_{3}\sin(3\tau)+a_{4}\cos(3\tau)$. The coefficients $a_{3}$ and $a_{4}$ are determined from equation for $\epsilon^{1}$ terms. We impose the initial condition $x^{\prime}(0)=0$ to ensure that the maximum amplitude is at [math], which implies $a_{1}=-3a_{3}$. Lastly, we determine $a_{2}$ by eliminating the secular terms from the equation for $\epsilon^{2}$ terms. Therefore, the second-order approximation of the amplitude of oscillations can be deduced from

[TABLE]

where the coefficients given below:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

To reproduce our numerical results from Section (4.2) - (4.3), set $\epsilon=1$ and $\Delta_{1}=\frac{1}{\epsilon}(\Delta-\Delta_{0})$, with $\Delta_{0}$ given by Equation (2.48). Note that in the equations above there is no presence of $\Delta_{2}$, because we have have set $\Delta_{2}=0$. There is no equation that determines $\Delta_{2}$ and $\Delta_{1}$ uniquely, and the only restriction is that $\Delta=\Delta_{0}+\epsilon\Delta_{1}+\epsilon^{2}\Delta_{2}$. Prior to choosing $\Delta_{2}$ to be 0, we experimented numerically with different combinations of $\Delta_{1}$ and $\Delta_{2}$, and determined that the pair $\Delta_{1}=\frac{1}{\epsilon}(\Delta-\Delta_{0})$ and $\Delta_{2}=0$ results in nearly the most accurate approximation.