Benjamin Legros
email: benjamin.legros@centraliens.net

Keywords: queueing systems, performance evaluation, M/G/1, threshold policy, strategic behavior Mathematics Subject Classification (2000) 90B22, 60K25, 68M20

published or not. The documents may come

M/G/1 queue with event-dependent arrival rates

realized event may however bias customers' decisions by inducing illusory correlation. For instance, if the last event was a service then a reduction of the queue size is observed. This event is either irrelevant or at least insufficient to evaluate the quality of the queueing system. Yet, one can conclude that this queue has a high speed of service since a customer has just left service. Another bias is to believe that a positive event is usually followed by negative ones. In this case, a customer may believe that since a service completion has just occurred, the next service times will be long.

Although the last observed event is not a rational indicator of the quality of a service system, [START_REF] Aksin | The effect of random waits on customer queue joining and reneging behavior: A laboratory experiment[END_REF] showed that the evolution of the queue size (increasing or decreasing) due to service completions or arrivals impacts strongly the decision of new arriving customers, by using laboratory experiments in which participants experience several observable queues with different characteristics in terms of queue length and service times.

Queues with workload-dependent arrival or service rates have already been widely studied. Therefore, we aim to investigate a queue with another feature for the customers' behavior which is the event-dependency.

In order to analyze the impact of this new feature, we neglect the other aspects of customer's behavior like the workload sensitivity. To the best of our knowledge, this paper is the first to model this behavior. In what follows we describe the queueing model studied in this article.

Model description. We consider a single-server first-come-first-served queueing model with a general service time distribution. We assume that the iid service times have a density, denoted by g. We denote by x the expected service time and by cv the coefficient of variation of the service time distribution; it is the ratio of the standard deviation divided by the expected value. Based on the last event which can either be an arrival or a service completion, the feature of the event-dependency is captured through two exponential arrival rates; λ + and λ -. More precisely, at a given time, if the last event was an arrival (respectively a service), the next customer will arrive in the system after a random time exponentially distributed with parameter λ + (respectively λ -).

Contributions.

The key contributions of this paper are the following.

-Performance Analysis. First, we consider the embedded Markov chain right after a service completion. This allows us to derive the probability generating function explicitly together with the stability condition and the probability of an empty system at departure instants. Contrary to the M/G/1 queue, the queue length distribution is not identical at arrival and arbitrary instants. Therefore, we also study the relation between these distributions. This leads to explicit expressions of the performance measures at arbitrary instants. Next, we consider three special cases of service time distributions (Exponential, Erlang and Deterministic). We show that our queueing model performs better than a system in which customers are insensitive to the last event if a service generates a higher arrival rate than an arrival (λ -> λ + ). Contrary to the M/G/1 queue, we show that the performance of our model can improve when the coefficient of variation of the service increases.

-Equilibrium Strategy. We next question if the event-dependency of the arrival rates can be the result of a rational strategy when customers only know the last realized event. For this purpose, we study the remaining service time distribution for an arriving customer given the last event. Next, we derive the expected waiting time of customers who arrive after an arrival and the one of customers who arrive after a service. We prove that the first one is always higher than the second one. Given a non-empty queue, a necessary and sufficient condition which involves the variability of the service time determines the comparison between the expected waiting times. As a conclusion, we show that λ -≥ λ + can be the result of a strategic behavior.

-Admission Control Problem. We investigate the historical admission control problem with eventdependency. A controller has to determine at the arrival of a customer whether to let this customer enter the queue or to reject this customer. The decision is based on the number of customers present in the system, the distribution of the remaining service time and the last realized event. The objective is to maximize the throughput of served customers with a service level constraint on the expected number of customers in the system. A two-thresholds policy depending on the last realized event and on the remaining service time of the customer in service is the only possible deterministic policy. Using a Markov decision process approach and approximating the service time distribution by a Coxian distribution allows us to find necessary conditions under which a deterministic threshold type policy is optimal.

These conditions are that λ -≥ λ + and that the departure rate out of each service phase of the Coxian distribution is decreasing in the number of remaining service phases. The relation between the optimal thresholds and the remaning service time might make the optimal policy complicate to implement in practice. We then propose a simplified version of this policy where the two thresholds depend only on the number of customers in the system. We finally develop an exact numerical method to obtain the performance measures under this policy at arbitrary instants and show that this policy outperforms the classical one-threshold policy.

The remainder of this paper is structured as follows. We conclude this section with a literature survey. In Section 2, we compute the performance measures and investigate the impact of the event-dependent arrival rates. In Section 3, we prove that the case λ -≥ λ + may be the result of a strategic behavior. In Section 4, we investigate the admission control problem. Finally, Section 5 opens on future research perspectives. The notations used are summarized at the end of the article.

Literature Review. Methodologically, the analysis of this paper is related to (i) queues with general service and state definition based on the residual service time [START_REF] Kerner | The conditional distribution of the residual service time in the Mn/G/1 queue[END_REF][START_REF] Sigman | Stationary remaining service time conditional on queue length[END_REF][START_REF] Boxma | Joint distribution of sojourn time and queue length in the M/G/1 queue with (in) finite capacity[END_REF], (ii) queueing systems with phase-type service time distributions [START_REF] Neuts | Matrix-Geometric Solutions in Stochastic Models: an Algorithmic Approach[END_REF] and (iii) Markov decision process approach [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF].

A stream of literature related to this paper is that of queueing systems with state-dependent parameters.

Single-server queues with workload-dependent arrival or service rates have been widely analyzed (e.g. see [START_REF] Kerner | The conditional distribution of the residual service time in the Mn/G/1 queue[END_REF][START_REF] Bekker | Queues with workload-dependent arrival and service rates[END_REF][START_REF] Boxma | On/off storage systems with state-dependent input, output, and switching rates[END_REF][START_REF] Harrison | The stationary distribution and first exit probabilities of a storage process with general release rule[END_REF]). Historically, Markovian models with general release rule have been considered as dam processes.

Later, another application came in packet-switched communication systems; the transmission rate of data connections can also be adapted based on the queue size. Other single-server queueing models have been proposed where the arrival or the service rates are depending on the waiting time of the customer in service [START_REF] Boxma | On queues with service and interarrival times depending on waiting times[END_REF][START_REF] Whitt | Queues with service times and interarrival times depending linearly and randomly upon waiting times[END_REF], the waiting time of the first customer in line [START_REF] Bekker | Queues with waiting time dependent service[END_REF], or on the remaining service time [START_REF] Legros | Stationary analysis of a single queue with remaining service time dependent arrivals[END_REF]. Other examples with workload dependent resources can be found via the slow server problem. With two servers (a slow and a fast one) [START_REF] Larsen | Control of a heterogeneous two-server exponential queueing system[END_REF], [START_REF] Lin | Optimal control of a queueing system with two heterogeneous servers[END_REF], [START_REF] Walrand | A note on optimal control of a queuing system with two heterogeneous servers[END_REF], and [START_REF] Koole | A simple proof of the optimality of a threshold policy in a two-server queueing system[END_REF] show that the fast server should be always used, and the slow server should be only used when the fast server is busy and the number of customers waiting in the queue exceeds a given threshold. Extension of these studies for more than two servers can be found in [START_REF] Legros | Routing in a queueing system with two heterogeneous servers in speed and in quality of resolution[END_REF], [START_REF] Rykov | Monotone control of queueing systems with heterogeneous servers[END_REF] and [START_REF] Luh | Threshold control policies for heterogeneous server systems[END_REF]. Our paper differs from these papers since the arrival-dependency is based on the last event instead of the observed workload or realized waiting time.

Performance Analysis

We investigate the impact of an event-dependent arrival process on the performance measures. In Section 2.1 we consider the embedded Markov chain right after a service completion. This approach, in line with the standard analysis of the M/G/1 queue, allows us to obtain explicitly the probability generating function.

Contrary to the M/G/1 queue, the queue length distribution is not identical at arrival and arbitrary instants. Therefore, in Section 2.2 we study the relation between these distributions in order to reach the performance measures at arbitrary instants. Finally in Section 2.3, we apply our results to particular service time distributions (Exponential, Erlang and Deterministic) to better understand the effect of the event-dependency of the arrival process. An alternative method to obtain the performance measures with a state definition based on the remaining service time is proposed in Section 1 of the online supplement.

The embedded Markov chain at service completion instants

Consider the system just after a customer has completed a service. The random variable X i represents the number of customers remaining in the system as the i th customer departs. We can write that

