domain: info math

text:
Raymond A Marie
email: raymond.marie@irisa.fr.

Performance Evaluation of Proactive Queues

Keywords: Performance evaluation, Analytical models, Stochastic modeling, Queueing models, proactive queue Evaluation de performance, Modèles analytiques, Modélisation stochastique, Files d'attente, Files proactives

Some service systems can be seen as proactive systems that anticipate future service demand. In this report, we present some basic proactive systems from the perspective of the performance criteria used for classical queues. The resources may anticipate services up to a given limit. This limit is a key factor with respect to the efficiency of the system under study. The performance variables are mainly service time, waiting time and response time; the metrics are mainly the probability distributions and the expectations of those, and also conditional distributions and expectations. According to the relative values of the limit of anticipation and of the number of servers, the results are more or less simple. We start this study with the case of the proactive single-server queue (M/M/1/<C>); the request arrival process is Poisson and the independent service times follow an exponential distribution; C denotes the limit of anticipation. Then, we study the proactive M/M/r/<C> queue; we have to consider two different situations, the easy case being when r is lower than (C + 1). Finally, we generalize the previous study by considering the case where the arrival rates and service rates depend on the state of the queue.

Introduction

Some service systems can be seen as proactive systems that anticipate future service demand. Perhaps the simplest paradigm is that of the pizza vendor who may decide to start production of one or more pizzas when he has no more orders pending. In a more professional context, let's consider an ad-hoc network of mini-sensors sending data in the form of messages transmitted through the network to an operating station. The energy is provided by a battery recharged by photovoltaic cells ( [START_REF] Doncel | Balancing energy consumption and losses with energy packet network models[END_REF]). The transmission of a message requires a quantum of energy and the fully charged battery contains a limited number of quanta. In case of strong activity of the sensor or in case of bad weather conditions, the battery discharges and it is necessary to wait to recover a quantum of energy to send the possible pending messages. In this example, the server is the photovoltaic unit that provides the quanta and the users are the messages to be transmitted.

In an industrial production context, a cell included in a Kanban chain can also be seen as providing a proactive service. The cell under consideration can contain several machines in parallel, implemented to satisfy the demand as well as possible, without exceeding a maximum level of anticipation C. The machines are the resources (the servers) and the clients are the requests of the cell located downstream. By persevering in this way, we can include storage systems; thus, even if he does not manufacture anything, the car dealer anticipates future demand by ordering cars whose future owners he does not yet know.

In the context of integrated logistical support (ILS), let us consider the example of a maintenance workshop where a sub-assembly of a given piece of equipment is repaired, for example the engine of an aircraft. Two policies (noted respectively A and B) are possible. With policy A, the engine taken in hand by a repair person will be reassembled on the original aircraft. With policy B, the engine is disconnected from the aircraft and the stripped aircrafts are placed in a discipline queue (first in, first out) and the first engine to be repaired is reassembled on the aircraft at the head of the queue. When technically practicable, this policy B makes sense if the maintenance center is initially allocated a stock of C engines in good condition. Note that these two policies do not provide the same performance in terms of aircraft unavailability. With policy A, the engine repair time corresponds to the service time of the customer (here the aircraft). If there are r repairers (or r repair teams), the whole set behaves like a classical M/M/r queue. The aircraft downtime is the sum of a possible waiting time for lack of available repairers and the engine repair time. With policy B, the set is proactive and the service seen by the customers (the airplanes) can be very short, or even almost zero if the time for the engine to be disengaged from the airplane is negligible.

The purpose of this research report is to present some basic proactive systems from the perspective of the performance criteria used for classical queues. We start by studying the case of the single-server queue (next Section); the request arrival process is Poissonian and the independent service times follow an exponential distribution. In Section 3, we study the proactive M/M/r queue. Then, in the section 4, we generalize the study of the previous section to the queue by considering the case where the arrival rates and service rates depend on the state of the queue. Finally, Section 5 expose our conclusions.

2 Single-server Queue M/M/1/<C>

Modelization

When the client queue is empty, the proactive server can produce up to C items/services. C is the value of the allowed anticipation window. The diagram in Figure 1 illustrates the operation of the queue. The queue has an operator of type "Join". The arrivals of the requests follow a Poissonian process of rate λ and the service times are independent and identically distributed (iid) according to an exponential distribution of rate µ. The parameter C corresponds to the maximum number of anticipated services. The variables i and j denote respectively the number of services waiting for a customer and the number of customers waiting for a service. Note that, by construction, the product of the values taken by the variables i and j is zero for any state of the queue. Both variables are simultaneously zero when the server anticipates a first service. If a client shows up before the end of this anticipated service, the client must wait until the end of the service. If j clients are present in the queue, they are served according to the discipline "first come -first served" (paps). These assumptions will allow us to consider a Markov model of the queue.

The variable k defined by the equality

k △ = (C -i + j) , (1)
makes it possible to determine in an ordered form the possible states of this waiting system. The values of k less than C corresponding to the existence of anticipated services and the values of k greater than C corresponding to customers waiting (or in service for the lead customer). The variable k equals zero when all the anticipable services are available (i = C). It is equal to C when both no customer is present (j = 0) and, on the other hand, no advance service is yet available (i = 0). For a given value k, the associated values i and j are characterized by the relations

i = (C -k) + (2)
j = (k -C) + (3)
where notation (x) + means max(0, x).

Inria

Performances Determination

Let us denote p k the asymptotic probability of state k, k = 0, 1, 2, ... . First we determine these asymptotic state probabilities (cf. Figure 2). It is easy to check, after having put ρ = λ/µ, that the steady state state probabilities are written p k = (1 -ρ)ρ k , k = 0, 1, 2, ... under the necessary ergodic condition ρ < 1.

In order to determine the expectation of the response time, let's start by determining the expectation of the number of requests present in the system:

E[N C ] = ∞ j=1 jp C+j = ∞ j=1 j(1 -ρ)ρ C+j , = (1 -ρ)ρ C+1 ∞ j=1 jρ j-1 = (1 -ρ)ρ C+1   ∞ j=1 ρ j   ′ , = (1 -ρ)ρ C+1 1 (1 -ρ) ′ = (1 -ρ)ρ C+1 1 (1 -ρ) 2 , = ρ C+1 (1 -ρ) .
Remembering that for a file M/M/1, the expectation E[N ] is equal to ρ/(1 -ρ), we notice that this last expectation is related to E[N C ] via the relation:

E[N C ] = ρ C E[N C ]| C=0 = ρ C E[N ] .
Let's use Little's formula to deduce the expectation of the response time:

E[R C ] = E[N C ] λ = 1 µ ρ C (1 -ρ) = ρ C E[R C ]| C=0 . (4)
With this queue, a customer will receive an immediate service if, upon arrival, i is strictly positive. If i and j are both zero, it will be enough for him to wait until the end of the (previously anticipated) service to exit the system. If the value of the variable j is strictly positive when it arrives, the client will have to wait for the server to finish serving the customer (s) already present before getting the server for himself. The probability that the wait for server availability is strictly positive is the probability that the client finds the system in a k state strictly greater than C; By complementarity, if W C denotes the waiting time to obtain server availability, we have:

P(W C = 0) = C j=0 p j , = (1 -ρ) C j=0 ρ j = (1 -ρ) 1 -ρ C+1 (1 -ρ) , = (1 -ρ C+1 ) . (5)
The probability P(S C = 0) that the response time is itself zero must be less than the probability P(W C = 0) since the event [S C = 0] is included in the event [W C = 0]. The sought probability is that the client finds the system in a state k such that the variable i is strictly positive:

P(S C = 0) = C-1 j=0 p j , = (1 -ρ) C-1 j=0 ρ j = (1 -ρ) 1 -ρ C (1 -ρ) , = (1 -ρ C ) . (6)
In Figure 3, we have represented the probabilities that a service or waiting time is zero as a function of C, for C = 0. 1, 2, ..., and 5, when ρ = 0.75 and µ = 0.5. These two probabilities are increasing according to the C level of anticipation. It can be verified, as the previous results have shown, that P(S C+1 = 0) = P(W C = 0). Given that the service time is exponentially distrbuted with rate µ, the expectation of the service time perceived by the customers is given by

E[S C ] = 1 µ P(S C > 0) = 1 µ (1 -(1 -ρ C )) = ρ C µ = ρ C E[S C ]| C=0 . (7)
Since the response time remains equal to the sum of the waiting time for the server (possibly zero time) and the service time perceived by the client (possibly zero time), we can determine the expectation of waiting

E[W C ] subtracting E[S C ] to E[R C ]: E[W C ] = E[R C ] -E[S C ] , = ρ C µ(1 -ρ) - ρ C µ = ρ C µ ρ (1 -ρ) -1 ,
= ρ C ρ µ(1 -ρ) , = ρ C E[W C ]| C=0 . (8)
In this single-server queue, the ρ parameter is necessarily less than 1. This causes ρ C to decrease as the C capacity increases. Thus, the expectations

E[N C ], E[R C ], E[W C
] and E[S C ] decreases according to the same proportion when increasing the C capacity.

In Figure 4, we have represented the evolution of the expectations E[W C ] and E[R C ] as a function of C , for C = 0. 1, 2, ..., and 5, when ρ = 0.75 and µ = 0.5. As we could expect, these two expectations are decreasing according to the C level of anticipation. For a given value C, the difference (

E[R C ] -E[W C ]
) is equal to the expectation of the service time E[S C ] that decreases when C increases (cf. the result 7) ; for C = 0, the difference is equal to 1/µ = 2. It can also be verified, as the previous results have shown, that

E[R C+1 ] = E[W C ].
Let's now look at the conditional expectation of waiting time knowing that the client did not find the server available when he arrived. For that, let's start from the following relation:

E[W C ] = E[W C |W C > 0]P(W C > 0) + E[W C |W C = 0]P(W C = 0) , which, since the conditional expectation E[W C |W C = 0] is null, allows to write: E[W C |W C > 0] = E[W C ] P(W C > 0)
. In the case of the M/M/1 queue, it is not difficult to find that this conditional expectation is

E[W |W > 0] = 1/µ(1 -ρ).
For the queue with a proactive server, we get:

E[W C |W C > 0] = E[W C ] P(W C > 0) , = ρ C ρ µ(1 -ρ) 1 ρ C+1 , = 1 µ(1 -ρ) , = E[W C |W C > 0]| C=0 = E[W |W > 0] .
We observe that this conditional expectation is identical to that of the classical M/M/1 queue. This result is not surprising because, if the client finds one or more clients when they arrive, everything happens as if he was in a non-proactive M/M/1 queue.

In the same way, since the probability P(R C = 0) must be equal to the probability P(S C = 0), we obtain from expressions 4 and 6 the expression of conditional expectation

E[R C |R C > 0]: E[R C |R C > 0] = E[R C ] P(R C > 0) = E[R C ] P(S C > 0) , = ρ C µ(1 -ρ) 1 ρ C = 1 µ(1 -ρ) , = E[R C |R C > 0]| C=0 = E[R C ]| C=0 , Inria = E[R] = E[W C |W C > 0] .
Note that this conditional expectation is identical to that of the classical M/M/1 queue for which the response time is positive with probability 1. It is also equal to the conditional expectation of the waiting time, for any value of C. Recall that, for the classical M/M/1 queue, the unconditional expectation E[W ] is related to the conditional expectation of the waiting time by the relation

E[W ] = ρ E[W |W > 0].
Let's go further and look for the probability distribution of W C . We already know that it is a mixed distribution with a mass at zero, equal to (1 -ρ C+1 ). We will consider any customer but to clarify the situation, we will give a name to this customer chosen purely randomly: August. Since arrivals follow a Poisson process, the probability that August finds j clients in front of him on arrival is equal, when j is strictly positive, to (1 -ρ)ρ C+j . Indeed, Poisson arrivals find the system in equilibrium (this is the property pasta : Poisson Arrivals See Time Averages). Knowing that he finds j clients in front of him, the waiting time experimented by August is distributed according to the Erlang distribution of order j and of parameter µ whose density is written:

f X (t) = µ j t j-1 (j -1)! e -µt , t ≥ 0 .
Which makes it possible to write the unconditional distribution:

dF W C (t) = (1 -ρ C+1 )δ(t) + ∞ j=1 (1 -ρ)ρ C+j µ j t j-1 (j -1)! e -µt , = (1 -ρ C+1 )δ(t) + (1 -ρ)ρ C e -µt ∞ j=1 λ j t j-1 (j -1)! , = (1 -ρ C+1 )δ(t) + (1 -ρ)ρ C e -µt λ ∞ i=0 (λt) i i ! , = (1 -ρ C+1 )δ(t) + (µ -λ) µ ρ C λe -(µ-λ)t , = (1 -ρ C+1 )δ(t) + ρ C+1 (µ -λ)e -(µ-λ)t .
where δ(t) denotes the Dirac function. Thus, the probability distribution of W C is a mixed probability distribution with a mass at zero, equal to (1-ρ C+1 ), and a continuous part over the domain ]0, +∞] corresponding to the exponential with rate (µ -λ) weighted by the coefficient ρ C+1 .

Note that if C = 0, we find the result known for the non-proactive M/M/1 queue. The conditional probability distribution of the waiting time knowing that the client has waited is therefore the exponential distribution with rate (µ -λ) as well for the proactive queue as for the M/M/1 standard queue.

the cumulative probability distribution of W C can be written:

F W C (t) = (1 -ρ C+1 )u(t) + ρ C+1 (1 -e -(µ-λ)t ) . (9)
In the same way, we can determine the probability distribution of the response time R C . This mixed probability distribution also has a mass in zero, equal to (1 -ρ C ), and the probability that August finds j clients in front of him on his arrival remains equal to (1 -ρ)ρ C+j . Given its own service, the response time is distributed according to the Erlang distribution of order (j + 1) and of parameter µ. Note that if j = 0, the response time corresponds to the end of the service started by the server before the August arrival time; but the service time being distributed according to the exponential distribution of rate µ, it is the same for this residual service time. We thus have:

dF R C (t) = (1 -ρ C )δ(t) + ∞ j=0 (1 -ρ)ρ C+j µ j+1 t j j! e -µt , = (1 -ρ C )δ(t) + ∞ i=1 (1 -ρ)ρ C-1+i µ i t i-1 (i -1)! e -µt , = (1 -ρ C )δ(t) + (1 -ρ)ρ C-1 e -µt ∞ i=1 λ i t i-1 (i -1)! = (1 -ρ C )δ(t) + (1 -ρ)ρ C-1 e -µt λ ∞ k=0 (λt) k k ! = (1 -ρ C )δ(t) + (µ -λ) µ ρ C-1 λe -(µ-λ)t = (1 -ρ C )δ(t) + ρ C (µ -λ)e -(µ-λ)t ,
Again, the continuous part of this mixed probability distribution corresponds to the exponential distribution of rate (µ -λ), weighted here by the factor ρ C . Let's note that by incrementing the allowed anticipation window from C to (C + 1), the probability distribution of the response time R C+1 is that of the waiting time W C . Thus the cumulative distribution function of R C can be written:

F R C (t) = (1 -ρ C )u(t) + ρ C (1 -e -(µ-λ)t ) . (10)
This cumulative distribution function is represented in Figure 5, as a function of C, C= 0, 1, 2, 3, and 4; with the same parameter values as in the previous figure: ρ = 0.75 and µ = 0.5. The complementary cumulative distribution function FR C (t) = P(R C > t) is equal to the decreasing function ρ C e -(µ-λ)t . For an arbitrary value t, we have FR C+1 (t) = ρ FR C (t) and as ρ is necessarily less than one, the gap between the curves is reduced when C increases. It follows from the relations 9 and 10 that F W C (t) = F R C+1 (t). Thus, Figure 5 also represents the distribution function F W C (t) for values C = 0, 1, 2 and 3 ; This explains why we did not previously use a figure in order to illustrate the distribution function F W C (t).

3 Multi-server Queue M/M/r/<C>

In the presence of several servers, the study of the queue becomes more complex but we will see that it is still possible to determine many performance metrics. The queue that we characterize by the quadruplet M/M/r/< C > where r denotes the number of servers and the last element still denotes the maximum number of anticipatable services, corresponds to a Markovian model and will be easier to study if r is not greater than (C + 1).

Modelization

We model this new queue by defining the variables i, j and k in a similar way to the previous case (see Figure 1 (page 12) as well as the relation 1) (page 4)). As there are several servers, it is necessary to specify the service procedure when the variable i is strictly positive, so that the maximum number of anticipated services does not exceed the limit C. We can realize the procedure thanks to C tokens representing as many authorizations to anticipate a service, the token being released when a client who presents himself can benefit from the anticipated service. Note that such a procedure is easily modeled using a timed stochastic Petri net. For a given value k, the associated values i and j are characterized by the same relations as in the mono-server case:

i = (C -k) + (11) j = (k -C) + (12)
where again notation (x) + means max(0, x). The diagram in Figure 6 illustrates the operation of the queue in this multiserver case where the number r of servers may be less than or greater than C. In order to further explore this alternative, let us note n a the number of active servers (busy with a proactive task or not) and n pa the number of proactive servers. Figure 7 illustrates, in the case where C is greater than or equal to (r-1), the values of variables i, j, n a and n pa based on k. The dashed curve represents the number n a of active servers; the number of proactive servers n pa corresponds to the curve in small dots. Figure 8 illustrates the equivalent results in the case where C is less than (r -1). In both cases, the number n a of active servers remains equal to min(k, r). The general expression of the number of servers occupied with a proactive task, more complex, can be written:

n pa = 1 (k≤max(r,C)) inf(k, min(r, C)) + 1 (max(r,C)<k<(r+C)) (r + C -k) ,
knowing that a server is no longer occupied with a proactive task if at least r clients are present (The task may have changed from the "proactive " status to the "active "). Note that if C is greater than or equal to (r -1), all the servers are active as soon as j is positive, i.e. also as soon as k is greater than C. Let's now consider the j clients that, at an instant, are in the buffer. These customers are positioned according to their order of arrival : the customer located in p-th position arrived after (resp. before) the customer located in (p -1)-th (resp. (p + 1)-th) position. When j is greater than r, the customer located in p-th position is on hold if p is greater than r. If p is less than or equal to r, the query of this client is already in progress. When r is greater than 1, we need to specify the intended service protocol: servers perform identical requests that satisfy the client at the top of the buffer. The probability distribution of the service of a client will depend on the state of the file seen by this client but also on relative values of the parameters r and C.

Our model is again a ctmc associated with the possible values of k, k = 0, 1, 2, .... Its transition graph is first represented in Figure 9 in the case where the anticipation C is greater than or equal to (r -1) and then represented in Figure 10 in the case where the anticipation C is lower than (r -1). Because of these ordered states, these particular ctmc are also called birth-death-processes (bdp). If C ≥ (r -1) (condition that can also be written r ≤ (C + 1)), the r servers are active as soon as k ≥ r and in particular the r servers are active when one or more clients are present. On the other hand, if C < (r -1) (a condition that can also be written r > (C + 1)), the r servers are all active only when k is greater than or equal to r, i.e. when the number j of present clients is greater than or equal to (r -C) + , in order to avoid violating the constraint on the maximum number of anticipated services.