X i+1 = { X i -1 + A i+1 , if X i > 0, A i+1 , if X i = 0, (1)
where A i+1 is the number of customer who arrived during the service time of the (i + 1)st customer. The service time of the (i + 1)st customer is independent of previous service times and the length of the queue, so we can denote by S this random variable without mentioning the index of the (i + 1)st customer. We now evaluate the distribution of A i+1 . Two cases should be considered.

Case 1.

If X i > 0, the service initiation time of Customer i+1 is the service completion time of Customer i. Therefore, the last event for the first customer who arrives during the service of Customer i + 1 is a service.

It is an arrival for all the other customers. Let us denote by N t the number of customers who arrive during a service of length t. The distribution of N t is given by the following set of differential equations:

           P (N 0 = 0) = 1, dP (Nt=0) dt = -λ -P (N t = 0), dP (Nt=1) dt = -λ + P (N t = 1) + λ -P (N t = 0), dP (Nt=n) dt = -λ + P (N t = n) + λ + P (N t = n -1), for n ≥ 2.
(

) 2
After some algebra, we obtain the solution of this system. It is given by P (N t = 0) = e -λ -t , and

P (N t = n) = λ - λ + -λ - ( λ + λ + -λ - ) n-1 ( e -λ -t -e -λ + t n-1 ∑ k=0 ((λ + -λ -)t) k k! ) ,
for n ≥ 1. This arrival process is a modified Poisson process where the first interarrival time follows a different distribution than the other interarrival times. Given that X i > 0, the number of customers who arrive during the service of the (i + 1)st customer is independent of the index i. We therefore simply denote this random variable by A. We hence have

P (A i+1 = n|X i > 0) = P (A = n) = ∫ ∞ 0 P (N t = n)g(t) dt = α n ,
for i, n ≥ 0.

Case 2.

If X i = 0, the service initiation of Customer i +1 is the arrival time of Customer i+ 1. Therefore, the last event for all customers who arrive during the service of Customer i+1 is an arrival. Thus, the number of customer's arrivals during a service of length t given X i = 0 follows a Poisson process with rate λ + t and is independent of the index i. We then denote by B the random variable which represents the number of customers who arrive during the service of Customer i + 1 given that X i = 0. We may write

P (A i+1 = n|X i = 0) = P (B = n) = ∫ ∞ 0 e -λ + t (λ + t) n n! g(t) dt = β n ,
for i, n ≥ 0. This corresponds to the transition probabilities in the embedded Markov chain of an M/G/1 queue.

As a conclusion Equation (1) can be simplified into

X i+1 = { X i -1 + A, if X i > 0, and, B, if X i = 0. ( 3
)
The definition of the process in Equation (3) allows us to define a discrete time Markov chain. The related matrix of transition probabilities is given by

M =          β 0 β 1 β 2 β 3 β 4 • • • α 0 α 1 α 2 α 3 α 4 • • • 0 α 0 α 1 α 2 α 3 • • • 0 0 α 0 α 1 α 2 • • • 0 0 . . . . . . . . . . . .          .
Assuming that steady state is reached, we let p n represent the stationary probability that n customers are in the system at a service departure instant. We now define the generating functions

P (z) = ∞ ∑ n=0 p n z n , A(z) = ∞ ∑ n=0 α n z n , B(z) = ∞ ∑ n=0 β n z n , ( 4
)
for |z| ≤ 1. In Theorem 1, we derive P (z) and p 0 . In addition, we give the condition of existence of P (z).

This condition is also the condition which ensures the stationary regime.

Theorem 1 Under the stability condition

x < 1 λ -(1 -G * (λ -)) + 1 λ + G * (λ -), we have P (z) = p 0 (1 -z) ((λ + -λ -)G * (λ -) -λ + zG * (λ + -λ + z)) (1 -z) ((λ + -λ -)G * (λ -) -λ + z) + λ -z(1 -G * (λ + -λ + z))
,

with p 0 = λ + (1 -xλ -) + (λ --λ + )G * (λ -) λ + + (λ --λ + )G * (λ -) . Proof. The vector (p 0 , p 1 , • • • ) is solution of (p 0 , p 1 , • • • ) = (p 0 , p 1 , • • • ) × M . Therefore, we have p 0 = β 0 p 0 + α 0 p 1 , p 1 z = β 1 p 0 z + α 1 p 1 z + α 0 p 2 z, p 2 z 2 = β 2 p 0 z 2 + α 2 p 1 z 2 + α 1 p 2 z 2 + α 0 p 3 z 2 ,
. . .

p n z n = β n p 0 z n + α n p 1 z n + α n-1 p 2 z n + • • • + α 0 p n+1 z n , . . .
By summing up these equations, we get

P (z) = p 0 B(z) + A(z)(p 1 + p 2 z + p 3 z 2 + • • • )
. This leads to P (z) = p 0 B(z) + A(z) z (P (z) -p 0 ). Finally, we deduce that

P (z) = p 0 A(z) -zB(z) A(z) -z . ( 5
)
Using L'Hôpital's rule, we obtain

P (1) = 1 = p 0 A ′ (1) -B(1) -B ′ (1) A ′ (1) -1 . ( 6
)
From the results of the M/G/1 queue, we know that B ′ (1) = λ + x (e.g., see [START_REF] Kleinrock | Queueing Systems, Theory[END_REF], page 185). Moreover, by definition, B(1) = 1. From [START_REF] Kleinrock | Queueing Systems, Theory[END_REF] page 184, we know that B(z) = G * (λ + -λ + z).

There remains to compute A(z) and A ′ (1). For this purpose, let us introduce

N (z, t) = ∞ ∑ n=0 P (N t = n) • z n ,
for t > 0 and |z| < 1. We have

A(z) = ∫ ∞ 0 N (z, t) • g(t) dt.
From Equation (2), we get the differential equation followed by N (z, t);

dN (z, t) dt = (1 -z) ( -λ + N (z, t) + P (N t = 0) • (λ + -λ -)
) .

Using P (N t = 0) = e -λ -t and N (z, 0) = 1, we can solve this differential equation using the method of variation of constants. After some algebra, we obtain

N (z, t) = (1 -z)(λ + -λ -) λ + (1 -z) -λ -e -λ -t - λ -z λ + (1 -z) -λ -e -λ + (1-z)t .

This leads to

A(z) = (1 -z)(λ + -λ -) λ + (1 -z) -λ -G * (λ -) - λ -z λ + (1 -z) -λ -G * (λ + (1 -z)),
where G * (.) is the Laplace-Stieltjes Transform (LST) of the service time; G * (s) = ∫ ∞ t=0 g(t)e -st dt. Moreover,

A ′ (z) = λ -(λ + -λ -) (λ + (1 -z) -λ -) 2 (G * (λ -) -G * (λ + (1 -z))) + λ + λ -z λ + (1 -z) -λ -G ′ * (λ + (1 -z)). Using G * (0) = 1 and G ′ * (0) = -x, leads to A ′ (1) = λ + -λ - λ - ( G * (λ -) -1 ) + λ + x.
Replacing this expression in Equation [START_REF] Bekker | Queues with waiting time dependent service[END_REF] gives the expression of p 0 and the stability condition. Finally, replacing the expression of A(z) in Equation (5) allows us to obtain P (z) as in the theorem. 2

In Corollary 1, we deduce the expected number of customers in the system at departure instants, E(Q d ).

In the case λ + = λ -, we obtain the Pollaczeck-Kinchin formulas for the M/G/1 queue. When λ + ̸ = λ -, the performance measures presented here not only depend on the first two moment of the service time but also on the LST of the service time at λ -.

Corollary 1

We have

E(Q d ) = λ + x(G * (λ -)((λ -) 2 -(λ + ) 2 ) + (λ + ) 2 (1 -λ -x)) (λ + + (λ --λ + )G * (λ -))((λ --λ + )G * (λ -) + λ + (1 -λ -x)) + (λ + ) 2 λ -(x) 2 (1 + cv 2 ) 2((λ --λ + )G * (λ -) + λ + (1 -λ -x))
.

Proof. The expected number of customers in the system is equal to P ′ (1). The expression of P (z) allows us to write P (z) = N (z) 1-z D(z) . Therefore,

P ′ (z) = N ′ (z) 1-z D(z) -N (z) D(z)+(1-z)D ′ (z) D 2 (z)
. Using L'Hôpital's rule, one can derive the limit of 1-z D(z) as z tends to 1. For D(z)+(1-z)D ′ (z)

D 2 (z)
, L'Hôpital's rule should be used twice to obtain the limit of this expression as z tends to 1. This explains why the second moment of the service time is in the expression of E(Q d ). The details of the computation are omitted. 2

Performance analysis at arbitrary instants

We now relate the stationary probabilities at arbitrary instants with those at departure instants in order to use the results of Section 2.1 to obtain the performance of the system. We denote by π 0 the probability of an empty system at arbitrary instants and by π n,+ and π n,-the probability of having n customers in the system after an arrival and after a service respectively (n ≥ 1) at arbitrary instants.

The queue length distribution is identical at departure instants and arrival instants. The reason is that only one event (an arrival or a service departure) occurs at a time. We therefore have

p 0 = λ -π 0 λ - ∞ ∑ k=0 π k,-+ λ + ∞ ∑ k=1 π k,+ , ( 7
)
p n = λ -π n,-+ λ + π n,+ λ - ∞ ∑ k=0 π k,-+ λ + ∞ ∑ k=1 π k,+ , for n ≥ 1. (8)
Due to flow conservation, one may write

λ - ∞ ∑ k=0 π k,-+ λ + ∞ ∑ k=1 π k,+ = 1 x (1 -π 0 )
. So, we deduce from Equation ( 7) that π 0 = p0 p0+λx . Using now the result of Theorem 1, we deduce that

π 0 = λ + (1 -xλ -) + (λ --λ + )G * (λ -) λ + + (λ --λ + )(1 + xλ -)G * (λ -) . ( 9
)
Let us denote by p t (n, r, +) and p t (n, r, -) the probability-densities of having n customers in the system after an arrival and after a service respectively, n ≥ 1 and a remaining service time of r, r ≥ 0, at time t (given some arbitrary initial distribution). We also define the limit values of these probabilities; p(n, r, +) = lim t→∞ p t (n, r, +) and p(n, r, -) = lim t→∞ p t (n, r, -), for n ≥ 1. In Lemma 1, we provide the differential equations for p(n, r, +) and p(n, r, -), r ≥ 0 and n ≥ 1.

Lemma 1 For all r ≥ 0, p(n, r, +) and p(n, r, -) obey the following differential equations

p ′ (1, r, +) = λ + p(1, r, +) -λ -π 0 g(r), ( 10
) p ′ (n, r, +) = λ + p(n, r, +) -λ -p(n -1, r, -) -λ + p(n -1, r, +), for n ≥ 2, ( 11
)
p ′ (n, r, -) = λ -p(n, r, -) -g(r)(p(n + 1, 0, +) + p(n + 1, 0, -)), for n ≥ 1, ( 12
)
where p ′ (n, r, .) = dp(n,r,.) dr .

Proof. We will start with the case where n = 1 (equation ( 10)). First, observe that

p t+dt (1, r, +) = (1 -λ + dt)p t (1, r + dt, +) + λ -π 0 g(r)dt.
Taking t → ∞ and dividing by dt leads to

p(1, r + dt, +) -p(1, r, +) dt = λ + p(1, r + dt, +) -λ -π 0 g(r).
Next, taking dt → 0, we obtain Equation [START_REF] Boxma | On queues with service and interarrival times depending on waiting times[END_REF]. Equation ( 11) is derived from

p t+dt (n, r, +) = (1 -λ + dt) • p t (n, r + dt, +) + λ -dt • p t (n -1, r + dt, -) + λ + dt • p t (n -1, r + dt, +),
and Equation ( 12) from

p t+dt (n, r, -) = (1 -λ -dt) • p t (n, r + dt, -) + g(r)dt(p(n + 1, 0, +) + p(n + 1, 0, -)),
with the same approach. 2

In Proposition 1, using Lemma 1 we relate π n,+ and π n,-for n ≥ 1. This proposition proves that the ratio π n,+ /π n,-is constant for n ≥ 1. This translates that given a non-empty system, the queue length and the last event are independent. This result also holds for the M/G/1 queue.

Proposition 1 We have π n,-= λ + λ - 1-G * (λ -) G * (λ -) π n,+ , for n ≥ 1.
Proof. Integrating both sides of Equations ( 10), ( 11) and ( 12) for r from 0 to ∞, we get

p(1, 0, +) = λ -π 0 -λ + π 1,+ , ( 13
) p(n, 0, +) = λ -π n-1,-+ λ + π n-1,+ -λ + π n,+ , for n ≥ 2, ( 14
) p(n, 0, -) = p(n + 1, 0, +) + p(n + 1, 0, -) -λ -π n,-, for n ≥ 1. ( 15
)
Summing up Equations ( 14) and ( 15) yields

p(n, 0, +) + p(n, 0, -) = λ -π n-1,-+ λ + π n-1,+ -λ + π n,+ + p(n + 1, 0, +) + p(n + 1, 0, -) -λ -π n,-,
for n ≥ 2. This equation is equivalent to

p(n, 0, +) + p(n, 0, -) -λ -π n-1,--λ + π n-1,+ = p(n + 1, 0, +) + p(n + 1, 0, -) -λ -π n,--λ + π n,+ , for n ≥ 2. Since the difference p(n, 0, +) + p(n, 0, -) -λ -π n-1,--λ + π n-1,+
is not a function of n and tends to zero when n tends to ∞, we get the identity

p(n + 1, 0, +) + p(n + 1, 0, -) = λ -π n,-+ λ + π n,+ , (16)
for n ≥ 1. Combining Equation ( 16) with Equation ( 15) leads to

λ + π n,+ = p(n, 0, -), ( 17
)
for n ≥ 1.

Equation ( 12) can be written as

e -λ -u (p ′ (n, u, -) -λ -p(n, u, -)) = -g(u)(p(n + 1, 0, +) + p(n + 1, 0, -))e -λ -u .
Since the left hand side here is (e -λ -u p(n, u, -)) ′ , integrating both sides for u from 0 to ∞ yields

p(n, 0, -) = G * (λ -)(p(n + 1, 0, +) + p(n + 1, 0, -)), for n ≥ 1, (18)
Combining Equations ( 16) and ( 18) leads to

p(n + 1, 0, +) + p(n + 1, 0, -) = p(n, 0, -) G * (λ -) = λ -π n,-+ λ + π n,+ . ( 19
)
Since

λ + π n,+ = p(n, 0, -), we obtain π n,-= λ + λ - 1-G * (λ -) G * (λ -) π n,+ , for n ≥ 1. 2
Using the result of Proposition 1 in Equation ( 8), we obtain

π n,+ = (1 -π 0 ) G * (λ -) λ + x p n , and π n,-= (1 -π 0 ) 1 -G * (λ -) λ -x p n , ( 20
)
for n ≥ 1. This allows us to derive the generating functions

P + (z) = ∞ ∑ n=1 π n,+ • z n and, P -(z) = ∞ ∑ n=0 π n,-• z n .
After some algebra, we get

P + (z) = λ -G * (λ -)((λ --λ + )G * (λ -) + λ + (1 -λ -x)) λ + (λ + + G * (λ -)(λ --λ + )(1 + λ -x)) × z(G * (λ + -λ + z) -1)(λ + (1 -z) + λ -) (1 -z) ((λ + -λ -)G * (λ -) -λ + z) + λ -z(1 -G * (λ + -λ + z)) , P -(z) = (1 -G * (λ -))((λ --λ + )G * (λ -) + λ + (1 -λ -x)) λ + + G * (λ -)(λ --λ + )(1 + λ -x) × (1 -z) ((λ + -λ -)G * (λ -) -λ + zG * (λ + -λ + z)) (1 -z) ((λ + -λ -)G * (λ -) -λ + z) + λ -z(1 -G * (λ + -λ + z))
.

In Corollary 2, we deduce the expected number of customers in the system, E(Q), and the expected waiting time, E(W ), at arbitrary instants. We also derive the expected remaining service time seen by an arriving customer at a non-empty system, E(R).

Corollary 2

We have

E(Q) x(λ + + G * (λ -)(λ --λ + )) (21)
= xλ + λ -[ (1 + cv 2 )(λ + + G * (λ -)(λ --λ + )) -2λ + ] + 2 [ (λ -) 2 G * (λ -) + (λ + ) 2 (1 -G * (λ -)) ] 2 [λ + + (λ --λ + )(1 + xλ -)G * (λ -)] [-xλ + λ -+ λ + + G * (λ -)(λ --λ + )] , E(W ) = xλ + λ -[ (1 + cv 2 )(λ + + G * (λ -)(λ --λ + )) + 2(λ --λ + ) ] -2λ + (λ --λ + )(1 -G * (λ -)) 2λ -[-xλ + λ -+ λ + + G * (λ -)(λ --λ + )]
x, and, (

E(R) = (1 + cv 2 )x(λ + + G * (λ -)(λ --λ + )) 2λ - - (λ --λ + )(1 -G * (λ -) -λ -x) (λ -) 2 . ( 22)
) 23
Proof. The expected number of customers in the system is

E(Q) = ∞ ∑ n=1 n(π n,+ + π n,-). Using Equation (20),
we deduce that

E(Q) = 1 -π 0 x ( G * (λ -) λ + + 1 -G * (λ -) λ - ) ∞ ∑ n=1 np n = 1 -π 0 x ( G * (λ -) λ + + 1 -G * (λ -) λ - ) E(Q d ).
This leads to the expression of E(Q). The throughput of served customers is 1

x (1-π 0 ). Applying Little's law, we get the expected time spent in the system for a given customer by dividing the expected number in the system by 1

x (1-π 0 ). Finally, subtracting x to this expression leads to the expected waiting time in the queue. As mentioned above, the queue length distribution is identical at departure instants and arrival instants. Hence, the expected waiting time can be written as a function of E(R) using the stationary probabilities at arrival instants, p n , for n ≥ 1;

E(W ) = ∞ ∑ n=1 p n ((n -1)x + E(R)) . ( 24
)
This leads to

E(R) = x+ E(W )-xE(Q d ) 1-p0
. Using the explicit expressions of E(W ), E(Q d ) and p 0 , the expression of E(R) can be derived. 2

Special cases

Exponential case

In Proposition 2, we give the performance measures associated to an exponential service time distribution. 1 We denote by a + and a -the products

a + = x • λ + and a -= x • λ -.
Proposition 2 Under the stability condition a -a + < 1, we have

π 0 = 1 -a -a + 1 + a -, ( 25
)
E(W ) = x a + (1 + a -) 1 -a -a + , and (26)
P (W > t) = a + 1 + a - 1 + a + e -t x 1-a -a + 1+a + , ( 27
)
for t > 0.

In Figure 1, we illustrate the impact of λ -and λ + on the waiting time distribution. We observe that when the arrival rate after an arrival is higher than the arrival rate after a service then the system performance deteriorates in comparison with an M/M/1 queue. Most of the observed monotonicity results are intuitive since they correspond to the impact of the arrival parameter for an M/M/1 queue. The impact of a + is more surprising.

The probability P (W > t) can be concave in a + for t > 0. Note that we have lim

t→0 P (W > t) = a + 1 + a - 1 + a + .
The variable a -is only in the numerator of this expression. This explains the almost linear form of the curves of P (W > t) as functions of a -when t is small (Figure 1(a)). 1 The performance measures in the exponential case could also be obtained using a Markov chain analysis.

Erlang case

In Proposition 3, we give the performance measures associated to an Erlang service time distribution with N exponential phases where each phase has an expected duration of x N . 2

Proposition 3 Under the stability condition a -a + < a -

( N a -+N ) N + a + ( 1 - ( N a -+N ) N
) , we have

π 0 = a + (1 -a -) + (a --a + ) ( N N +a - ) N a + + (a --a + )(1 + a -) ( N N +a -
) N , and

E(W ) = x a + a - ( ( 1 N + 1 ) ( (a --a + ) ( N N +a - ) N + a + ) + 2 (a --a + ) ) -2a + (a --a + ) ( 1 - ( N N +a - ) N ) 2a - ( (a --a + ) ( N N +a - ) N -a + a -+ a +
) .

The impact of N on the stability of the system depends on the difference λ + -λ -. Consider the function

F (N ) = ( N a -+N ) N (a --a + ) + a + (1 -a -). We have F (N ) > 0 if and only if a -a + < a - ( N a -+N ) N + a + ( 1 - ( N a -+N ) N
) . So, the system is stable if and only if F (N ) > 0. Since the function

( N a -+N ) N is decreasing in N , F is decreasing in N if and only if λ + < λ -.
In other words, increasing N stabilizes the system if and only if λ + > λ -.

Finally, we can conclude that the system is stable for all Erlang service time distributions with expected

duration x if and only if { a -a + < 1 a + ≥ a -or { e -a -> (a --1)a + a --a + a + < a - .
These relations are obtained either by choosing N = 1 (exponential distribution) or by letting N → ∞ (deterministic distribution).

We now evaluate the impact of the number of phases on E(W ) for different values of λ + and λ -(Figure 2). We observe that the number of phases can deteriorate the expected waiting time when λ -> λ + (Figure 2(a)). This result is surprising since it is in contradiction with the improvement which should result from the reduction of the variability in the service process when the number of phases increases. The reason is related to the presence of the term G * (λ -) in the expression of E(W ) which involves the service time distribution. The improvement related to the decreasing of the coefficient of variability when N increases can be compensated by an increasing number of arrivals during service.

The expected number of customers arriving during a service has already been computed in Section 2.1.

It is either equal to A ′ (1) or to B ′ (1). The expression of A ′ (1) explains the position of the curves in Figure 2. Since G * (λ -) decreases as N increases, the expected number of arrivals during service increases as N increases in the case λ -> λ + . This is consistent with the conclusion derived above for the stability region;

if λ -> λ + , increasing N reduces the stability region.

2 Another way to compute the performance measures in the Erlang case is to use a Matrix geometric approach.

Deterministic case

By letting r → ∞ in Proposition 3 we obtain the performance measures in the case of a deterministic service with duration x.

Proposition 4 Under the stability condition a -a + < a + + e -a -(a --a + ), we have

π 0 = a + (1 -a -) + (a --a + )e -a - a + + (a --a + )(1 + a -)e -a -
, and

E(W ) = xa + (a -) 2 + (a --a + )(a --2 + (a -+ 2)e -a -) 2a - ( (a --a + ) e -a --a + a -+ a +
) .

Equilibrium Strategy

An interesting question is to determine if the event-dependency of the arrival process can be the result of an individual rational strategy. 3 For this purpose, we are interested in the service level differentiation after an arrival or after a service. The idea is to determine how does the nature of the last event give an indication on the expected waiting time encountered by an arriving customer. If the last event is a service, an arriving customer may arrive in an empty system. This cannot be the case if the last event is an arrival. This gives an advantage to the customers who arrive after a service. However, given a non-empty system, if the last event is a service, an arriving customer is the first to arrive during a service time. Customers who arrive after an arrival arrive later during the service time. The order of arrival during a service may also influence the expected remaining service time seen by an arriving customer. For instance, with a deterministic service time, the first customer who arrives during a service has a longer expected remaining service time than the next customers. In order to compare the expected waiting time after an arrival or a service, we first evaluate in Section 3.1 the expected remaining service times after an arrival or a after service seen by an arriving customer. In Section 3.2, we derive the expected waiting time of a customer who arrives after an arrival, and the one of a customer who arrives after a service, and compare between them.

Expected remaining service time

We are interested in the expected remaining service of an arriving customer. This metric is function of the last realized event and of the number of customers present in the system. We denote by r n,+ and r n,-the expected remaining service time seen by a customer who arrives after an arrival or a service respectively when n customers are present in the system, for n ≥ 1. We denote by P * (n, s, +) and P * (n, s, -)

the Laplace-Stieltjes Transform (LST) of p(n, r, +) and p(n, r, -); P * (n, s, +) = ∫ ∞ 0 p(n, r, +)e -sr dr, and P * (n, s, +) = ∫ ∞ 0 p(n, r, +)e -sr dr. We have

r n,+ = 1 π n,+ ∫ ∞ 0 rp(n, r, +) dr = -lim s→0 P ′ * (n, s, +) π n,+ ,
and,

r n,-= 1 π n,- ∫ ∞ 0 rp(n, r, -) dr = -lim s→0 P ′ * (n, s, -) π n,- ,
for n ≥ 1. In Proposition 5, we provide a recursive formula for the conditional distributions of the residual service times. These relations allows us to compute r n,+ and r n,-for n ≥ 1. In addition, we show in Section 1 of the online supplement how these relations can be used to obtain the performance measures already found in Section 2.

Proposition 5 Under the stability condition

x < 1 λ -(1 -G * (λ -)) + 1 λ + G * (λ -),
we have the following initial relations:

π 1,+ = λ - λ + (1 -G * (λ + ))π 0 , ( 28
)
P * (1, s, +) = λ + 1 -G * (λ + ) G * (λ + ) -G * (s) s -λ + π 1,+ . ( 29
)
Next, for n ≥ 1,

P * (n, s, -) = λ - 1 -G * (λ -) G * (λ -) -G * (s) s -λ - π n,-, (30)
π n+1,+ = π n,+ G * (λ -) - λ - λ + P * (n, λ + , -) -P * (n, λ + , +), (31)
P * (n + 1, s, +) = λ -P * (n, λ + , -) -P * (n, s, -) s -λ + + λ + P * (n, λ + , +) -P * (n, s, +) s -λ + . ( 32
)
Proof. First, Equation ( 10) can be written as

e -λ + u (p ′ (1, u, +) -λ + p(1, u, +)) = -λ -e -λ + u π 0 g(u).
Since the left hand side here is (e -λ + u p(1, u, +)) ′ , integrating both sides for u from r to ∞ leads to

p(1, r, +) = λ -π 0 e λ + r ∫ ∞ u=r e -λ + u g(u) du. ( 33
)
Inserting r = 0 in (33) leads to

p(1, 0, +) = λ -π 0 G * (λ + ). ( 34
)
Combining now Equation ( 13) of Section 2 and Equation (34), leads to

π 1,+ = λ - λ + (1 -G * (λ + ))π 0 (Equation (28)).
From Equation [START_REF] Boxma | On queues with service and interarrival times depending on waiting times[END_REF], one may write

∫ ∞ r=0 e -sr p ′ (1, r, +) dr = λ + ∫ ∞ r=0 e -sr p(1, r, +) dr -λ -π 0 ∫ ∞ r=0 e -sr g(r) dr.
This equation leads to

sP * (1, s, +) -p(1, 0, +) = λ + P * (1, s, +) -λ -π 0 G * (s). Since p(1, 0, +) = λ -π 0 G * (λ + ) (Equation (34)), we obtain P * (1, s, +) = λ -G * (λ + )-G * (s) s-λ +
sP * (n, s, -) -p(n, 0, -) = λ -P * (n, s, -) -(p(n + 1, 0, +) + p(n + 1, 0, -))G * (s). ( 35
)
Combining Equation ( 19) and Proposition 1 leads to p(n + 1, 0, +) + p(n + 1, 0, -)

= λ + πn,+ G * (λ -)
. Using now Equation ( 17) with Proposition 1 in Equation (35), we obtain Equation [START_REF] Walrand | A note on optimal control of a queuing system with two heterogeneous servers[END_REF].

With the same approach as for Equations ( 33) and (34), we obtain

p(n, r, +) = e λ + r ∫ ∞ u=r e -λ + u (λ -p(n -1, u, -) + λ + p(n -1, u, +)) du, for n ≥ 2, ( 36
) p(n, 0, +) = λ -P * (n -1, -, λ + ) + λ + P * (n -1, +, λ + ), for n ≥ 2, ( 37
)
using Equation [START_REF] Harrison | The stationary distribution and first exit probabilities of a storage process with general release rule[END_REF]. Combining Equation (37) and Equation ( 14) leads to

λ + π n,+ = λ -π n-1,-+ λ + π n-1,+ -(λ -P * (n -1, -, λ + ) + λ + P * (n -1, +, λ + )),
for n ≥ 2. With Equation ( 19), λ -π n-1,-+ λ + π n-1,+ = p(n, 0, +) + p(n, 0, -). Combining Equation [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF] with Equation ( 17) leads to p(n, 0, +) + p(n, 0, -) = p(n-1,0,-)

G * (λ -) = λ + πn-1,+ G * (λ -)
. The aforementioned relations prove Equation [START_REF] Whitt | Queues with service times and interarrival times depending linearly and randomly upon waiting times[END_REF].

From Equation (11), we have

∫ ∞ r=0 e -sr p ′ (n, r, +) dr =λ + ∫ ∞ r=0 e -sr p(n, r, +) dr -λ - ∫ ∞ r=0 e -sr p(n -1, r, -) dr -λ + ∫ ∞ r=0 e -sr p(n -1, r, +) dr,
for n ≥ 2. The equality is equivalent to

sP * (n, s, +) -p(n, 0, +) = λ + P * (n, s, +) -λ -P * (n -1, s, -) -λ + P * (n -1, s, +).
Finally, Equation (37) leads to Equation (32).

2

The consequence of Equation ( 30) is that the expected remaining service time of a customer who arrives in a non-empty system after a service does not depend on the system size. This property allows us to explicitly derive in Corollary 3 the expected remaining service time seen at arrival for a customer who arrives after a service, E(R -), and after an arrival, E(R + ).

Corollary 3

We have

E(R -) = x 1 -G * (λ -) - 1 λ -, and, ( 38
)
E(R + ) = λ -x(1 + cv 2 )(λ + (1 -G * (λ -)) + λ -G * (λ -)) -2λ + (λ -x + G * (λ -) -1) 2(λ -) 2 G * (λ -) . ( 39
)
Proof. We have r n,-= -lim s→0 P ′ * (n,s,-) πn,-, for n ≥ 1. Using Equation (30), we get r n,-=

x 1-G * (λ -) -1 λ -, for n ≥ 1. Therefore, the expected remaining service time seen by an arriving customer after a service does not depend on n. Hence, E(R -) = r n,-. The same property does not hold for r n,+ . However, we can compute E(R + ) using the following decomposition:

E(W ) = ∞ ∑ n=1 ((n -1)x + E(R + ))p n,+ + ∞ ∑ n=1 ((n -1)x + E(R -))p n,-, (40)
where p n,+ and p n,-are the stationary probabilities to have n customers in the system at arrival instants after an arrival or a service (p n = p n,+ + p n,-, for n ≥ 1). This leads to

E(R + ) = x + E(W ) -xE(Q d ) -(E(R -) -x) ∞ ∑ n=1 p n,- ∞ ∑ n=1 p n,+
.

All parts of this expression are know except

∞ ∑ n=1 p n,-and ∞ ∑ n=1 p n,+ .
In what follows we derive these metrics.

We have

p n,-= λ -πn,- λ - ∞ ∑ n=0 πn,-+λ + ∞ ∑ n=1 πn,+
. The denominator of this expression is 1-π0 x due to flow conservation.

Next, from Equation (20), we obtain

∞ ∑ n=1 π n,-= (1 -π 0 ) 1 -G * (λ -) λ -x ∞ ∑ n=1 p n .
Therefore,

∞ ∑ n=1 p n,-= λ -(1-π0) 1-G * (λ -) λ -x (1-p0) 1-π 0 x
. Using the expression of p 0 in Theorem 1, we obtain

∞ ∑ n=1 p n,-= λ + λ -x(1 -G * (λ -)) λ -G * (λ -) + λ + (1 -G * (λ -))
.

With the same approach, we obtain

∞ ∑ n=1 p n,+ = λ + λ -xG * (λ -) λ -G * (λ -)+λ + (1-G * (λ -))
. This finishes the proof. 2

Remark. The independence between the system size and the remaining service time for a customer who arrives after a service is not surprising. As mentioned above, a customer who arrives after a service in a non-empty system is the first to arrive during a service irrespective of the number of customers already present. For a customer who arrives after an arrival, the number of customers present in the system has an impact. For instance, if a customer arrives after an arrival when one customer is already in the system, then it means that the customer in service arrived in an empty system. So, r 1,+ is the expected remaining service time of the first customer who arrives during a service. As expected, with Equation (29), we obtain

r 1,+ = x 1-G * (λ + ) -1 λ +
. This is exactly the expression of E(R -) by replacing λ -by λ + . Consider now the same situation when two customers are already in the system. In this case the arrived customer cannot be the first to arrive during a service, so the expression of r 2,+ may differ from the one of r 1,+ .

Expected waiting time

In Theorem 2, we compare between the expected waiting times after an arrival and after a service. Let us denote by E(W + ), E(W -), and E(W -|non-empty) the expected waiting times of a customer who arrives after an arrival, after a service, and after a service given a non-empty system.

Theorem 2

The following holds.

E(W + ) ≥ E(W -|non-empty) if and only if

cv 2 +1 2 ≥ E(R -) x , 2. E(W + ) ≥ E(W -).
Proof. Let us start with the first statement. Proposition 1 proves that the ratio π n,+ /π n,-is constant

for n ≥ 1. Since p n,+ = λ + πn,+ λ - ∞ ∑ n=0 πn,-+λ + ∞ ∑ n=1 πn,+ and p n,-= λ -πn,- λ - ∞ ∑ n=0 πn,-+λ + ∞ ∑ n=1 πn,+
, the ratio p n,+ /p n,-is also constant for n ≥ 1. This translates that given a non-empty system, the queue length at customers' arrival and the last event are independent. This result is important for the comparison between the expected waiting times after an arrival or after a service given a non-empty system. It means that they only differ in their remaining service times. This can also be shown via the following decomposition:

E(W + ) = ∞ ∑ n=1 p n,+ ((n -1)x + E(R + )) ∞ ∑ n=1 p n,+ = E(R + ) -x + E(Q d ) λ + (1 -G * (λ -)) + λ -G * (λ -) λ + λ -
, and,

E(W -|non-empty) = ∞ ∑ n=1 p n,-((n -1)x + E(R -)) ∞ ∑ n=1 p n,- = E(R -) -x + E(Q d ) λ + (1 -G * (λ -)) + λ -G * (λ -) λ + λ - .
Therefore, by comparing E(R + ) and E(R -) using their expressions in Corollary 3, we get the condition of the first statement.

Let us now consider the second statement. The order of arrival during a service determines the comparison between E(W + ) and E(W -|non-empty) (first statement). The comparison between E(W + ) and E(W -) also involves the probability to arrive in an empty system. This makes the comparison more complex since two phenomena are involved. In this case, the explicit expressions of the expected waiting times are required for the comparison. These expressions are computed in a similar way as E(R + ) and E(R -) in Corollary 3. We have

E(W + ) = [ λ + (1 -G * (λ -)) + λ -G * (λ -) ] [ λ -x(1 + cv 2 )(λ + (1 -λ -x) + G * (λ -)(λ --λ + + λ + λ -x)) -2λ + (1 -λ -x)(λ -x + G * (λ -) -1) ] 2(λ -) 2 G * (λ -)(λ + (1 -λ -x) + (λ --λ + )G * (λ -)) , E(W -|non-empty) = [ λ + (1 -G * (λ -)) + λ -G * (λ -) ] [ λ -λ + x 2 (1 -G * (λ -))(1 + cv 2 ) + 2(1 -λ + x)(λ -x + G * (λ -) -1) ] 2λ -(1 -G * (λ -))(λ + (1 -λ -x) + (λ --λ + )G * (λ -))
, and,

E(W -) = λ + x [ λ + (1 -G * (λ -)) + λ -G * (λ -) ] [ λ -λ + x 2 (1 -G * (λ -))(1 + cv 2 ) + 2(1 -λ + x)(λ -x + G * (λ -) -1) ] 2(λ + (1 -λ -xG * (λ -)) + (λ --λ + )G * (λ -))(λ + (1 -λ -x) + (λ --λ + )G * (λ -)) .
One then may write

E(W + ) -E(W -) = (λ + (1 -G * (λ -)) + λ -G * (λ -))(λ -x(1 + cv 2 )(λ + (1 -G * (λ -)) + λ -G * (λ -)) + 2λ + (1 -λ -x -G * (λ -))) 2(λ -) 2 G * (λ -)(λ + (1 -G * (λ -)) + λ -G * (λ -)(1 -λ + x))
The sign of the denominator depends on the sign of λ

+ (1 -G * (λ -)) + λ -G * (λ -)(1 -λ + x). This expression can be rewritten as λ + λ -x(1 -G * (λ -)) + λ + (1 -xλ --G * (λ -)) + λ -G * (λ -). Since G * (λ -) ≤ 1, λ + λ -x(1 - G * (λ -)) ≥ 0. Next, we have 1 -xλ --G * (λ -) = ∫ ∞ 0 g(t)(1 -λ -t -e -λ -t ) dt. Let us define ϕ 1 (t) = 1 -λ -t -e -λ -t , for t ≥ 0. We have ϕ ′ 1 (t) = -λ -(1 -e -λ -t ) ≤ 0. So, ϕ 1 (t) is decreasing in t. Moreover, ϕ 1 (0) = 0 then ϕ 1 (t) ≤ 0, for t ≥ 0. This proves that 1 -xλ --G * (λ -) ≤ 0. Hence, λ + (1 -xλ --G * (λ -)) + λ -G * (λ -) is decreasing in λ + . The stability condition in Theorem 1 is equivalent to λ + < λ -G * (λ -) xλ -+ G * (λ -) -1 . ( 41
)
As λ + tends to

λ -G * (λ -) xλ -+G * (λ -)-1 , λ + (1 -xλ --G * (λ -)) + λ -G * (λ -) tends to 0. This proves that λ + (1 - G * (λ -)) + λ -G * (λ -)(1 -λ + x) ≥ 0. Consider now the numerator. Since 0 ≤ G * (λ -) ≤ 1, λ + (1 -G * (λ -)) + λ -G * (λ -) ≥ 0. Next, the inequality λ -x(1 + cv 2 )(λ + (1 -G * (λ -)) + λ -G * (λ -)) + 2λ + (1 -λ -x -G * (λ -)) ≥ 0 is equivalent to cv 2 + 1 2 ≥ λ + (G * (λ -) + λ -x -1) λ -x(λ -G * (λ -) + λ + (1 -G * (λ -)) . ( 42
)
We define the function in λ + , ψ(λ

+ ) = λ + (G * (λ -)+λ -x-1) λ -x(λ -G * (λ -)+λ + (1-G * (λ -)) . We have ψ ′ (λ + ) = G * (λ -)(G * (λ -)+λ -x-1) x(λ -G * (λ -)+λ + (1-G * (λ -)) 2 ≥ 0, since we have G * (λ -) + λ -x -1 ≥ 0. So ψ is increasing in λ + . As λ + tends to λ -G * (λ -) xλ -+G * (λ -)-1 (upper bound for λ + , see Equation (41)), ψ(λ + ) tends to G * (λ -)+λ -x-1 (λ -x) 2 . Therefore, ψ(λ + ) ≤ G * (λ -)+λ -x-1 (λ -x) 2 . Consider now cv 2 +1 2 -G * (λ -)+λ -x-1 (λ -x) 2 = 1-λ -x+ (λ -) 2 2 E(S 2 )-G * (λ -) (λ -x) 2
, where E(S 2 ) = ∫ ∞ 0 t 2 g(t) dt. We have

1 -λ -x + (λ -) 2 2 E(S 2 ) -G * (λ -) = ∫ ∞ 0 ( 1 -λ -t + (λ -) 2 2 t 2 -e -λ -t ) g(t) dt. We define ϕ 2 (t) = 1 -λ -t + (λ -) 2 2 t 2 -e -λ -t , for t ≥ 0. We have ϕ ′ 2 (t) = -λ -ϕ 1 (t) ≥ 0. So ϕ 2 (t) is increasing in t. Moreover, ϕ 2 (0) = 0. So, ϕ 2 (t) ≥ 0, for t ≥ 0 and 1 -λ -x + (λ -) 2 2 E(S 2 ) -G * (λ -) ≥ 0. This proves that 0 ≤ cv 2 +1 2 -G * (λ -)+λ -x-1 (λ -x) 2 ≤ cv 2 +1
2 -ψ(λ + ). Hence, Inequality (42) holds in all cases. This finishes the proof of the second statement.

2

For the first statement, the condition cv 2 +1

2

≥ E(R -)
x reveals the importance of the service variability in the performance comparison. With high variability, the indication that the last event was an arrival is a signal that the expected waiting time may be longer than if the last event was a service.

Examples.

1. The service time follows an exponential distribution. Then,

E(R -) = x and cv = 1. So E(W + ) = E(W -|non-empty).
2. The service time follows a deterministic distribution. Then E(R -)

x = 1
1-e -λ -x -1 λx and cv = 0. We can show that

1 1-e -λ -x -1 λ -x ≥ 1/2. This proves that E(W + ) ≤ E(W -|non-empty). 3.
The service time follows a particular hyperexponential distribution for which a customer is either served with an exponential duration with rate µ 0 with probability q or is instantaneously served with probability

1 -q. We have G * (λ -) = qµ0 µ0+λ -, x = q µ0 and cv 2 = 2 q -1. Then, cv 2 +1 2 -E(R -) x = (1-q)( µ 0 λ -+1) 2 q(1+ µ 0 λ -(1-q)) ≥ 0. So, E(W + ) ≥ E(W -|non-empty).
The second statement indicates that in all cases, the expected waiting time of an arriving customer after an arrival is longer than the expected time after a service. This might be surprising. Given the first statement, one could imagine that some counterexamples could be found for instance in cases with a low service variability and a high workload situation.

We now explain how the model studied in the article can be the result of a strategic behavior. Consider an initial unobservable system with a potential arrival rate λ. At arrival, a customer can decide to join the queue or not to join. Although the system is unobservable, an arriving customer is informed about the last realized event; an arrival or a service completion. We define the net benefit for a customer who joins by the value of service, B, minus the cost of waiting proportional to a waiting cost per time unit, C. Given that the expected waiting time of an arriving customer is different whether the last event was an arrival or a service, a strategy after an arrival or after a service can be described by two different probabilities of joining p + and p -, such that λ + = p + • λ and λ -= p -• λ. The possible values for p + and p -are such that the stability condition in Theorem 1 is satisfied. The net benefit for a customer who arrives after and arrival is therefore

B -C • E(W + ) and it is B -C • E(W -)
for a customer who arrives after a service completion. We have the following cases:

-Case 1: E(W + ) ≤ B C for p + = p -= 1. In this case, since E(W -) ≤ E(W + ), we also have E(W -) ≤ B C
for p + = p -= 1. In this case, even if all potential customers join after an arrival or a service, they all enjoy a non-negative benefit. Therefore, the strategy of joining with probability p + = p -= 1 is an equilibrium strategy. 4

-Case 2: E(W + ) ≥ B C as p + tends to 0. As p + tends to 0, E(W + ) tends to x 1+cv 2 2 and E(W -) tends to 0. 5 Even if no other customer joins after an arrival, the net benefit of a customer who joins after an arrival is non-positive. Therefore, the strategy of joining after an arrival with probability p + = 0 is an equilibrium strategy and no other equilibrium is possible after an arrival. If all customers join after a service completion, they all enjoy a non-negative benefit since E(W -) = 0. Therefore the strategy of joining with p + = 0 and p -= 1 is an equilibrium strategy.

-Case 3: E(W + ) > B C for p + = p -= 1 and E(W + ) < B C as p + tends to 0. In this case if p + = 1 then a customer who joins after an arrival suffers a negative benefit. This cannot be an equilibrium strategy. If p + = 0, then all customers balk after an arrival. Yet, a customer who joins after an arrival would get a positive benefit. This contradiction shows that p + = 0 cannot be an equilibrium strategy. There exists a unique equilibrium where λ

+ solves E(W + ) = B C for p -= 1. In this case E(W -) ≤ B C , therefore if 4 In this case E(W + ) = x 1+cv 2 2 ( 1 G * (λ) + λx 1-λx ) + 1-λx λG * (λ) -1 λ , and E(W -) = x ( 1+cv 2 2 (λx) 2 (1-G * (λ)) (1-λx)(1-λxG * (λ)) + G * (λ)+λx-1 1-λxG * (λ)
)

. 5 The limit of E(W + ) corresponds to the expected remaining service time in an M/G/1 queue.

all customers join after a service, they all enjoy a non-negative benefit. Therefore the strategy of joining with 0 < p + < 1 and p -= 1 is also an equilibrium strategy.

As a conclusion, the case λ -≥ λ + can be seen as the result of a strategic behavior.

The Admission Control Problem

We question here the possibility to take different decisions after a service or after an arrival to answer a classical routing problem in the queueing theory. This problem is referred in the literature as the admission control problem (Section 1 in [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF]).

The optimisation problem

We propose to solve this problem under event-dependency. As mentioned in the previous section, the eventdependency may be the result of a customer strategic behavior when only the information of the last event is given. This leads to λ -≥ λ + . However, it might be interesting for the system to better control the arrival process by accepting or rejecting customers based on the system size. Let us specify that the arrival of a customer who was rejected is not seen as an arrival for the purpose of the rate changing.

A controller has to determine at arrival of a new customer whether we allow this new customer to enter the system or whether we reject this customer from the system. The optimization problem may be written as

{ Maximize T S , subject to E(Q) ≤ E(Q), ( 43
)
where T S is the throughput of served customers, E(Q) is the expected number of customers in the queue and, E(Q) is the service level constraint on E(Q). We are restricting the class of admissible policies to the class of deterministic policies. In real system, for instance in a shop, deterministic policies are easier to implement than non deterministic ones which may require randomization. Yet, deterministic policies are not necessarily optimal. In order to saturate the constraint, it may be useful to randomize between two or more deterministic policies. This however may only improve the class of deterministic policies. It does not lead to the optimal policy. 6 In Section 4.2, we give conditions under which deterministic policies (or randomization between a number of them) are optimal.

Let us now specify the nature of a deterministic policy. Any policy within the class of deterministic stationary policies is equivalent to a two-thresholds policy for a given remaining service time. More precisely, for a given remaining service time r and k customers in the system, a two-thresholds policy is defined by two

thresholds k + r and k - r (k + r , k - r ≥ -1) such that,
To prove that the optimal policy has a threshold structure, we need the conditions λ -≥ λ + , and, (

(1 -r x )µ x ≥ (1 -r x+1 )µ x+1 , for 1 ≤ x ≤ N. ( 45)
) 46
In Proposition 6, under Conditions (45) and ( 46), we prove by induction on the value function that the optimal policy is of threshold type based on the number of customers in the system for a given number of remaining phases. The proof follows a standard MDP methodology where structural properties of the value function are proven by induction. The complete detailed proof is given in Section 2 of the online supplement. Condition ( 45) is required to show the convexity and the supermodularity properties of the value function.

Condition (46) means that the departure rate out of a given service phase increases with the number of elapsed phases of service. This condition allows the value function to be increasing in the number of remaining service phases. This monotonicity property of the value function is required to prove that the optimal policy has a threshold structure.

Without this condition, the value function can be non-increasing in the number of remaining phases of service. Consider for instance a situation with a Coxian service time distribution with two phases. The first phase of service has an expected duration of 1 and the second one has an expected duration of 100. The probability to end service after the first waiting phase is 90%. So, if at a customer arrival there is 1 customer in the system being in the second phase of service, then the expected waiting time of an arriving customer is 100. If at a customer arrival there are 2 customers in the system and the one in service being in the first phase of service, then the expected waiting time of an arriving customer is 2 × (1 + 0.1 × 100) = 22 < 100.

Therefore in this case, the congestion of the system is not positively bound to the number of remaining phases of service and the value function may not be increasing.

Numerical illustration. We consider a particular Coxian distribution for the service time with two exponential phases with rate µ 1 and µ 2 and with probabilities r 1 = 0 and r 2 = 1. 7 From proposition 6, the optimal policy is determined by the thresholds k + 1 , k - 1 , k + 2 , and k - 2 . The computation of the performance measures can be done using a Matrix geometric approach (e.g., see [START_REF] Neuts | Matrix-Geometric Solutions in Stochastic Models: an Algorithmic Approach[END_REF]). In Figures (3(

= k - 1 + k - 2 and k + = k + 1 + k 2 2
for three different values for k + -k -. The parameters k -and k + represent the thresholds on the number of remaining service phases. In this example, the optimal thresholds to answer the optimization problem are k + = 3 and k -= 2. This means that an arriving customer should be rejected after an arrival if there is strictly more than 2 customers in the system or if there is two customers in the system and the customer in service still has 2 phases of service to achieve. A customer should be rejected after a service if there is strictly more than one customer in the system. This illustrates a case where the two-thresholds policy achieves a better solution to the optimization problem than the one-threshold policy.

Performance analysis under the two-thresholds policy

The result of Proposition 6 may be difficult to implement in practice. First, it might be difficult to find the appropriate Coxian distribution which well approximates the considered service time distribution. A field of research is dedicated to this problem (e.g., see [START_REF] Van Der Heijden | On the three-moment approximation of a general distribution by a coxian distribution[END_REF], [START_REF] Asmussen | Fitting phase-type distributions via the em algorithm[END_REF], [START_REF] Horváth | Phfit: A general phase-type fitting tool[END_REF]). Second, the optimal policy depends on the remaining number of service phases or more generally on the remaining service time. This information might also be complicate to obtain in practice.

To overcome these difficulties, we propose a numerical analysis to obtain the performance measures under a two-thresholds policy where the thresholds do not depend on the remaining service time. Although the proposed policy is not optimal, it is simple to implement and may lead to enhanced performance compared to a one-threshold policy based on the number of customers in the system. Moreover, compared to the Coxian approximation followed by a Matrix geometric approach, the method developed here leads to the exact performance measures for any service time distribution.

We denote by k + and k -the thresholds on the system size after an arrival or a service. After an arrival (respectively a service) customers are rejected if there is strictly more than

k + (respectively k -) customers      β 0 β 1 β 2 β 3 β 4 • • • • • • • • • 1 - k + -1 ∑ n=0 β n α 0 α 1 α 2 α 3 α 4 • • • • • • • • • 1 - k + -1 ∑ n=0 α n 0 α 0 α 1 α 2 α 3 • • • • • • • • • 1 - k + -2 ∑ n=0 α n 0 0 α 0 α 1 α 2 • • • • • • • • • 1 - k + -3 ∑ n=0 α n . . . . . . . . . . . . . . . . . . • • • • • • . . . 0 • • • • • • • • • α 0 α 1 • • • • • • 1 - k + -k - ∑ n=0 α n 0 • • • • • • • • • 0 1 0 • • • 0 0 • • • • • • • • • 0 0 1 • • • 0 . . . . . . . . . . . . . . . . . . . . • • • . . . 0 • • • • • • • • • • • • • • • • • • 1 0                                
.

The first k -+ 1 lines of the matrix are identical to those in the matrix in Section 2.1 except that the number of customers after a service completion is bounded by k + . From line k -+ 2 to line k + + 1, there cannot be any arrival during a service because the first arrival after the service completion cannot occur. So, the number of customers is forced to be reduced by one at the next service completion. One can then solve the system P = P • M , where P = (p 0 , p p n = 1 to obtain the stationary distribution of the system size at departure instants. In Proposition 7, we relate these probabilities to those at arbitrary instants.

Proposition 7

We have

π n,+ = (1 -π 0 ) G * (λ -) λ + x p n , and π n,-= (1 -π 0 ) 1 -G * (λ -) λ -x p n , for 1 ≤ n ≤ k -, and
π n,+ = 1 -π 0 λ + x p n , π n,-= (1 -π 0 )p n , and π k + +1,+ = 1 -   π 0 + k + ∑ n=1 π n,+ + π n,-   , for k -+ 1 ≤ n ≤ k + .
Proof. Contrary to Section 2.2, the system is not identical at arrival instants and departure instants but it is identical at arrival instants of customers who join and departure instants. Therefore, we have

p 0 = λ -π 0 λ - k - ∑ k=0 π k,-+ λ + k + ∑ k=1 π k,+ , ( 47
)
p n = λ -π n,-+ λ + π n,+ λ - k - ∑ k=0 π k,-+ λ + k + ∑ k=1 π k,+ , for 1 ≤ n ≤ k -, and (48)
p n = λ + π n,+ λ - k - ∑ k=0 π k,-+ λ + k + ∑ k=1 π k,+ , for k -+ 1 ≤ n ≤ k + . ( 49
)
Due to flow conservation, one may write λ -

k - ∑ k=0 π k,-+ λ + k + ∑ k=1 π k,+ = 1 x (1 -π 0 )
. So, we deduce from Equation (47) that π 0 = p0 p0+λx . The throughput of served customers is hence λ - p0+λx . Similarly to the proof of Lemma 1, p(n, r, +) and p(n, r, -) obey the following differential equations:

p ′ (1, r, +) = λ + p(1, r, +) -λ -π 0 g(r), ( 50
) p ′ (n, r, +) = λ + p(n, r, +) -λ -p(n -1, r, -) -λ + p(n -1, r, +), for 2 ≤ n ≤ k -+ 1, (51)
p ′ (n, r, -) = λ -p(n, r, -) -g(r)(p(n + 1, 0, +) + p(n + 1, 0, -)), for 1 ≤ n ≤ k -, (52)
p ′ (n, r, +) = λ + p(n, r, +) -λ + p(n -1, r, +), for k -+ 2 ≤ n ≤ k + , ( 53
) p ′ (n, r, -) = -g(r)(p(n + 1, 0, +) + p(n + 1, 0, -)), for k -+ 1 ≤ n ≤ k + -1, ( 54
) p ′ (k + + 1, r, +) = -λ + p(k + , r, +), ( 55
) p ′ (k + , r, -) = -g(r)p(k + + 1, 0, +). ( 56
)
By integrating Equations (51) and (52) for r from 0 to ∞ and summing up the two obtained equations, we deduce that p(n, 0, +) + p(n, 0, -) -λ + π n-1,+ -λ -π n-1,-is a constant for 2 ≤ n ≤ k -+ 1. With the same approach with Equations (53) and (54), we deduce that p(n, 0, +) + p(n, 0, -) -λ + π n-1,+ is equal to the same constant for k -+ 2 ≤ n ≤ k + . Using Equations ( 53) and (56), proves that this constant is also equal to p(k + + 1, 0, +) -λ + π k + ,+ . Finally, Equation (55), leads to p(k + + 1, 0, +) = λ + π k + ,+ . So, the aforementioned constant is 0. As a conclusion, we may write p(n, 0, +) + p(n, 0, -) = λ + π n-1,+ + λ -π n-1,-, for 2 ≤ n ≤ k -+ 1, p(n, 0, +) + p(n, 0, -) = λ + π n-1,+ , for k -+ 2 ≤ n ≤ k + , p(k + + 1, 0, +) = λ + π k + ,+ .

Combining this last set of equations with Equations ( 52) and (54), we get p(n, 0, -) = λ + π n,+ , for 1 ≤ n ≤ k -.

With the same approach as in Proposition 1 with Equation (52), we then deduce that π n,-= λ + λ - 8 Yet, we do not claim that the two-thresholds policy is optimal in this case. Expected throughput of served customers E(W + ), E(W -)

Notations

Expected waiting time of a customer who arrives after an arrival or a service E(R)

Expected remaining service time seen by an arriving customer at a non-empty system E(R + ), E(R -) Expected remaining service time seen by a customer who arrives after an arrival or a service Markov decision process

V + n (x), V - n (x)
Value function depending on the state of the system µ j Exponential rate of remaining service phase j in the Coxian distribution (1 ≤ j ≤ N ) r j Probability to enter remaining service phase j -1 after leaving remaining service phase j (1 ≤ j ≤ N ) in the Coxian distribution

Fig. 1

 1 Fig. 1 Impact of λ + and λ -on P (W > t) (x = 1)

Fig. 2

 2 Fig. 2 Impact of λ + and λ -on E(W ) (x = 1)

π 0 .

 0 Finally, Equation (28) leads to Equation (29). From Equation (12), we deduce that ∫ ∞ r=0 e -sr p ′ (n, r, -) dr =λ - ∫ ∞ r=0 e -sr p(n, r, -) dr -(p(n + 1, 0, +) + p(n + 1, 0, -)) ∫ ∞ r=0 e -sr g(r) dr, for n ≥ 1. Thus,

  a)) and (3(b)), we represent the performances measures as a function of k -

Fig. 3

 3 Fig. 3 Performance measures (λ+ = 1, λ -= 2, µ 1 = 5, µ 2 = 10, r 1 = 0, r 2 = 1, E(Q) = 0.6)

1 -( 2 Fig. 4

 124 Fig.[START_REF] Asmussen | Fitting phase-type distributions via the em algorithm[END_REF] Performance measures (λ + = 0.8, λ -= 0.1, q = 1/3, µ 1 = 0.5, µ 2 = 2, E(Q) = 0.3)

Proposition 6

 6 Under the condition, λ -≥ λ + , and(1-r x )µ x ≥ (1-r x+1 )µ x+1 , for 1 ≤ x ≤ N , the optimal

	admission policy has a threshold structure. More precisely, there exists two thresholds k + x and k -x such that,
	for k ≥ 0 and 1 ≤ x ≤ N ,
	-if an arrival occurs in state (kN + x, +), it is optimal to reject this customer if k > k + x , otherwise this
	customer is accepted,
	-if an arrival occurs in state (kN + x, -), it is optimal to reject this customer if k > k -x , otherwise this
	customer is accepted.

Table 1 Notations

 1 Exogenous parametersλ + , λ - Arrival rates after an arrival and after a service S Random variable which represents the service time duration x Expected service time cv Coefficient of variation of the service time distribution; it is the ratio of the standard deviation divided by the expected value g(.)Probability density function of the service time G * (.)Laplace-Stieltjes Transform (LST) of the service time;G * (s) = ∫ ∞ 0 g(t)e -st dt Decision parameters k + , k -Thresholds on the system size or on the remaining number of service phases to accept or reject customers at arrival after an arrival or a service p + , p - Probabilities that an arriving customer accepts to join after an arrival or a service Probabilities p t (n, r, +), p t (n, r, -) Probability-density of having n customers in the system, n ≥ 1 and a remaining service time of r, r ≥ 0, at time t after an arrival or a service p(n, r, +), p(n, r, +) p(n, r, +) = lim -) for n ≥ 1 π n,+ , π n,- Stationary probability to have n customers in the system at arbitrary instants after an arrival or a service (π n,+ = ∫ ∞ r=0 p(n, r, +) dr, π n,-= ∫ ∞ r=0 p(n, r, -) dr for n ≥ 1) π nStationary probability to have n customers in the system at arbitrary instants,π n = π n,+ + π n,-for n ≥ 0 p n,+ , p n,- Stationary probability to have n customers in the system at departure instants after an arrival or a service for n ≥ 0 p n Stationary probability to have n customers in the system at departure instants, p n = p n,+ + p n,-for n ≥ 0 r + n , r - Expected remaining service time seen by a customer who arrives after an arrival or a service with n customers present in the system, n ≥ 1. α n , β n Probability that n customers arrive during a service if the service is initiated by a service completion or by an arrival, n ≥ 0

		Performance measures
	E(Q d ), E(Q)	Expected number of customers in the system at departure and arbitrary instants
	E(Q)	Service level objective on E(Q)
	E(W )	Expected waiting time at arbitrary instants
	T s

t→∞ p t (n, r, +), p(n, r, -) = lim t→∞ p t (n, r, n

We refer the reader to the book of[START_REF] Hassin | To queue or not to queue: Equilibrium behavior in queueing systems[END_REF] for an overview on equilibrium behavior in queueing systems.

[START_REF] Altman | Non-threshold equilibrium for customers joining an M/G/1 queue[END_REF] shows an example where deterministic policies are not optimal for the admission control in an M/G/1 queue.

This particular Coxian distribution is an hypoexponential distribution.

if an arrival occurs after an arrival, then this customer is rejected if k > k + r , otherwise this customer is accepted,

if an arrival occurs after a service, then this customer is rejected if k > k - r , otherwise this customer is accepted.

The two-thresholds policy can be used to improve the classical one-threshold policy where the same decision is taken after an arrival or a service. By allowing different decisions after an arrival or after a service, a larger range of values may be reachable for E(Q). This may lead to a better solution for the optimization problem.

Optimal policy

We propose to formulate the routing problem as a Markov decision process (MDP) and next use the value iteration technique to prove the threshold structure of the optimal policy. We choose to approximate the service time duration by a Coxian distribution. Since Coxian distributions are dense in the field of all nonnegative distributions [START_REF] Schassberger | Warteschlangen[END_REF], the obtained results apply for a general distribution. Consider a Coxian random variable which represents the service duration. It is defined by the parameters µ j (µ j > 0, 1 ≤ j ≤ N ), and r j (r j ∈ [0, 1], 0 ≤ j ≤ N ) with r 1 = 0. The quantity r j is the probability to enter the remaining phase j -1 after leaving remaining phase j and the parameter µ j is the rate of the exponential distribution describing the random duration spent at remaining phase j. However, our system is not a standard Markov decision process (MDP). The form of the problem makes it a constrained MDP; maximize the throughput of served customers with a constraint on the expected number of customers in the system. Constrained MDP's can be solved using various techniques. Here we use one that introduces the constraint in the objective using a Lagrange multiplier. Under weak conditions it can be seen that the optimal policy for a certain Lagrange multiplier is optimal for the constrained problem if the value of the constraint under this policy attains exactly E(Q). From the theory on constrained MDP's it follows that this policy is stationary and randomizes in at most 1 state. For this and other results on constrained MDP's, see the book of [START_REF] Altman | Constrained Markov decision processes[END_REF]. The optimization problem may then be rewritten as min(E(Q) -c • T S ), where the coefficient c (c ≥ 0) is the Lagrange multiplier which translates the relative importance given, by the system manager, to the throughput of served customers (T S ) compared to expected number of customers in the system (E(Q)).

Let us denote by (x, y) a state of the system where x is the number of potential remaining phases of work for the server, x ≥ 0 and y is the nature of the last event; an arrival (y = +) or a service (y = -). We denote the transition rate from state (x, y) to state (x ′ , y ′ ) by q (x,y),(x ′ ,y ′ ) . Hence for x, x ′ ≥ 0 and y, y ′ ∈ {+, -}, we have

(1 -r j )µ j , if x = kN + j, x ′ = kN, y ′ = -for k ≥ 0, 1 ≤ j ≤ N, and y ∈ {+, -}, 0, otherwise, which corresponds to arrivals and service departures.

We choose to discretize our continuous-time model. This is possible because it is uniformizable (Section 11.5.2. in [START_REF] Puterman | Markov Decision Processes[END_REF]). We formulate a 2-step value function, in order to separate transitions and actions. We define the dynamic programming value functions W - n (x), W + n (x), V - n (x), and V + n (x) over n ≥ 0 steps, depending on the state of the system. We choose W -,+ 0 (x) = 0 and V -,+ 0 (x) = 0 for x ≥ 0. We assume without loss of generality that λ + + λ -+ ∑ N j=1 µ j = 1, such that the rate out of each state is equal to 1; thus we can consider the rates to be transition probabilities. We then may write for 1 ≤ x ≤ N , and k ≥ 0,

The operator W n represents the decision to accept or to reject a new customer from the system. After an arrival we have

, and after a service, we have

For each n > 0 and every state, there is a minimizing action at customer's arrival: accept this customer or reject this customer. One way of obtaining the long-run average optimal actions is to use the value iteration technique introduced by [START_REF] Bellman | Dynamic programming[END_REF] and [START_REF] Howard | Dynamic programming and Markov processes[END_REF], by recursively evaluating V n using Equation (44), for n ≥ 0. To prove the form of the optimal policy we want to establish structural properties of the value function. In particular, it would be interesting to obtain conditions for which threshold policies based on the number of customers in the system are optimal.

Future Research

Many questions following this study are open for future research. For instance, it would be interesting to include other features for the customer's behavior like abandonment or workload-dependency. We could also consider a multi-server queue instead of a single-server one. This however would not change the results since the observation of a change in the queue size would only occur when all agents are busy. Another extension of the model is the possibility of a customer's decision not only based on the last event but on a larger finite number of events. Finally, it could also be interesting to consider the symmetrical case where the server adapts its service rates to the last realized event.