0 1 C-1 C C+1 r r+1 • • • • • • • • • µ λ Cµ λ (

RR n 9496

Note that in the case where C is greater than or equal to (r -1), the r servers are all active when j is strictly positive, unlike when C is less than (r -1). As we will see later, some performance metrics will have a simpler expression in the first case than in the second.

Performance Determination

Since now r is eventually greater than 1, let redefine ρ as ρ = λ rµ . So, ρ represent the loading rate of this proactive multi-serveur queue. It is easy to verify that the state probabilities of the stationary regime satisfy the relations:

p k = ρ r k p k-1 if k ≤ r , ρp k-1 if k ≥ r .
This allows us to express all probabilities in terms of the probability p 0 :

p k =      ρ k r k k ! p 0 if k ≤ r , ρ k r r r ! p 0 if k ≥ r .
Subject to satisfying the necessary and sufficient condition of ergodicity ρ < 1, the probability p 0 is obtained through the normalizing equation:

p 0 = r-1 i=0 ρ i r i i! + r r r! ρ r (1 -ρ) -1 , (13)
Note that these state probabilities p k , k = 0, 1, ..., do not depend on C. Of course, this is not the case for the performances we are looking for.

Case

C ≥ (r -1)
First, we consider the situation where C equals or exceeds (r -1) ; this is the simplest situation. Let's place ourselves at the moment of arrival of an arbitrary client (to whom we have already given the nice name of August in order to distinguish him), finding the queue in the state (i, j). Several cases are possible:

-case i ≥ 1 (and therefore j = 0): the service is instantaneous (the service distribution is the Dirac's distribution δ(0)),

-case i = 0 and j = 0: August takes the leading position in the buffer and his service time follows an exponential distribution of rate rµ.

-case 1 ≤ j < r (and i = 0): August takes the position (j + 1) and his service time follows an Erlang distribution of order (j + 1) and of parameter rµ.

-case j ≥ r (and i = 0): August takes the position (j + 1) and starts his wait which ends when he reaches the position r; his service time will follow an Erlang distribution of order r and parameter rµ.

Inria

Let us first determine the expectation of the number of requests waiting in the queue:

E[ν C ] = ∞ j=r+1 (j -r)p C+j , = ∞ k=C+r+1 (k -C -r)p k , = ∞ k=C+r+1 (k -C -r) ρ k r r r ! p 0 = p 0 r r r ! ∞ i=1 iρ i+C+r , = ρ C+r+1 p 0 r r r ! ∞ i=1 iρ i-1 = ρ C+r+1 p 0 r r r ! 1 (1 -ρ) 2 , = ρ C p 0 r r r ! ρ r+1 (1 -ρ) 2 , ( 14
) = ρ C E[ν C ]| C=0 . (15)
Note that, given the initial expression of E[ν C ], this result is valid regardless of the value of C (that C is less than, equal to or greater than (r -1)). Thanks to Little's formula, we immediately deduce from this result the expectation of the waiting time:

E[W C ] = E[ν C ] λ , = ρ C p 0 r r r ! ρ r rµ(1 -ρ) 2 , = ρ C E[W C ]| C=0 , (16)
a result also valid regardless of the value of C. Let's now determine the expectation of the number of requests present in the system (in the case where C ≥ (r -1)):

E[N C ] = ∞ j=1 jp C+j , = ∞ k=C+1 (k -C)p k = ∞ k=C+1 (k -C) ρ k r r r ! p 0 , = p 0 r r r ! ∞ i=1 iρ i+C = ρ C+1 p 0 r r r ! ∞ i=1 iρ i-1 , = ρ C+1 p 0 r r r ! 1 (1 -ρ) 2 . ( 17
)
Note that this expression is simpler than the one obtained for the non-proactive queue. It will be quite the opposite in the second case where C will be less than (r -1)... unfortunately !

We deduce the expectation of the response time using Little's formula:

E[R C ] = E[IN C ] λ , = ρ C p 0 r r r ! 1 rµ(1 -ρ) 2 . ( 18
)
Comparing expressions 14 and 17, we notice that, when

C ≥ r, E[N C ] can express itself according to E[ν C ]: E[N C ] = E[ν C-r ] C ≥ r . ( 19
)
From this, it is immediately deduced that

E[R C ] can express itself according to E[W C ]: E[R C ] = E[W C-r ] C ≥ r . (20)
In the boundary case where

C = r -1, we obtain E[N r-1 ] = ρ -1 E[ν 0 ] and E[R r-1 ] = ρ -1 E[W 0 ]. Let's now determine the expectation E[N S C
] of the number of in-service requests using results 14 and 17:

E[N S C ] = E[N C ] -E[ν C ] , = ρ C p 0 r r r ! 1 (1 -ρ) 2 ρ -ρ r+1 , = ρ C+1 p 0 r r r ! (1 -ρ r ) (1 -ρ) 2 . ( 21
)
Let's use Little's formula to obtain the expectation of service time:

E[S C ] = E[N S C ] λ , = 1 r µ ρ C p 0 r r r ! (1 -ρ r ) (1 -ρ) 2 . ( 22
)
It should be noted that this unconditional expectation of service time takes into account customers whose service time is zero. Let us recall here that these last two results giving the expressions of E[N S C ] and E[S C ] are valid only if C is equal to or greater than (r -1) (and thus not valid when

C = 0 if r > 1!).
Let us determine the probability P(S C = 0) that the service time of a customer is zero. This probability is the one that the client finds the system in a state k such that the variable i is strictly positive.

P(S C = 0) = C-1 j=0 p j , = ∞ j=0 p j - ∞ j=C p j = 1 -p 0 r r r! ∞ i=C ρ i , Inria = 1 -p 0 ρ C r r r! ∞ j=0 ρ j , = 1 -ρ C r r r! p 0 (1 -ρ) . ( 23
)
Note that this result is also valid if there is only one server since the probability p 0 is then equal to (1 -ρ) (cf the expression 6). By associating this last result 23 with the expression of E[N C ] (cf. 17), it is easy to show that:

E[N C ] = ρ (1 -ρ) P(S C > 0) . ( 24
)
This expression for the expectation E[N C ] is worth commenting on. First note that the probability P(S C > 0) represents the percentage of customers who cannot pass through the station without stopping.

If Z(t) denotes the state of the cmtc at time t (the one shown in Figure 9), let us respectively note B(t) and I(t) the two processes defined using the Z(t) process as:

B(t) = 1 if Z(t) > C , 0 else . I(t) = 1 if Z(t) < C , 0 else .
Thus, B(t) is 1 if at least one client is present in the queue (j > 0 and i = 0) and I(t) is 1 if at least one anticipated service is available in the queue (i > 0 and j = 0).

Figure 11 represents possible partial trajectories over a time interval between t 0 and at t f . At time t 0 , we see that B(t 0 ) = 1 and I(t 0 ) = 0; that means that at least one client is present in the queue.

I(t) B(t) t 0 t 1 t 2 t 3 t 4 t 5 t 6 t 7 t 8 t 9 t 10 t 11 t 12 t 13 t 14 t f U 1 D 1 U 2 D 2 U 3 D 3
Figure 11: Partial trajectories of B(t) and of I(t) over a time interval between t 0 and t f .

We see that there are 14 changes of state that we can interpret in the following way:

-At time t 1 , the service of the last customer in the queue ends: B(t + 1 ) = 0 and I(t 1 ) = 0; the next event will be either a new arrival or an early service provisioning.

-At time t 2 , an anticipated service occurs: B(t 2 ) = 0 and I(t + 2 ) = 1. -At time t 3 , a customer arrival ends a possible series of instantaneous services performed on the time interval of lenght U 1 = (t 3 -t 2 ): B(t 3 ) = 0 and I(t + 3 ) = 0; the next event can be either a new arrival or an anticipated service provisioning.

-At time t 4 , a customer arrival occurs (before the provision of an anticipated service):

B(t + 4 ) = 1 and I(t 4 ) = 0.
-At time t 5 , the last of the customers who arrived on the interval [t 4 , t 5 [ leaves the queue, his service terminated: B(t + 5 ) = 0 and I(t 5 ) = 0; the next event can be either a new arrival or a provision of an anticipated service.

-At time t 6 , a customer arrival occurs (before the provision of an anticipated service):

B(t + 6 ) = 1 and I(t 6 ) = 0.
-At time t 7 , the service of the last customer in the queue ends: B(t + 7 ) = 0 and I(t 7 ) = 0; the next event can be either a new arrival or an early service provision.

-At time t 8 , an early service provisioning occurs: B(t 8 ) = 0 and I(t + 8 ) = 1, and so on... As can be seen in the figure, time can be decomposed into a sequence of intervals U i and D i , i = 1, 2, 3, .... During each interval U i , i = 1, 2, 3, ..., the queue is empty because the services are instantaneous.

The state k = C of the cmtc is a regenerative point for the process Z(t) (and thus for the process N (t), defined as the number of clients present in the queue (i.e. N (t) = max(0, (Z(t) -C))) as well as for the processes B(t) and I(t). This means that the random variables U i , i = 1, 2, 3, ..., are iid. The same is true for the random variables D i , i = 1, 2, 3, .... Moreover, the latter are independent of the former.

During each interval U i , i = 1, 2, 3, ..., the client queue is empty because the services are instantaneous. While during each interval D i , i = 1, 2, 3, ..., there is no anticipated service and the queue behaves like a M/M/1 queue with arrival rate rλ and service rate rµ; , i.e. with a load rate equal to ρ (because C ≥ (r -1)). Considering the queue in stationary regime, this allows us to write the expressions of the following conditonal expectations:

E[N C |I(t) = 1] = 0 , (25)
E[N C |I(t) = 0] = ρ (1 -ρ) , ( 26
)
where t is any instant of the queue in steady state. By deconditioning, we retrieve the relation 24:

E[N C ] = 0 + ρ (1 -ρ) × E[D i ] E[U i ] + E[D i ] , = ρ (1 -ρ) P(S C > 0) . (27)
We can now determine the conditional expectation of service time knowing that it is not zero:

Inria E[S C |S C > 0] = E[S C ] P(S C > 0) , = 1 r µ ρ C p 0 r r r ! (1 -ρ r ) (1 -ρ) 2 r! ρ C r r (1 -ρ) p 0 , = (1 -ρ r ) rµ(1 -ρ) . ( 28
)
We can find this result without knowing the unconditional expectation E[S C ] noticing that, knowing that the newcomer finds j clients present, his service time expectation will be equal to (j + 1) r µ for j = 0, 1, ...(r -2) and to r r µ for j ≥ (r -1):

P(S C > 0)E[S C |S C > 0] = r-2 j=0 (j + 1) r µ p C+j + ∞ j=r-1 r r µ p C+j , = r-1 j=1 j r µ p C+j-1 + r r µ ∞ i=0 p C+r-1+i , = p 0 ρ C-1 r µ r r r!    r-1 j=1 jρ j + rρ r ∞ i=0 ρ j    , = p 0 ρ C-1 r µ r r r! ((r -1)ρ (r+1) -rρ r + ρ) (1 -ρ) 2 + rρ r (1 -ρ) , = p 0 ρ C-1 r µ r r r! ρ(1 -ρ r ) (1 -ρ) 2 , = p 0 ρ C r r r! (1 -ρ r ) rµ(1 -ρ) 2 .
Which, by dividing by the probability P(S C > 0) gives us back the expression 28. The wait of a customer is null if, on his arrival, the service which he will benefit is already started (possibly already available). The probability that a customer's wait is zero is then expressed as:

P(W C = 0) = C-1+r j=0 p j , = ∞ j=0 p j - ∞ j=C+r p j = 1 -p 0 r r r! ∞ i=C+r ρ i , = 1 -p 0 ρ (C+r) r r r! ∞ j=0 ρ j , = 1 -ρ (C+r) r r r! p 0 (1 -ρ) . ( 29
)
Note that this formula is valid when r = 1 since the probability p 0 is then equal to (1 -ρ) (cf. expression 5). This result makes it possible to write the complementary probability:

P(W C > 0) = ρ (C+r) r r r! p 0 (1 -ρ) . ( 30
)
For the non-proactive M/M/r queue, it's easy to show that:

P(W > 0) = ρ r r r r! p 0 (1 -ρ) . ( 31
)
This allows to write a relation between these probabilities in the form:

P(W C > 0) = ρ C × P(W C > 0)| C=0 = ρ C P(W > 0) , (32)
as well as:

P(W C+1 > 0) = ρ P(W C > 0) . (33)
Note also, thanks to the relations (30) and ( 23), that

P(W C > 0) = ρ r P(S C > 0) = P(S C+r > 0) . (34)
We can now determine the conditional expectation of waiting knowing that it occurs:

E[W C |W C > 0] = E[W C ] P(W C > 0) , = ρ C p 0 r r r ! ρ r rµ(1 -ρ) 2 r! ρ (C+r) r r (1 -ρ) p 0 , = 1 rµ(1 -ρ) . (35)
As in the case where the server is unique, we verify here that this conditional expectation of waiting for the proactive queue is identical to that of the non-proactive queue.

Recall here that if the probability p 0 does not depend on the variable C, it depends on the variable r. It is therefore not a good idea to compare the results obtained by varying the number of servers without taking it into account.

Let's now look for the probability distribution of W C , always when C ≥ (r -1). We already know that W C is a mixed random variable with a mass at zero, equal to P(W C = 0) (given by the expression 29). Since arrivals follow a Poisson process, the probability that August finds j clients in front of him on arrival is equal, when j is strictly positive, to p (C+j) (property pasta). August suffers a wait if j is greater than (r -1). In this case, its waiting time is distributed according to the Erlang probability distribution of order (j -r + 1) and of parameter rµ whose density is written as:

f X (t) = (rµ) (j-r+1) t j-r (j -r)! e -rµt , t ≥ 0 .
Inria What, through a deconditioning, allows to write:

dF W C (t) = P(W C = 0)δ(t) + ∞ j=r p C+j (rµ) (j-r+1) t j-r (j -r)! e -rµt , = P(W C = 0)δ(t) + p 0 ρ C+r r r r! rµe -rµt ∞ j=r ρ j-r (rµ) (j-r) t j-r (j -r)! , = P(W C = 0)δ(t) + p 0 ρ C+r r r r! rµe -rµt ∞ j=r λ (j-r) t j-r (j -r)! , = P(W C = 0)δ(t) + p 0 ρ C+r r r r! rµe -rµt ∞ i=0 (λt) i i ! , = P(W C = 0)δ(t) + p 0 ρ C+r r r r! rµe -rµt e λt , = P(W C = 0)δ(t) + p 0 (1 -ρ) ρ C+r r r r! (1 -ρ)rµe -(rµ-λ)t , = P(W C = 0)δ(t) + p 0 (1 -ρ) ρ C+r r r r! (rµ -λ)e -(rµ-λ)t , = P(W C = 0)δ(t) + (1 -P(W C = 0))(rµ -λ)e -(rµ-λ)t . (36)
where δ(t) denotes the Dirac's function. The probability distribution of W C is a mixed one with a mass in zero, equal to the probability P(W C = 0) whose expression is given by the relation 29, and a continuous part on the domain ]0, +∞] corresponding to the exponential distribution of rate (rµ -λ) weighted by the factor (1 -P(W C = 0)).

Let's now look for the probability distribution of R C , again when C ≥ (r -1). We already know that it is a mixed distribution with a mass in zero, equal to P(S C = 0) (given by the expression 23). If on arrival, August finds less than r clients present while i = 0, he is taken care of immediately and his service follows an Erlang distribution of order (j + 1) (j being the number of clients present at the time of its arrival) and with parameter rµ because all the servers are active (we are in the case C ≥ (r -1)). If j ≥ r, August must wait until he finds himself in the r th position; his waiting time follows the Erlang distribution of order (j -r + 1) and of parameter rµ. Then the client goes on his service following the Erlang distribution of r order r and with parameter rµ. Waiting times and service being independent, the probability distribution of the conditional response time knowing j (and i = 0) is the Erlang distribution of order (j + 1) and of parameter rµ, result identical to the previous case where j was less than r.

So we have:

dF R C (t) = P(R C = 0)δ(t) + ∞ j=0 p C+j (rµ) (j+1) t j (j ! e -rµt , = P(R C = 0)δ(t) + p 0 ρ C r r r! rµe -rµt ∞ j=0 ρ j (rµ) j t j j ! , = P(R C = 0)δ(t) + p 0 ρ C r r r! rµe -rµt ∞ j=0 (λt) j j ! , = P(R C = 0)δ(t) + p 0 ρ C r r r! rµe -rµt e λt , = P(R C = 0)δ(t) + p 0 (1 -ρ) ρ C r r r! (1 -ρ)rµe -(rµ-λ)t , = P(R C = 0)δ(t) + (1 -P(R C = 0))(rµ -λ)e -(rµ-λ)t , (37)
Again, the continuous part of this mixed probability distribution corresponds to the exponential distribution of rate (rµ -λ), weighted here by the probability P(R C > 0); the latter being equal to the probability P(S C > 0) that the service is not zero.

Cas C < (r -1)

Let's now examine the situation when r is greater than (C + 1). This situation is more complex than the previous one because now the number n a of activated servers is limited to the quantity n a = min((C + j), r)s where j corresponds to the number of clients present in the buffer. Indeed, one should not anticipate more than C services. However, some previous results remain valid when r is greater than (C+1). This is the case of expectations E[ν C ] (cf. ( 14) and ( 15) 16)). But this is not the case for the expectation E[N C ] whose expression will become more complex. This is because even though this last expectation can always be written as when C was greater than or equal to (r + 1), in the form:

) et E[W C ] (cf. (
E[N C ] = ∞ k=C+1 (k -C)p k ,
the expressions of the first values of p k as a function of p 0 do not respect the general form obtained when k becomes greater than (r -1).

Here, considering what we already know the expression of E[ν C ], it seemed preferable to calculate first the expectation E[N S C ] of the number of requests during service (in the present case where C < (r -1)):

We show in appendix A that E[N S C ] can be written as:

E[N S C ] = p 0 r-C i=1 i (rρ) i+C (i + C)! + p 0 r r r ! ρ r+1 (1 -ρ C ) (1 -ρ) 2 + (r -C) (1 -ρ) . ( 38
)
Note that for C = 0, this expression 38 gives the expectation of the number of busy servers in queue M/M/r (i.e., the value rρ):

E[N S C ]| C=0 = p 0 r i=1 i (rρ) i (i)! + rp 0 r r r ! ρ r+1 (1 -ρ) , Inria = rρ p 0   r-1 j=0 (rρ) j j ! + r r r ! ρ r (1 -ρ)   , = rρ p 0 [p 0 ] -1 , = rρ . ( 39
)
We can still rewrite the right part of the expression 38 in the following form,

E[N S C ] = p 0 r-C i=1 i (rρ) i+C (i + C)! + p 0 r r r ! ρ r+1 1 (1 -ρ) 2 + (r -C) (1 -ρ) -p 0 r r r ! ρ C+r+1 (1 -ρ) 2 = p 0 r-C i=1 i (rρ) i+C (i + C)! + E[ν C ]| C=0 [1 + (r -C)(1 -ρ)] -E[ν C ] (40)
Thus we get an expression of the expectation of N C :

E[N C ] = p 0 r-C i=1 i (rρ) i+C (i + C)! + E[ν C ]| C=0 [1 + (r -C)(1 -ρ)] . (41)
Recall that the expression of E[ν C ] is given by the formula 14. Note also that this last expression 41 also gives the good result if we give the zero value to C:

E[N C ]| C=0 = rρ + E[ν C ]| C=0 . (42)
Similarly, we obtain a new expression for the probability that the response time is zero valid for the case C < (r -1):

P(S C = 0) = C-1 j=0 p j = p 0 C-1 i=0 (rρ) i i ! . ( 43
)
This last expression is not as nice as the previous formula 23 obtained in the case C ≥ (r -1). Note nevertheless that we can express the probability P(S C = 0) in a slightly different and potentially useful form if the number of servers r is important and if the parameter C is close to r:

P(S C = 0) = C-1 j=0 p j , = r-1 j=0 p j - r-1 j=C p j , RR n 9496 = P(S C = 0)| C=r -p 0 r-1 j=C (rρ) j j! , ( 44
) = 1 -ρ r r r r! p 0 (1 -ρ) -p 0 r-1 j=C (rρ) j j! . ( 45
)
It seems interesting to find the expression of E[N C ] (41) using an approach similar to that used to justify the formula 24 (page 17). For this we consider the bdp obtained from the one on Figure 10 by eliminating the states for which I(t) = 1. This gives us the bdp represented on Figure 12. Let us note q k the asymptotic probability of the state k of this new BDP, k = 0, 1, .... These probabilities are classically expressed in terms of the asymptotic probability q 0 :

0 1 r-C r-C+1 • • • • • • (C+1)µ λ λ (C + 2)µ
q k =            ρ k r k C ! (C + k) ! q 0 si k ≤ r -C , ρ k C ! r r-C r ! q 0 si k ≥ r -C .
Note that the transition rates retained in this new BDP have not been changed. Which allows us to state that

q i q j = p C+i p C+j , ∀i, j .
This means that the distribution q 0 , q 1 , q 2 , ... is a homothety of the partial distribution p C , p C+1 , p C+2 , .... Noting furthermore that (cf. 43):

∞ j=0 q j = 1 and ∞ j=0 p C+j = 1 - C-1 k=0 p k = P(S C > 0) ,
we deduce that

p C+j = q j (1 - C-1 k=0 p k ) = q j P(S C > 0) , ∀j . Inria Let us note E[ N C
] the expectation of the number of customers present in this pseudo queue. We have

E[ N C ] = ∞ j=1 jq j , = r-C j=1 jq j + ∞ j=r-C+1 jq j , = C ! q 0 r-C j=1 j (rρ) j (C + j) ! + C ! r r-C r ! q 0 ∞ j=r-C+1 jρ j , = C ! q 0 r-C j=1 j (rρ) j (C + j) ! + C ! r r-C r ! q 0 ∞ i=1 (i + r -C)ρ (i+r-C) , = C ! q 0 r-C j=1 j (rρ) j (C + j) ! + C ! r r-C r ! q 0 ρ (r-C+1) ∞ i=1 iρ i-1 + (r -C)ρ (r-C+1) ∞ i=1 ρ i-1 , = C ! q 0 r-C j=1 j (rρ) j (C + j) ! + C ! (rρ) r-C r ! q 0 ρ (1 -ρ) 2 + ρ(r -C) (1 -ρ) . (46)
Multiplying E[ N C ] by P(S C > 0), we obtain:

E[ N C ] × P(S C > 0) =   C ! q 0 r-C j=1 j (rρ) j (C + j) ! + C ! (rρ) r-C r ! q 0 ρ (1 -ρ) 2 + ρ(r -C) (1 -ρ)   (1 - C-1 k=0 p k ) . (47) Since q 0 (1 -C-1 k=0 p k ) is equal to p C and that p C = ρ C r C C ! p 0 , we get E[ N C ] × P(S C > 0) = p 0 r-C j=1 j (rρ) j+C (j + C) ! + (rρ) r r ! p 0 ρ (1 -ρ) 2 + ρ(r -C) (1 -ρ) , = p 0 r-C j=1 j (rρ) j+C (j + C) ! + p 0 r r r ! ρ r+1 (1 -ρ) 2 {1 + (r -C)(1 -ρ)} , = p 0 r-C j=1 j (rρ) j+C (j + C) ! + E[ν C ]| C=0 {1 + (r -C)(1 -ρ)} , = E[N C ]. ( 48
)
This confirms that E[N C ] is indeed equal to the product E[ N C ] × P(S C > 0)...

RR n 9496

We obtain the value of the expectation of the response time

E[R C ] by dividing E[N C ] by the arrival rate λ. E[R C ] = E[N C ] λ = 1 µ p 0 r-C i=1 i (rρ) i+C-1 (i + C)! + E[W C ]| C=0 [1 + (r -C)(1 -ρ)] . (49)
In figure 13, we have represented the evolution of the expectations E[W C ] and E[R C ] of a queue with 3 servers as a function of C, for C = 0. 1, 2, ..., 5, when ρ = 0.75 and µ = 0.5. As we could expect, these two expectations are decreasing according to the C level of anticipation. The expectation E[W C ] is given by the formula 16 whereas the expectation E[R C ] is given by the formula 49 if C < (r -1) and by the formula 18 if C ≥ (r -1). We observe that the difference (E[R C ] -E[W C ]), which is equal to the expectation of the service time E[S C ], decreases when C increases. For C = 0, the difference is equal at 1/µ = 2.

We can also verify, as the previous results have shown, that for a given value Concerning the probability P(W C = 0) that the waiting time is zero, the conditions of its determination are here the same as in the case C ≥ (r -1); we must therefore find the Inria same result as the one obtained page 19:

C ≥ r, E[R C ] = E[W C-r ].
P(W C = 0) = 1 -ρ (C+r) r r r! p 0 (1 -ρ) . ( 50
)
It is the same for the conditional expectation of waiting knowing that it exists:

E[W C |W C > 0] = 1 rµ(1 -ρ) . (51)
Again, this conditional expectation of waiting for the proactive queue is identical to that of the non-proactive queue.

On Figure 14, we have represented the probabilities that a service or waiting time is zero as a function of C, C = 0.1, 2, ..., 8. The number of servers is 3 while ρ = 0.75 and µ = 0.5. These two probabilities are increasing according to the C level of anticipation. The probabilities P(W C = 0) are given by the formula 50 whereas the probabilities P(S C = 0) are given respectively by the formula 43 if C < (r -1) and by the formula 23 if C ≥ (r -1). We can verify, as the previous results have shown, that P(S C = 0) = P(W C-r = 0) when C is greater than or equal to r (here equal to 3). Thanks to Little's formula, we determine the expectation of service time:

E[S C ] = E[N S C ] λ . ( 52
)
Remember that this unconditional expectation of service time takes into account customers whose service time is zero.

The determination of the conditional expectation of service time knowing that it is not zero is very easily done using the relation:

E[S C |S C > 0] = E[S C ] (1 -P(S C = 0)) , (53)
and values obtained from the expressions 43 and 52.

On reading the above, we find that the expressions of several performance indicators become more complex when we consider the case where the number of servers is greater than (C + 1). Unfortunately, this will be even more true when we consider the distribution of the sojourn time of any customer (that we will continue to call August). However, if we only consider the distribution of August's waiting time, the expression remains simple. Indeed, as long as at least one client is waiting, the overall service rate of the queue is equal to rµ, regardless of the value of C and therefore, knowing that August found j clients in front of him on arrival, j higher than (r -1), his waiting is distributed according to the Erlang distribution of order (j -r + 1) and of parameter rµ; of which the density of this Erlang distribution is written:

f X (t) = (rµ) (j-r+1) t j-r (j -r)! e -rµt , t ≥ 0 . ( 54
)
As in the case where C was greater than or equal to (r -1), the distribution of W C is a mixed distribution with a mass P(W C = 0) when t = 0. The computation of the distribution dF W C (t) is similar to that made page 21, and provides the result:

dF W C (t) = P(W C = 0)δ(t) + (1 -P(W C = 0))(rµ -λ)e -(rµ-λ)t . ( 55
)
where δ(t) denotes the Dirac's function. The only difference with the case where C was greater than (r -1) is that, at constant value of the r number of servers, the mass P(W C = 0) is here weaker, in the sense that if C 1 ≥ (r -1) and C 2 < (r -1), then P(W C 2 = 0) < P(W C 1 = 0) ( cf the equality 50). Now let's take a look at August's customer service distribution. If, on his arrival i ≥ 1 (and so j = 0), August's waiting is null and his service is instantaneous (the distribution of service is Dirac's distribution δ(0)). In a second step, let's examine the conditional service distribution knowing that upon arrival, August finds the system in the state i = 0 and j = 0, state in which C of the r servers are active. Since r is here strictly greater than (C + 1), the arrival of August causes the activation of an additional server, the number n a of active servers becomes equal to (C + 1) and the service rate becomes (C + 1)µ because August will benefit of the first service completed. The probability that August's service ends before a new client arrives is (C + 1)µ/((C + 1)µ + λ). If none of the (C + 1) active servers completes the execution of its request before the arrival of another client, a new server is activated (r is here greater than (C + 1) per hypothesis) and the overall service Inria rate that August will benefit becomes equal to (C + 2)µ. The probability that August's service ends before the arrival of a second new client is (C + 2)µ/((C + 2)µ + λ). There is no limit to the number of customers that August can see arriving before he leaves, but after (r -C -1) customer arrivals, August's service rate will equal rµ and will not change anymore. August's service distribution is here a particular phase-type distribution: it is the Coxian distribution of order (r -C). A schematic representation of this distribution is given on Figure 15. Recall that a Cox distribution can be seen as the distribution of a random sum of independent random variables distributed according to exponential distributions. The diagram of this figure 15 makes it easy to enter the semantics. Noting Λ i = λ + iµ, the service starts with a first phase whose duration is distributed according to an exponential distribution with rate Λ (C+1) . Then, eather the service is terminated (with probability (C + 1)µ/Λ (C+1) equal to (1 -λ/Λ (C+1) )), or it continues (with probability λ/Λ (C+1) ) by a second phase distributed according to the exponential distribution with rate Λ (C+2) ; and so more ... The sequence will have at most (r -C) phases; the last possible phase is of rate rµ. It is therefore a distribution of order (r -C) for the case C < (r -1). Noting successively:

a i △ = λ Λ (C+i)
, i = 1, 2, ..., r -C -1, and a r-C = 0,

A i △ = i-1 j=1 a j , j = 2, 3, ..., r -C, and A 1 = 1, B i △ = A i (1 -a i ), j = 1, 2, ..., r -C,
it is possible to show (cf. [START_REF] Marie | Introduction aux probabilités -Applications en informatique et télécommunications, fiabilité et gestion d'incertitudes[END_REF], appendix A.15) that the expectation of the rv S C following this Coxian probability distribution can be written:

E[S C ] = r-C i=1 A i Λ (C+i) (56)
whereas the Laplace transform of the distribution of this rv S C can be written:

RR n 9496 L S C (s) = r-C i=1 B i   i j=1 Λ (C+j) (Λ (C+j) + s)   .
It remains however to decompose these rational fractions into several simple elements to find the original, i.e., the density function of this Coxian distribution of order k. Now consider the case where, upon arrival, August finds himself in the p position, such that 1 < p ≤ r -C, he finds (C + p -1) active servers and his arrival causes the activation for an additional server, the number n a of active servers becomes equal to (C + p) and the service rate becomes equal to (C + p)µ. Before he can be served and thus leave the queue, August must wait until the (p -1) customers who precede him are themselves served. During all this time, the overall service rate of the queue may decrease following departures or increase to the maximum rate rµ due to arrivals. The August's service distribution is here a phase-type distribution of which we denote the states by the pairs (n, m) where n denotes August's position in the queue and where m denotes: eather the number of clients in the queue behind August when m is less than (r -C -1), or a number of clients in the queue behind August greater than or equal to (r -C -1) when m = (r -C -1).

Indeed, as soon as at least (r -C -1) customers are present in the queue behind August, the service rate of the queue will remain equal to rµ at least until the end of August's service. The initial state of this phase-type distribution is the state (p, 0). A schematic representation of this distribution is given in Figure 16.

The semi-formal expression of this phase-type distribution can be obtained by using the algorithm ace (cf. [START_REF] Marie | Transient analysis of acyclic Markov chains[END_REF]) 1 . It is also possible to obtain the formal expression of this phase law using a formal calculation software (for a reasonable number of states of the phase distribution). As an example, consider the special case where r = 6, C = 3, when August is in p = 2 position on arrival (cf. Figure 17). For didactic purposes, we will look for the formal expression of the phase-type distribution without using a formal calculation software.

After renumbering the states, we obtain the acyclic ctmc of the Figure 18. This enumeration of states has been chosen so that the infinitesimal generator A = (a ij ) of the Markov chain is upper triangular. This acyclic ctmc has 7 states of which the first six are transitory states and the seventh is the absorbing state.

Let us note P i (t) the probability that the ctmc is in state i at time t, i = 1, 2, ..., 7. The distribution function of the desired service distribution is equal to the probability P 7 (t). Let us recall that the diagonal elements a ii of A are defined by a ii =j =i a ij . Note a i = -a ii , i = 1, 2, ..., 7. a i is the total output rate of the state i. We use the integral form of the Kolmogorov equation which, knowing that the initial state is the state 1 with probability 1, is written:

1,0 1,1 1,2 1,r-C-2 1,r-C-1 2,0 2,1 2,2 2,r-C-2 2,r-C-1 p-1,0 p-1,1 p-1,2 p-1,r-C-2 p-1,r-C-1 p,0 p,1 p,2 p,r-C-2 p,r-C-1 • • • • • • • • • • • • • • • • • • • • • • • • • • • λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ (C + 1)µ (C + 2)µ (C + 3)µ (r -1)µ rµ 1 (C + 2)µ (C + 3)µ (C + 4)µ rµ rµ (C + 3)µ (C + 4)µ (C + 5)µ rµ rµ (C + p -1)µ (C + p)µ (C + p + 1)µ rµ rµ (C + p)µ (C + p + 1)µ (C + p + 2)µ rµ rµ
P i (t) = 1 {i=1} exp - t 0 a 1 du + t 0   k =i P k (x)a ki exp - t x a i du   dx . (57)
Using again the notations Λ i = λ + iµ, for i = 4, 5, 6, we get sucessively:

P 1 (t) = e -Λ 5 t . ( 58
) P 2 (t) = t 0 P 1 (x)λe -Λ 6 (t-x) dx , = λe -Λ 6 t t 0 e -Λ 5 x e +Λ 6 x dx = e -Λ 6 t λ (Λ 5 -Λ 6 ) -e -(Λ 5 -Λ 6 )x t 0 , = e -Λ 6 t λ (Λ 5 -Λ 6 ) 1 -e -(Λ 5 -Λ 6 )t = - λ µ e -Λ 6 t 1 -e -µt , = λ µ e -Λ 5 t -e -Λ 6 t .
The determination of the formal expressions of the probabilities P i (t), i = 3, 4, 5, 6, is presented in the appendix B (page 57). These formal expressions can be written: Figure 18: New denomination of the states of the studied acyclic ctmc.

P 3 (t) = λ (λ -µ) e -6µt - λ 2 µ(λ -µ) e -Λ 5 t + λ µ e -Λ 6 t . P 4 (t) = 5 e -Λ 4 t -e -Λ 5 t . P 5 (t) = 5 λ µ e -Λ 4 t -11 λ µ -λt e -Λ 5 t + 6 λ µ e -Λ 6 t . P 6 (t) = 6µλ (λ -µ) t - 5λ 2 (λ -µ) 2 + 6 + 5λ 2 µ(λ -2µ) - 11λ 2 µ(λ -µ) + 6λ µ e -6µt + - 5λ 2 µ(λ -2µ) e -Λ 4 t + 5λ 2 (λ -µ) 2 + 11λ 2 µ(λ -µ) - λ 2 t (λ -µ)
e -Λ 5 t -6 1 + λ µ e -Λ 6 t .

Inria

By not factoring more in the expression of the coefficients, it is quite easy to check that the state probability P 6 (t) is zero at the time origin.

The probability P 7 (t) which is none other than August's distribution of service time can be determined using the integral form of the Kolmogorov equation which, knowing that the state 7 is absorbing, is written:

P 7 (t) = + t 0 6 k=4 P k (x)a k7 dx , (59)
or, taking into account that at any time the ctmc is in a single state:

P 7 (t) = 1 - 6 k=1 P k (t) . ( 60
)
The beginning of August's service time distribution function is represented in Figure 19, when t ∈ [0, 1]. The different transient probabilities (P i (t), i=1, 2, .., 6) of the states of the phase-type distribution are also represented; these states are the transient states of the ctmc while state 7 is the absorbing state. The small time interval chosen for this figure (t ∈ [0, 1]) makes it possible to identify the various functions P i (t). The Figure 20 giving the answers when t ∈ [0, 3] shows that the distribution function of the service time tends to 1 while the probabilities of transient states tend to zero. Note that beyond an example of this size, it seems reasonable to use some software to calculate this distribution function of service time...

Let us now consider the case where, when he arrives, August is in the p position such as r -C < p ≤ r, he finds the r servers already active. He will take advantage of one of the services already in progress. The overall service rate remains, at least initially, equal to rµ. The distribution function of August's service time is still a phase-type distribution whose states are characterized by the pairs (n, m) as above. The initial state of this phasetype distribution is the state (p, 0), r -C < p ≤ r. From a state (n, m) such that n > 1 and m < (r -C -1), two transitions are possible: one transition to the state (n -1, m) with the rate min(r, (n + m + C))µ and another one to the state (n, m + 1) with the rate λ. When n = 1, the transition with the rate (1 + m + C)µ leads to the absorbing state signifying the end of the service and the departure of August. To know the distribution function of the service time, it is necessary to study this new acyclic ctmc to deduce the formal expression of this distribution function (as we did in the previous example) or the semi-formal expression using software. For example, taking the previous file where r = 6 and C = 3 and assuming this time that at his arrival, our August is in position p = 4, the diagram defining the phase service distribution is that in Figure 21. In this new situation, it is clear that the expectation of August s' service time is higher than the previous one. Let us now specify that, in this situation where p is greater than (r -C), an initial simplification of the phase distribution is possible. This simplification is presented hereafter with the study of the case where p is greater than r. Now, if upon arrival August finds himself in the p position such that p > r, he suffers a wait until the (p -r) customers at the front of the queue are served. Only then

P 7 (t)
Figure 19: Transient state probabilities of the states of the ctmc as a function of time t, t ∈ [0, 1], for p = 2, r = 6 and C = 3, with ρ = 0.8, λ = 3 and so µ = 0.625. will August's service begin. But while he was waiting, some customers were able to arrive and stand behind him in the queue. These arrivals will allow for a longer overall maximum service rate because August's service will start in a state (r, m) such that m > 0. To identify the distribution function of the service time, one must know what is the probability that August begins his service in a state (r, m), m = 0, 1, ..., r -C -1.

For this purpose, we consider the ctmc consisting of states of the type (n, m) such that {(n, m)|n = r + 1, r + 2, ..., p et m = 0, 1, ..., r -C -1.} and (r -C) absorbing states a m , m = 0, 1, ..., r-C -1, whose transition graph is represented in Figure 22. This ctmc is not irreducible, it contains many transient classes and (r -C) permanent classes corresponding to the absorbing states. The asymptotic state probability of an absorbing state a m is the probability that August's service begins in the state (r, m). In view of the number of transient classes, such a calculation would be very tedious; but we can circumvent this difficulty by considering an irreducible and ergodic cmtc which reproduces indefinitely the execution of the distribution function of August's waiting (cf. the graph of transition represented in Figure 23). Since the overall service rates are the same whatever the value of m when August leaves the (r + 1) position to end up in r position and start serving, we need only to know the relative values of the asymptotic state probabilities p r+1,m of the positive recurent states (r + 1, m), m = 0, 1, ..., r -C -1,. Indeed, if we know these asymptotic probabilities p r+1,m up to a multiplicative constant K, then the probability α m that August's service starts after a transition from state (r + 1, m) to state (r, m) will be equal to: Let us consider the partition of the state space of the cmtc represented in figure 23 such that:

α m = p r+1,m r-C-1 l=0 p r+1,l = Kp r+1,m r-C-1 l=0 Kp r+1,l . ( 61
) r+1,0 r+1,1 r+1,2 r+1,r-C-2 r+1,r-C-1 a0 a1 a2 ar-C-2 ar-C-1 r+2,0 r+2,1 r+2,2 r+2,r-C-2 r+2,r-C-1 p-1,0 p-1,1 p-1,2 p-1,r-C-2 p-1,r-C-1 p,0 p,1 p,2 p,r-C-2 p,r-C-1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ
E i = {(r + i, m)|m = 0, 1, ..., C -1} i = 1, 2, ..., (p -r) . (62)
Thus, we partition the states into (p -r) classes such that each class contains the states corresponding to a given position of August.

Let us introduce the (p -r) state cmtc Z(t) associated to the previous cmtc and such that Z(t) goes from state i to (i -1) when the initial cmtc transits from a state of class E i to a state of class E i-1 , if i = 2, ...., (p -r) and transitions from state 1 to state (p -r) when i = 1. Since all transition rates from states of class E i to states not belonging to class E i cause the cmtc to pass into a state of class E i-1 and are all equal to rµ, the cmtc Z(t) transits from state i to state (i -1) with the rate rµ, for i = 2, ...., (p -r). The same is true for the transition from state 1 to state (p -r). The fact that all the transition rates taking the initial chain from one class E i to E i-1 are all equal allows to obtain an exact aggregation. Examination of the transition graph of the Z(t) allows us to deduce that the stationary state probabilities π i of this irreducible chain with (p -r) states are all equal to (p -r) -1 . This leads, at the level of the cmtc of the Figure 23, that the sum of probabilities

Inria r+1,0 r+1,1 r+1,2 r+1,r-C-2 r+1,r-C-1 r+2,0 r+2,1 r+2,2 r+2,r-C-2 r+2,r-C-1 p-1,0 p-1,1 p-1,2 p-1,r-C-2 p-1,r-C-1 p,0 p,1 p,2 p,r-C-2 p,r-C-1 • • • • • • • • • • • • • • • • • • • • • • • • • • • λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ
r-C-1 l=0 p i,l is such that: r-C-1 l=0 p i,l = (p -r) -1 i = 1, 2, ..., (p -r) . ( 63
)
At the level of the cmtc Z(t), the duration of a cycle between two entries in state 1 (or between two entries in any given state) follows an Erlang distribution of order (p -r) and rate rµ. Note that we already know that the distribution of August's waiting knowing that he starts his waiting at position p, follows such a distribution (cf. the relation 54).

But the use of the present approach concerns the determination of the vector α in order to determine the initial probabilities of the service distribution of August, knowing that he starts his waiting at position p. Let's develop the approach in the particular case where, as before, r -C = 3 and where p = r + 7. Rename the states as noted on the representation of the graph of transition in Figure 24. Note p h the asymptotic probability of the state h, h, h = 1, 2, ..., 21. We consider the coefficients c h equal to a multiplicative constant close to asymptotic probabilities. The c h coefficients also verify the Chapmann-Kolmogorov equations. For example, we have:

Λ r c 2 = λc 1 , rµc 3 = λc 2 , Λ r c 4 = rµc 1 , Λ r c 5 = rµc 2 + λc 4 , rµc 6 = rµc 3 + λc 5 ,
Inria and so on. Note c j , j = 0, 1, ..., 6, the line vectors c j = (c 1+3j , c 2+3j , c 3+3j ) , j = 0, 1, ..., 6 .

(64)

Recall that ρ is equal to λ/rµ and note that

c 4 = 1 (1 + ρ) c 1 , c 5 = 1 (1 + ρ) c 2 + ρ (1 + ρ) c 4 = ρ (1 + ρ) 2 c 1 + 1 (1 + ρ) c 2 , c 6 = c 3 + ρc 5 = ρ 2 (1 + ρ) 2 c 1 + ρ (1 + ρ) c 2 + c 3 .
In matrix form, we get:

c 1 = c 0 A , (65)
where

A =       1 (1 + ρ) ρ (1 + ρ) 2 ρ 2 (1 + ρ) 2 0 1 (1 + ρ) ρ (1 + ρ) 0 0 1       (66)
Given the repetition of the transition rate pattern on the graph (all rates from a state i to a state i + 1 are equal to λ and all rates from a state i to a state i + 3 are equal to rµ), we obtain:

c j = c j-1 A , j = 1, ..., 6 . ( 67
)
or still c j = c 0 A j , j = 1, ..., 6 .

Let us introduce the unit column vector e, transpose of the row vector (1, 1, 1) and note that for any row vector (x 1 , x 2 , x 3 ), the matrix A has the property:

(x 1 , x 2 , x 3 )Ae = (x 1 , x 2 , x 3 )e = x 1 + x 2 + x 3 . ( 69
)
This gives rise to the property

c j e = c 0 e , j = 1, ..., 6 . ( 70
)
This result is not surprising. It follows from the identity expressed by the equation (63).

Moreover the vector c 0 can be written

c 0 = (c 1 , c 2 , c 3 ) = c 1 (1, ρ/(1 + ρ), ρ 2 /(1 + ρ)) , . (71)
This results in c 0 e = c 1 (1 + ρ) and more generally:

c j e = c 1 (1 + ρ) , j = 0, 1, ..., 6 . ( 72
)
By fixing arbitrarily the constant K equal to the inverse of the probability q 1 , the coefficient c 1 becomes equal to 1. Having computed the vector c 6 = c 0 A 6 , we obtain

α 0 = c 19 (c 19 + c 20 + c 21 ) , α 1 = c 20 (c 19 + c 20 + c 21 ) , α 2 = c 21 (c 19 + c 20 + c 21 ) . ( 73
)
where we already know that the denominator (c 19 + c 20 + c 21 ) is equal to (1 + ρ).

In this case where p = r + 7, the vector α is equal to: α = (0.0163, 0.0508, 0.9328) .

We notice that if August arrives in seventh position in the list of customers whose service has not yet started, he has a probability equal to 0.9328 that his service follows with certainty an Erlang distribution of order 6 and rate 6µ. Moreover, being seventh of the waiting customers when he arrives, the probability that August starts his service with at least one customer behind him is greater than 0.98.

More generally, assume that August starts his wait at position p, p > r. In this case, the vector α p = (α 0 , α 1 , α 2 ) of the initial probabilities of the service distribution is obtained using the vector c p-r-1 given by the relation

c p-r-1 = c 0 A p-r-1 , (75)
and from the expression of the vector sought as a function of the vector c p-r-1

α p = 1 (1 + ρ) c p-r-1 . (76)
Finally, for a given p value, the vector α p is calculated using the relation:

α p = 1 (1 + ρ) c 0 A p-r-1 . (77)
For example, in the case where p = r + 1 the vector α r+1 is collinear with the vector c 0 (because (p -r -1) = 0) and takes the value:

α r+1 = (0.5556, 0.2469, 0.1975) , (78)
which leads to a slightly higher expectation of service time than the Erlang distribution of the previous case (but in this case, the wait undergone by August will be a priori much shorter ...). For p = r + 2 and p = r + 3, the vectors α r+2 and α r+3 take the values:

α r+2 = (0.3086, 0.2743, 0.4170) , (79)
α r+3 = (0.1715, 0.2286, 0.5999) , (80)
values already very different from that of the vector α r+1 .

Inria

1,0 1, 1 1,2 2,0 2,1 2,2 3,0 3,1 3,2 5,0 5,1 5,2 6,0 6,1 6,2 a • •
• • • • • • • λ λ λ λ λ λ λ λ λ λ 5µ 6µ 6µ 6µ 6µ 6µ 6µ 6µ 6µ 6µ 6µ 6µ 6µ 6µ 6µ 4µ 5µ 6µ α0 α1 α2
Figure 25: Phase-type service time distribution when p > r, for a given p, in the particular case when r = 6 and C = 3. Before simplification.

At this point, the corresponding distribution is the phase type distribution with three potential initial states shown in Figure 25. The existence of several possible initial states distinguishes this distribution from the one already studied in the case where p was not greater than r and illustrated in Figure 18. But there is still a significant simplification to be made to simplify the calculation of the distribution function of this new distribution (simplification promised on page 33). To do this, let us take the Figure 25 and consider a first part of the states of the phase type distribution corresponding to the subset

{(n, m)|n > r -C -1 et m = 0, 1, ..., C -1} . (81)
Since the transition rates from any state (n, m) to a state (n-1, m) are all equal to rµ, this subset can be replaced by a linear chain comprising the (C + 1) states r, (r -1), ..., (r -C) provided that one can determine the transition probabilities from the state (r -C) to each of the states of the set

{(r -C -1, m)|m = 0, 1, ..., C -1} , (82)
i.e., to be able to determine the new vector α which we will note α + p . The transition rates from a new state j to a new state (j -1), j = r, (r -1), ..., (r -C + 1) are obviously equal to rµ. If August arrives at position p in the list of clients whose service has not yet begun, the vector α + p we are looking for is obtained by multiplying the previous vector (relative to the beginning of the service) by the matrix A C+1 . This gives the new vector α + p :

α + p = 1 (1 + ρ) c 0 A p-r+C . (83)
We thus obtain a phase type distribution with fewer states as shown in Figure 26 associated with the previously studied example. The number of states in the phase type distribution has thus been reduced from 18 to 10 states. Let us notice that in our example, the vectors α + r+1 , α + r+2 , α + r+3 , that take the values:

α + r+1 = (0.0529, 0.1176, 0.8295) , (84) α
+ r+2 = (0.0294, 0.0784, 0.8922) , (85)
α + r+3 = (0.0163, 0.0508, 0.9328) ,

values largely different from the vectors α r+1 , α r+2 and α r+3 calculated previously. This is not surprising if we consider the simple relationship between α + p and α p :

α + p = α p+C+1 . (87)
As we announced on page 33, this simplification is applicable in the situation where our client August is in position p, such that r -C < p ≤ r. In this case, the Erlang distribution is of order (p -r + C + 1) and the vector α + p provides the transition probabilities from the state (r -C) to the states ((r -C -1), l), l = 0, 1, 2. The relation 87 still allows to determine this vector α + p . Let us add that if p = (r -C), the arrival of August will cause the activation of the last inactive server and will increase the global service rate to rµ, which allows a small simplification: the Erlang distribution is of order 1 (it is the exponential distribution of rate rµ) and the α + p vector is equal to c 0 /(1 + ρ). Note that if α p does not depend on C (cf. 77), α + p does depend on C. Depending on the number of states, the distribution function of this new service distribution can be determined according to one of the methods already mentioned.

When p tends to infinity, August's service distribution tends to the limiting distribution corresponding to the situation where at least (r -C) customers are present behind him at the beginning of his service. In this limiting case, the service distribution is the Erlang distribution of order r and rate rµ. If we now look for the distribution of August's service time knowing only that he waited (i.e., p > r), an event of probability P(W > 0), we need to determine the vector α conditional on this event. Knowing that the unconditional probability that August is in the queue in p-th position is equal to p (C+p-1) , we can give an expression for the vector α conditional on the event [p > r], a vector denoted α att :

α att = 1 P(W > 0) 1 (1 + ρ) ∞ p=r+1 p (C+p-1) c 0 A p-r-1 . (88)
Or, using the expression of P(W > 0) given in 30:

α att = 1 P(W > 0) 1 (1 + ρ) ∞ p=r+1 ρ (C+p-1) r r r ! p 0 c 0 A p-r-1 , = 1 P(W > 0) r r p 0 r !(1 + ρ) ∞ p=r+1 ρ (C+p-1) c 0 A p-r-1 , = 1 P(W > 0) r r p 0 ρ C+r r !(1 + ρ) c 0 ∞ p=r+1 ρ p-r-1 A p-r-1 , = (1 -ρ) (1 + ρ) c 0 ∞ l=0 (ρA) l , = (1 -ρ) (1 + ρ) c 0 ∞ l=0 B l ,
where B = If the spectral radius of the matrix B is less than 1 (i.e. if the eigenvalues of the matrix B are less than unity in modulus), we can still write

α att = (1 -ρ) (1 + ρ) c 0 (I -B) -1 . (89)
Let us write matrix B using the version of A given by 66:

B =       ρ (1 + ρ) ρ 2 (1 + ρ) 2 ρ 3 (1 + ρ) 2 0 ρ (1 + ρ) ρ 2 (1 + ρ) 0 0 ρ      
The determinant of the matrix (B -λI

) is equal to (ρ -λ) ρ (1 + ρ) -λ 2 .
The characteristic equation of the matrix B thus admits the simple root λ 1 = ρ and the double root

λ 2 = ρ (1 + ρ)
. These roots are by definition positive and less than 1 (the spectral radius of the matrix B is equal to ρ). The matrix (I -B) -1 exists. We can therefore compute it by starting from:

(I -B) =       1 (1 + ρ) - ρ 2 (1 + ρ) 2 - ρ 3 (1 + ρ) 2 0 1 (1 + ρ) - ρ 2 (1 + ρ) 0 0 1 -ρ      
The determinant of the matrix (I -B) is equal to (1 -ρ)

1 (1 + ρ) 2 .
We finally obtain:

Inria (I -B) -1 =         (1 + ρ) ρ 2 ρ 3 (1 -ρ) 0 (1 + ρ) ρ 2 (1 -ρ) 0 0 1 (1 -ρ)        
We can therefore use the relation 89 to determine the vector α att , knowing just that August waited. We obtain:

α att = ((1 -ρ), ρ(1 -ρ), ρ 2 ) , (90)
If we want to look for the service time distribution of August knowing that he waited, it is better to use the simplified representation and the vector α + att computable thanks to the relation:

α + att = α att A C+1 , (91)
Recall that if α att does not depend on C (cf. 89), of course, α + att depends on C. In the numerical example considered above, ρ is equal to 0.8. This gives the vectors: α att = (0.2000, 0.1600, 0.6400) α + att = (0.0191, 0.0491, 0.9318) .

In the Figures 27 and28, we represent respectively the beginning and the end side of the distribution function of the service time random variable in the following situations: p = r + 1, W > 0 (we just know that p is greater than r), as well as the limit when p tends to infinity. We can see that the three corresponding functions are close (hence the use of two distinct figures). The values of the parameters of the queue are those of the previous example r = 6 and C = 3. These functions were determined using the semi-formal software ace. It is also possible to evoke the response time distribution but, in order to restrict (and simplify) the discourse, we will present the approach by considering only the example already considered, the one where r = 6 and C = 3.

Let's start by looking for the distribution of the response time knowing that August starts his stay by waiting (p > r); we already know (cf. relation 55) that the conditional Inria distribution of the waiting is an exponential distribution of rate (rµ -λ). The conditional distribution of the response time knowing that August starts his stay by waiting can thus be calculated as a phase distribution starting with an initial phase of rate (rµ -λ) and continuing with the first phase of the phase distribution given in Figure 26. The diagram of the phase distribution corresponding to the response time knowing that August begins his stay with a wait (knowing W > 0 and thus p > r) is represented in Figure 29. The coefficients α + att,i , i=1,2,3, are the components of the vector α + att . In the Figure 30, we represent the response time distribution function knowing W > 0, with r = 6 and C = 3. In this figure, we have added (in dotted line) the distribution function of the Erlang limit service distribution of order r and rate rµ in order to appre-ciate the experienced wait. Again, these functions were determined using the semi-formal software ace.

Let us specify that the graphical representation of this distribution function of the residency distribution knowing W > 0 in the case C = 4 (and all data being equal) would be almost undistinguishable from that where C = 3. Indeed, for this conditional distribution, the parameter C only intervenes on the service distribution which, thanks to the value of the vector α + att , is already very close to the Erlang-r distribution. In the present case where C > (r -1), the representation of the (unconditional) residence distribution F R C (t) would be tedious but not infeasible. We already know that the distribution of R C admits a mass at t = 0 equal to P(S C = 0). When the service is no longer immediate and August finds himself in position p, p ≤ r, with probability p C-1+p , his stay time is equal to his service time, whose distribution function for any value of p, p = 1, 2, ..., r, we have seen previously. Finally, it remains to add to this sum the conditional distribution function F C,W >0 (t) that we have just determined (see Fig. 30), weighted by the probability P(W > 0). This can be summarized as:

F R C (t) = P(S C = 0)u(t) + r q=1 p C-1+q F R C,p=q (t) + P(W > 0)F R C,W >0 (t) . ( 92
)
Inria 4 Multi-server Queue with state-dependent rates

In this section, we study proactive queues where the arrival rate λ(i, j) and the service rate µ(i, j) can depend on the variables i and/or j. For these queues, we will use the notation λ(i, j)/µ(i, j)/ < C >. The study of their performances requires the writing of algorithms but before presenting them, let us give some concrete examples of this family of models. The first example concerns the study of a finite population of N individuals where the arrival rate λ(i, j) depends only on the number j of customers present in the queue:

λ(j) = (N -j)λ e j ≤ N . (93)
The second example corresponds to the case where the service is not carried out on the place of contact with the customers but is carried out by a distant (and possibly complex) process on request of the site receiving the customers. Without addressing the technique of obtaining it for the moment, we can assume that we are able to determine a sequence of service rates {µ(i, j)} depending only on the number (C + j -i) of requests issued and not satisfied. The variable k defined by the relation 1 of the page 4 can still be used to characterize the state of the queue since the knowledge of k makes it possible to find the couple (i, j) thanks to the relations 2 and 3 (page 4). In order not to make the expressions too long, we will respectively note λ(k) and µ(k) the rates λ(i(k), j(k)) and µ(i(k), j(k)) associated to the value of k.

The study of the performance of these queues starts with the determination of the steady state probabilities. These probabilities satisfy the relations:

p k = λ(k -1) µ(k) p k-1 k > 0. ( 94
)
Let us set

f (k) △ = λ(k -1) µ(k) and H(k) △ = k n=1 λ(n -1) µ(n) .
This allows us to express all probabilities in terms of the probability p 0 .:

p k = H(k)p 0 k > 0. ( 95
)
The constraint k > 0 can be ignored provided that H(0) △ = 1. If the value of the arrival rate λ(k) becomes zero for k ≥ K lim , the number of states in the queue is equal to (K lim + 1) and the probability p 0 is obtained thanks to the normalization equation:

p 0 = K lim k=0 H(k) -1 .
When j is strictly positive, the value of i is always zero and the arrival rate λ(k) can only depend on the number j of customers present in the queue. A priori, the subsequence {λ(k)} is then decreasing, eventually with a null limit value for j ≥ J lim . On the other hand, the sub-sequence {µ(j)}, j > 0, is assumed here to be non-decreasing. This assumption is (naturally) satisfied in all real systems modeled by queues of infinite capacity and admitting an asymptotic state probability distribution (in theory, one could also accept a sub-sequence decreasing to a sufficiently high lower bound value to satisfy the ergodicity conditions).

If the limit value of λ(k) is not zero and does not have a simple formal expression, things get complicated because it is out of the question to calculate the function H(k) until infinity. We must then accept a slight approximation. We suppose that there exists a value α, α < 1, and a value K α such that

λ(k -1) µ(k) ≤ α k ≥ K α . (96)
Then, we consider a first queue of limited capacity such that k = 0, 1, ..., K α and a second queue of infinite capacity whose function similar to the function f (k), denoted f (k), is defined as a function of f (k) as follows

f (k) △ = f (k) if k ≤ K α , f (K α ) if k ≥ K α .
The previous definition allows us to obtain the expression of the function similar to the function H(k), noted H(k):

H(k) △ = H(k) if k ≤ K α , H(K α )f (K α ) k-Kα if k ≥ K α .
The present approximation relies on the fact that, if the boundary value of λ(k) is not zero, the value α is chosen to be sufficiently small so that K α is sufficiently large and the majoration of f (k) by f (K α ), k ≥ K α , does not lead to a significant error on the computed numerical results.

For the queue of limited capacity, the probability, noted p 0 , a upper bound of p 0 is obtained by the normalization equation:

p 0 = Kα k=0 H(k) -1 .
It is easy to see that p 0 is strictly greater than p 0 . Note that to avoid an overflow during its computation, this approximate probability p 0 can be obtained by the recursion

p 0 (n) = p 0 (n -1) 1 + H(n)p 0 (n -1) n = 1, 2, ...K α avec p 0 (0) = 1. (97)
Note that the sequence p 0 (n) is strictly decreasing since (1 + H(n)p 0 (n -1)) is strictly greater than 1.

For the approximate queue of infinite capacity, the probability p 0 , a lower bound of p 0 , is written:

p 0 =   Kα k=0 H(k) + ∞ k=Kα+1 H(k)   -1 , =   Kα k=0 H(k) + ∞ k=Kα+1 H(K α ) (f (K α )) k-Kα   -1 , =   Kα k=0 H(k) + H(K α )f (K α ) ∞ k=Kα+1 (f (K α )) k-(Kα+1)   -1 , = Kα k=0 H(k) + H(K α )f (K α ) (1 -f (K α )) -1 because f (K α ) < 1. Note that if f (k) = f (K α ) for k ≥ K α , then p 0 = p 0 .
In the opposite case ( H(k) > H(k)), the probability p 0 is strictly greater than p 0 . Let us now examine the ratio of p 0 to p 0 :

p 0 p 0 = Kα k=0 H(k) + H(K α )f (K α ) (1 -f (K α )) -1 × Kα k=0 H(k) , =   1 + H(K α )f (K α ) Kα k=0 H(k) (1 -f (K α ))   -1 = 1 + p 0 H(K α )f (K α ) (1 -f (K α )) -1 , = (1 -f (K α )) (1 -f (K α )) + p 0 H(K α )f (K α ) , (98)
= 1 - p 0 H(K α )f (K α ) (1 -f (K α )) + p 0 H(K α )f (K α ) .
Let us put:

ε(K α ) △ = p 0 H(K α )f (K α ) (1 -f (K α )) + p 0 H(K α )f (K α ) (99)
and let us illustrate these approximations on a proactive M/M/r/ < C > queue subjected to a load factor of 0.7 and having 5 servers. For this queue, the function f (k) becomes constant very quickly (for k ≥ r) and as soon as K α is chosen to be equal to or greater than r, the minor value is equal to the true value p 0 . When K α = 40, we obtain p 0 = 0.0259 and ε(40) = 1.0019 × 10 -6 . When K α = 50, we obtain p 0 = 0.0259 and ε(50) = 2.8302 × 10 -8 ; while H(40) = 1.6580 × 10 -6 and H(50) = 4.6835 × 10 -7 . We notice that the two values of K α correspond to two very close p 0 values. The observation of the formula 97 and the very low values of H(40) and H(50) allow to understand the reason of this very strong proximity of the values. We can say that these approximations provide near-accurate results with low computational cost.

The method of calculation of the upper bound p 0 consists in calculating simultaneously p 0 (n), via the recurrence formula 97 and ε(n) in terms of p 0 (n) according to the formula

ε(n) △ = p 0 (n)H(n)f (n) (1 -f (n)) + p 0 (n)H(n)f (n) ,
until ε(n) becomes lower than a chosen bound (eg, 10 -5 ).

The relation 98 allows then to obtain the value of the lower bound p 0 according to p 0 . Once the probabilities p 0 and p 0 of the approximate distributions have been calculated, we can calculate the other probabilities p k and p k

p k = H(k)p 0 , k = 1, 2, ..., K α et p k = H(k)p 0 , k = 1, 2, 3, ... .
By construction of the probabilities p 0 and p 0 , we have indeed:

Kα k=0 p k = 1 , ∞ k=0 p k = 1 . (100)
Let us also note that

p k p k = p 0 p 0 , k = 1, 2, ..., K α and p k p k = p 0 p 0 , k = 1, 2, 3, ... . (101)
Note that in the very special case where H(k) = H(k), k ≥ K alpha , the exact probability distribution {p k }is easy to calculate since p k = p k , k ≥ 0. In the general case where H(k) > H(k), k > K α , note that there exists a value K cross where the probability p k becomes greater than the probability p k :

K cross = min k {p k > p k } .
Before turning to performance indicators, let's recall the notion of stochastic dominance.

4.1.

Let X and Y be two real rv admitting F X (.) and F Y (.) as distribution function. The rv X is stochastically dominated by the rv Y if

F X (u) ≥ F Y (u) , ∀u ∈ IR . (102)
This stochastic dominance (of the first order) is generally noted X ≺ st Y . This dominance means that small values are more likely for rv X than for rv Y . If the expectations of X and Y exist, it is easy to verify that the first order stochastic dominance leads to:

Inria E[X] ≤ E[Y ] ; (103)
It is sufficient to use the fact that here the expectation E[X], resp. E[Y ], is also equal to the integral of the complementary function of its distribution function (i.e., (1 -

F X (u))), resp. (1 -F Y (u)))
It is also possible to show that if an X is stochastically dominated to the first order by a rv Y and if g is a increasing function, then:

E[g(X)] ≤ E[g(Y )] . (104)
Let us specify that if the inequality is strict at the level of distribution functions (relation 102), it also becomes strict at the level of expectations (relations 103 and 104).

We will also need to compare expectations E[g(.)] when g is a decreasing function. In such a case, letting us put s △ = -g, s is an increasing function and we obtain, if the rv X is stochastically dominated at first order by the rv Y

, E[s(X)] ≤ E[s(Y )] and thus E[-g(X)] ≤ E[-g(Y )] or -E[g(X)] ≤ -E[g(Y )], that is: E[g(X)] ≥ E[g(Y )]
(case when g is decreasing) .

Let us return to the relationship 101. In the general case where H(k) > H(k), k > K α , this one implies that p k < p k , k = 0, 1, 2, ..., K α and that p k < p k , k = 0, 1, 2, ... this implies that:

k j=0 p j > k j=0 p j , ∀k ∈ IN et k j=0 p j > k j=0 p j , ∀k ∈ IN . (106)
In other words, the distribution {p k } k≥0 stochastically dominates (to first order) the distribution {p k } k≥0 as well as the exact (but not calculated) distribution {p k } k≥0 . As it is also easy to verify that the exact distribution stochastically dominates the {p k } k≥0 , we finally obtain:

{p k } k≥0 ≺ st {p k } k≥0 ≺ st {p k } k≥0 . (107)
Let us insist on the fact that the stochastic dominance of {p k } k≥0 on {p k } k≥0 exists only in the general case where H(k) > H(k), k > K α since in the special case, the distribution {p k } k≥0 corresponds to the exact distribution. Let's now look at the performance indicators, and among these, first of all those that take the form of:

IP i = ∞ k=k inf Y i (k)p k , (108)
where Y i (k) is an increasing function of k. This is the case for several indicators such as

E[ν C ], E[N S C ], E[N C ] ...
If IP i, (p 0 ) , resp. IP i, (p 0 ) , denotes the performance indicator IP i calculated with the distribution {p k } k≥0 , resp. {p k } k≥0 , the relations 104 and 107 allow us to write IP i, (p 0 ) ≤ IP i ≤ IP i, (p 0 ) . But this search for bounds is only useful for the general case where H(k) > H(k), k > K α . And in this general case, the strict inequalities obtained by the relation 106 allow to write here:

IP i, (p 0 ) < IP i < IP i, (p 0 ) (general case) . (109)
Let us apply the method by considering the evaluation of the expectation of the number of pending requests in such a proactive queue.

Let us denote

E[ν C ] (p 0 ) (resp. E[ν C ] (p 0 )
) the expectation of the number of pending requests in the approximate queue with limited (resp. infinite) capacity. The lower bound of E[ν C ] is obtained from the relation 14 of the page 15:

E[ν C ] (p 0 ) = Kα k=C+r+1 (k -C -r)p k = p 0 Kα k=C+r+1 (k -C -r)H(k) , (110)
The upper bound of E[ν C ] is obtained from the same relation 14:

E[ν C ] (p 0 ) = ∞ k=C+r+1 (k -C -r) p k , = Kα k=C+r+1 (k -C -r) H(k) p 0 + ∞ k=Kα+1 (k -C -r) H(k) p 0 , = p 0   Kα k=C+r+1 (k -C -r)H(k) + ∞ k=Kα+1 (k -C -r)H(K α ) (f (K α )) k-Kα   , = p 0   1 p 0 E[ν C ] (p 0 ) + H(K α ) ∞ k=Kα+1 (k -K α ) (f (K α )) k-Kα +(K α -C -r)H(K α )f (K α ) ∞ k=Kα+1 (f (K α )) k-Kα-1   , = p 0 1 p 0 E[ν C ] (p 0 ) + H(K α )f (K α ) (1 -f (K α )) 2 + (K α -C -r) H(K α )f (K α ) (1 -f (K α )) , = p 0 p 0 E[ν C ] (p 0 ) + p 0 H(K α )f (K α ) (1 -f (K α )) 1 (1 -f (K α )) + (K α -C -r) . (111)
Let us illustrate the method by considering the evaluation of the expectation of the number of pending requests in a proactive queue M/M/10/ < 5 > with an arrival rate subject to a discouraging distribution. More precisely, the clients that show up at the entrance of the queue are accepted with a probability depending on the number of clients already present. If j denotes this number of customers, the acceptance probability is given Inria here by 1/(1+β(j -r) + ) where β is chosen less than 1. The initial arrival rate is poissonian (noted λ b ) and can be greater than 1/rµ. The sequence of arrival rates is thus written: This simple example shows that this approach allows to estimate very precisely the different performance indicators with a low calculation cost.

λ(k) = λ b 1 + β(k -C -r) + k ≥ 0. ( 112

Conclusion

In this paper, just as we know the formulas for determining the performance of non proactive queues, we sought to determine the set of results that would allow us to evaluate the performance of proactive queues. We have seen that in the end, some results are simply obtained while others are more difficult. The fundamental difference between the proactive and non proactive models is that the services of the proactive queues are potentially usable by any client. This window of anticipation of C services allows lucky customers to be served instantly and allows all customers to benefit from shorter response times. Nevertheless, we determined the performance in the more complex case where the number of servers is greater than (C + 1). We also proposed an efficient approach to evaluate a queue where the arrival and service rates depend on the number of clients in the queue. We conclude by reminding that proactive models can be used in different industries such as telecommunications, industrial production and integrated logistical support. (1 -ρ) ,

= 1 {r>C+1} p 0 r-C-1 i=1 i (rρ) i+C (i + C)! + p 0 r r r ! ρ r ∞ i=0 (i + r -C)ρ i -ρ C ∞ i=1 (i + r)ρ i + +rp 0 r r r ! ρ C+r+1 (1 -ρ) , = 1 {r>C+1} p 0 r-C-1 i=1 i (rρ) i+C (i + C)! + p 0 r r r ! ρ r (r -C) + (1 -ρ C ) ∞ i=1 iρ i + (r -C -rρ C ) ∞ i=1 ρ i + +rp 0 r r r ! ρ C+r+1 (1 -ρ) , = 1 {r>C} p 0 r-C i=1 i (rρ) i+C (i + C)! + p 0 r r r ! ρ r+1 (1 -ρ C ) (1 -ρ) 2 + (r -C -rρ C ) (1 -ρ) + rp 0 r r r ! ρ C+r+1 (1 -ρ) , = p 0 r-C i=1 i (rρ) i+C (i + C)! + p 0 r r r ! ρ r+1 (1 -ρ C ) (1 -ρ) 2 + (r -C) (1 -ρ) . ( 114
)
Inria B Appendix : Determination of the formal expressions of the probabilities P i (t), i = 3, 4, 5, 6, related to the phase type distribution in the particular case where p = 2, r = 6 and C = 3 (situation where 1 < p ≤ r -C)

Recall the expression of the already calculated probabilities P 1 (t) and P 2 (t):

P 1 (t) = e -Λ5t , (115)
P 2 (t) = λ µ e -Λ5t -e -Λ6t .
By continuing to use the Kolmogorov equation (cf. 57) and the notation Λ i = λ+iµ, we obtain successively: t -6λ 2 (λ -µ) 2 1 -e -(λ-µ)t + 6 1 -e -λt + + 5λ 2 µ(λ -2µ)

P 3 (t) =
1 -e -(λ-2µ)t -11λ 2 µ(λ -µ)

1 -e -(λ-µ)x + + 6λ µ 1 -e -λtλ 2 (λ -µ) 2 (1 + (λ -µ)t)e -(λ-µ)t +

λ 2 (λ -µ) 2 , = 6µλ (λ -µ) t - 5λ 2 (λ -µ) 2 + 6 + 5λ 2 µ(λ -2µ) - 11λ 2 µ(λ -µ) + 6λ µ e -6µt + - 5λ 2 µ(λ -2µ) e -Λ4t + 5λ 2 (λ -µ) 2 + 11λ 2 µ(λ -µ) - λ 2 t (λ -µ) e -Λ5t -6 1 + λ µ e -Λ6t .
By not further factoring the expression of the coefficients, it is quite easy to verify that the state probability P 6 (t) is zero at the origin.

Inria

C Appendix : A reminder of the ace method Let us denote Q = {q ij } the infinitesimal generator of the acyclic cmtc and q i the global departure rate of a state i. Let us call these global departure rates poles because they correspond to the poles of the Laplace transforms of the expressions of the transient state probabilities of the acyclic cmtc states. If there are N different poles in this chain, the list (γ 1 , γ 2 , ..., γ n , ..., γ N ) designate these different poles and γ(n) designate the value of the γ n pole. The goal of the presented method is to determine the transient state probability of a state j in the following semi-formal form:

P j (t) = n e γ(n)t k a(j, n, k)t k , (116)
where the coefficients a(j, n, k) are expressed in numerical form. Let us note Γ(j) the set of poles present in the expression of the probability P j (t). We can write more precisely this probability:

P j (t) = N n=1 1 {γn∈Γ(j)} e γ(n)t Kj (n) k=0 a(j, n, k)t k , (117)
where the integers K j (n) are assumed to be known. We assume that the states are ordered so that the matrix Q is upper triangular. Let J(i) be the set of "parent" states of state i. A state j is parent of a state i if q ji > 0. Let us note Γ(J(i)) the set of poles that appear in at least one transient state probability of the "parent" states of state i.

Using the integral form of the Kolmogorov equation, the expression for the probability of state i is written:

where γ * = q i . As the state probabilities will a priori share more and more common poles, we reverse the double sum:

1 {γn∈Γ(j)} 1 {k≤Kj (n)} a(j, n, k)q ji , (120)
for k = 0, 1, ..., K i (n), where

K i (n) = max{K j (n)|j ∈ J(i)} . (121)
This allows us to rewrite the expression of P i (t) as P i (t) = P i (0)e γ * t + x k e (γ(n)-γ * )x dx , (122)

Figure 1 :

 1 Figure 1: Schematization of the proactive single-server queue.

Figure 2 :

 2 Figure 2: cmtc associated to the proactive single-server queue.

Figure 3 :

 3 Figure 3: Probabilities of null service times (x) or of null waiting times (+) as a function of C, for C=0, 1, 2, ..., 5; with ρ = 0.75 and µ = 0.5.

Figure 4 :

 4 Figure 4: Expectations of waiting times (x) and of response times (+) as a function of C, for C=0, 1, 2, ..., 5; with ρ = 0.75 and µ = 0.5.

Figure 5 :

 5 Figure 5: Cumulative distribution function of R C , F R C (t),as a function of C, for C=0. 1, 2, 3, and 4; with ρ = 0.75 and µ = 0.5.

Figure 6 :

 6 Figure 6: Schematization of the proactive multi-server queue.

Figure 7 :

 7 Figure 7: Case C ≥ (r -1). Illustration of relations between i, j, n a and n pa as a function of k.

Inria

Figure 8 :

 8 Figure 8: Case C < (r -1). Illustration of relations between i, j, n a and, n pa as a function of k.

Figure 9 :

 9 Figure 9: ctmc associated to the proactive multi-serveur queue (case C ≥ (r -1)).

Figure 10 :

 10 Figure 10: ctmc associated to the proactive multi-serveur queue (case C < (r -1)).

Figure 12 :

 12 Figure 12: BDP associated with the proactive multi-server pseudo-file (case C < (r -1)).

Figure 13 :

 13 Figure 13: Expectations of waiting times (+) and of response times (x) as a function of C, for C = 0, 1, 2, ..., 5. Case with 3 servers, when ρ = 0.75 and µ = 0.5.

Figure 14 :

 14 Figure 14: Probabilities of null service times (x) and of waiting times (+) as a function of C, for C = 0, 1, 2, ..., 8. Case with 3 servers, when ρ = 0.75 and µ = 0.5.

rµ 1 Figure 15 :

 115 Figure 15: Coxian distribution of order (r -C).

Figure 16 :

 16 Figure 16: Phase-type service time distribution, case 1 < p ≤ r -C.

Figure 17 :

 17 Figure 17: Phase-type service time distribution in the particular case where p = 2, r = 6 and C = 3.

Figure 20 :Figure 21 :

 2021 Figure 20: Transient state probabilities of the states of the ctmc as a function of time t, t ∈ [0, 3], for p = 2, r = 6 and C = 3, with ρ = 0.8, λ = 3 and so µ = 0.625.

Figure 22 :

 22 Figure 22: Non irreducible ctmc modeling August's waiting time distribution function.

Figure 23 :

 23 Figure 23: Irreducible ctmc modeling the execution of an infinity of independent realizations of the August's waiting time distribution function.

Figure 24 :

 24 Figure 24: Irreducible cmtc; case p = r + 7.

Inria

Figure 26 :

 26 Figure 26: Phase-type service time distribution when p > r, for a given p, in the particular case when r = 6 and C = 3. After simplification.

Figure 27 :

 27 Figure 27: Start of the service time distribution function, r = 6, C = 3; with ρ = 0.8, λ = 3 and so µ = 0.625.

Figure 28 :

 28 Figure 28: End side of the service time distribution function, r = 6, C = 3; with ρ = 0.8, λ = 3 and so µ = 0.625.

1 .

 1 6µ α + att,2 .6µ α + att,3 .6µ

Figure 29 :

 29 Figure 29: Phase-type response time distribution knowing W > 0, in the particular case where r = 6 and C = 3.

Figure 30 :

 30 Figure 30: Distribution function F R C,W >0 (t) of the response time random variable knowing W > 0, when r = 6, C = 3; with ρ = 0.8, λ = 3 and thus µ = 0.625. The dotted curve represents the distribution function of the Erlang limit service ctmc of order r and rate rµ.

)

  With the data r = 10, C = 5, µ = 0.2, λ b = 3 and β = 0.2, we obtain:4, 7686 < E[ν C ] < 4, 7703 et K α = 36 avec la consigne ε ≤ 10 -4 , 4, 7701 < E[ν C ] < 4, 7702 et K α = 39avec la consigne ε ≤ 10 -5 . (113)

A

  Appendix : Computation of the expectation E[N S C ] (situation where 1 < p ≤ r -C) E[N S C ] = E[N C ] -E[ν C ] , = ∞ j=C+1 (j -C)p j -∞ j=C+r+1 (j -C -r)p j , = C+r j=C+1 (j -C)p j + ∞ j=C+r+1 (j -C -(j -C -r))p j , r -C)ρ i+r -

t 0 P 2 0 P 1 = 5µe -Λ4t t 0 e 0 , 0 P 2 0 P 4 0 ++5λte -Λ5t e µx µ t 0 - 0 P 3 0 P 5

 0201000204000305 (x)λe -6µ(t-x) dx , -e -Λ6t e +6µx dx , (x)5µe-Λ4(t-x) dx , -Λ5x e +Λ4x dx = e -Λ4t 5µ (Λ 5 -Λ 4 ) -e -(Λ5-Λ4)x t = e -Λ4t 5µ (Λ 5 -Λ 4 ) 1 -e -(Λ5-Λ4)t = 5e -Λ4t 1 -e -µt ,= 5 e -Λ4t -e -Λ5t .P 5 (t) = t (x)6µe -Λ5(t-x) dx + t (x)λe -Λ5(t-x) dx , = 6λe -Λ5t t 0 e -Λ5x -e -Λ6x e Λ5x dx + +5λe -Λ5t t 0 e -Λ4x -e -Λ5x e Λ5x dx , Λ5t -6λe -Λ5t -e -µx µ t 5λte -Λ5t , = 6λte -Λ5t -6 λ µ e -Λ5t1 -e -µt + +5 λ µ e -Λ5t (e µx -1) -5λte -Λ5t , (x)6µe -6µ(t-x) dx + t (x)λe -6µ(t-x) dx ,

P

  i (t) = P i (0)e γ * t + , n, k)x k   q ji e γ * (t-x) dx ,

P

  i (t) = P i (0)e γ * t + , n, k)x k   q ji e γ * (t-x) dx , (119) Let's define ã(i, n, k) as the sum ã(i, n, k) = j∈J(i)

N n=1 1

 1 {γn∈Γ(J(i))} Ki(n) k=0 e γ * t ã(i, n, k) t 0

RR n 9496

A reminder of this method is presented in appendix C Inria

Inria

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

For the calculation of the integral, we must distinguish two possible cases. If γ(n) = γ * , we immediately obtain:

whereas if γ(n) = γ * , we obtain after having put α = (γ(n) -γ * ):

x k e αx dx = e αt k m=0

We obtain the new expression for P i (t):

.

To simplify the calculation algorithm, it is useful to note that

This allows us to calculate the new coefficients relative to a pole γ(n) = γ * in the following way:

In order to determine the coefficient a(i, n * , 0) relative to the γ * pole, note first that, according to the relation 126,

and so, that according to the relation 125,

It immediately follows that:

Inria

If γ * is already in the set of poles of the parents of the state i, the order of the polynomial expression associated to the pole increases by one unit (according to the relation 125) and, according to the relation 123, its new coefficients can be obtained by the recurrence:

Note that for poles γ(n

A schema of the algorithm is proposed table 1 and table 2. The presented algorithm uses many input data. First of all, the infinitesimal generator of the acyclic cmtc and the vector of initial state probabilities. The data N S and N P give the number of states and the number of poles respectively. Note that the states are ordered so that the matrix T Q = {q ij } is upper triangular. The vector of initial state probabilities is noted T P robEnt. The algorithm also uses the vectors T GamStar and T V P ole. The component T GamStar(i) provides the pole number of the state i while the value of this pole is provided by the component T V P ole(i). The algorithm also uses the Boolean matrices BP ole(N S, N P ), BP ar(N S, N S) and BP olP ar(N S, N P ), matrices defined as follows: A second algorithm allowing to fill the matrices T GamStar, T V P ole, BP ole, BP ar and BP olP ar using only the knowledge of the matrix T Q is (fortunately) proposed in table 3.

The matrices BP ole (provided by the algorithm of the table 3), a(i, n, k), and Kmax(i, n) (provided by the algorithm of the tables 1 and 2) are used to compute the values of the various state probabilities P i (t). ---------------------------BEGIN INPUT DATA : N E ; N P ; T Q(i, j) ; T P robEnt(i) ; BP oles(i, n) ; BP ar(j, i) ; T V P ole(n) ;

BP olP ar(i, n) ; T GamStar(i) ; i, j = 1, 2, 3, ..., N E; n = 1, 2, ..., N P . OUTPUT DATA : a(i, n, k), Kmax(i, n) . % Initialization state 1 : a(1,1,0) := TProbEnt(1) Kmax(1,1) := 0 % Loop over the states :

Sum := 0.0 For j = 1 : (i -1) Do V alKmj := Kmax(j, n) If ((BP ar(j, i)) and (BP olP ar(i, n)) and

% Ending of the computation of coefficients ã(i, n, k) relative to the node i. --------------------------- ----------------------------% Start computation of a(i, n, k) relative to the poles of the node i exept pole γ * GS := T GamStar(i) ---------------------------