Two approximation formulas for the mean waiting time of a system with multi-type jobs and multi-type servers

Bachelor Final Project

Abstract 3

1 Introduction 4

2 A system with two servers and two job types 6

2.1 Model description and product-form stationary distribution 6

2.2 Waiting time distributions and approximation formulas 7

2.2.1 Approximation based on weighted average residual service 8

2.2.2 Approximation based on light traffic-heavy traffic interpolation 10

3 The Multidimensional Model 13

3.1 A system with three servers and three job types 13

3.2 Main result 16

3.3 Quality assessment 18

References 19

Abstract

In this report we continue studying the model described by Visschers, Adan and Weiss [1]. They considered a memoryless single station service system with servers M={'〖' m'〗'_1,…,m_K}, and with job types C={a,b,c,…}. Service is skill-based, so that server m_i can serve a subset of job types C(m_i ). Waiting jobs are served on a first-come-first-served basis, while arriving jobs that find multiple idle servers are assigned to a feasible server randomly. Visschers, Adan and Weiss showed the existence of assignment probability distributions under which the system has a product-form stationary distribution, and obtained explicit expressions for it. They also derived waiting time distributions in steady state. Using these results, we derive two approximation formulas for the mean waiting time for general service requirements, one based on the weighted average residual service and one based on light traffic-heavy traffic interpolation, and assess the quality using simulation.

Introduction

We continue studying the model used by Visschers, Adan and Weiss [1]. They studied a service system with K servers (machines), labelled'〖' m'〗'_1,…,m_K, which serves several types of jobs, labelled a,b,c,…. The set of job types is denoted by C and the set of machines by M. Jobs arrive at the system in independent Poisson streams with rates λ_i,i'∈'C, and have independent service requirements. Each machine is capable of handling a specific subset of job types. Machine m_i can only handle jobs from the set C(m_i )'⊂'C and works at rate μ_(m_i ).

The service discipline is a combination of First-Come-First-Served (FCFS) and random assignment. Arriving jobs which find no feasible available machine wait in a single queue, and are processed in an FCFS order as long as it is possible. This means that as soon as a machine finishes a job it takes the first job in the queue that it can process, possibly skipping several jobs that it cannot process. Jobs which upon arrival find multiple available feasible machines are assigned to one of them randomly and go into service immediately.

We assume that an arriving job will choose a feasible machine from the set of idle machines according to a specified assignment probability distribution which depends on the job type and on the set of idle machines. There is one assignment probability distribution for each type of job and for each subset of idle machines which contains at least one feasible machine for that type of job. We treat these assignment probability distributions as control parameters for the system. Figure 1 shows two examples of such systems, which we will refer to as System A and System B.

Fig. 1 Some service systems with multi-type jobs and multi-type servers

To fully describe the system, we list all the jobs in the system in their order of arrival, including jobs which are being processed, and imagine that the machines are situated in the queue on the position of the job that they are processing. To illustrate we consider System B in Fig. 1, in which there are three job types, a,b,c, and three machines, m_1,m_2,m_3, with C(m_1 )={a,b},C(m_2 )={a,c} "and" C(m_3 )={a}.

Fig. 2 depicts a possible state of System B. The jobs are denoted by circles and the machines by rectangles. Jobs in service have a rectangle drawn around them with the identity of the machine inside. There are 12 jobs in the system and all three machines are busy. Machine m_1 is processing the first job in line, which must therefore be either of types a or b. Following the line, machine m_2 is processing the first job in the line which it can process, which is job 5 in the line, and must therefore be either of types a or c. Machine m_3 is processing the first job in the line (apart from jobs 1 and 5) which it can process, which must be of type a. There are three jobs waiting between machines m_1 and m_2. These cannot be processed by either machine m_2 or by machine m_3, so they must be of type b. There are four jobs waiting between machines m_2 and m_3. Those cannot be processed by machine m_3, so they must be of types b or c. Finally there are two jobs at the back of the queue, behind machine m_3, which may be either of types a,b, or c.

Fig. 2 A possible state for System B

We can aggregate some of the states in this description to simplify the model. We will retain the identity and location of the busy machines, but we will not specify the type of job which they are working on, and we will only record the number of jobs between the busy machines, and not specify the string of job types. Thus the state of Fig. 2 will be denoted as (2,m_3,4,m_2,3,m_1). With this reduced description the system is still Markovian.

Visschers, Adan and Weiss showed that, for exponentially distributed service requirements, there exist choices of the assignment probability distributions, such that the system has a product-form stationary distribution. This product-form stationary distribution is unique, and we obtain it explicitly. Furthermore, for exponentially distributed service requirements, they derived the waiting time distributions of jobs of various types.

In this report we derive two approximation formulas for the mean waiting time of jobs of various types for general service requirements, one based on the weighted average residual service and one based on light traffic-heavy traffic interpolation. These formulas are based on the results obtained by Visschers, Adan and Weiss [1] and are given in Section 3.2.

To motivate and illustrate both formulas we will first analyse System A of Fig. 1, with two job types, C={a,b}, and two servers, M={m_1,m_2 }, where C(m_1 )={a,b} and C(m_2 )={a}. Here only jobs of type a have a choice of machines, and that happens only when a job of type a arrives to find both machines idle; in other words, when the system is empty. Hence the assignment probability distributions reduce simply to the probability η of assigning an arriving job of type a to machine m_1 when the system is empty. In Section 2 we give the correct assignment probability η and the product-form solution derived by Visschers, Adan and Weiss and summarized in Theorem 1 of [1] (Sect. 2.1), and derive the approximation formulas for the mean waiting time for general service requirements based on these results (Sect. 2.2).

We then formulate and derive our main result in Section 3. We will first derive both formulas for System B of Fig. 1 (Sect. 3.1), then for the multidimensional model (Sect. 3.2), and conclude by assessing the quality using simulation (Sect. 3.3).

Note that the model used in this report is formulated as a manufacturing model. However, it should find as much use also to describe, for example, skill-based routing of calls to operators in call centres, routing of wireless messages to ad hoc nodes, processing of chips, multiprocessor scheduling and mounting of printed circuit boards; see in particular [2-5].

A system with two servers and two job types

In this section we analyse System A of Fig. 1. Adan, Foley and McDonald [6] found the product-form solution for this system for a special value of η. Visschers, Adan and Weiss [1] derived the waiting time distributions in steady state. Here we restate the product-form solution and the waiting time distributions for this particular system and develop an intuition for the approximation formulas, to illustrate and motivate the general result in Sect. 3.

Note that service times will be exponential at first, before moving to general service requirements. For ease of presentation, the notation used in this section will slightly deviate from the one introduced in Section 1.

Model description and product-form stationary distribution

System A has two machines, m_1,m_2, and two job types, a,b, where machine m_1 can serve both types, and machine m_2 can serve only type-a jobs. Arrivals are Poisson, with rates λ_a,λ_b, and service times are exponential with machine-dependent rates μ_i,i=1,2. Service is on an FCFS basis, and there is one control parameter:

η is the probability that a type-a job that arrives to find an empty system will be assigned to machine m_1.

Clearly, for ergodicity it is necessary that both

λ_b/μ_1 <1, (λ_a+λ_b)/(μ_1+μ_2 )<1

(1)

should hold. Visschers, Adan and Weiss showed that these conditions are also sufficient.

In our Markovian description the system can be in the following states:

(n,m_2,m,m_1), where m,n≥0, with m+n+2 jobs, machine m_1 is working on the first job followed by m jobs waiting between machine m_1 and machine m_2, all of which must be of type b, followed by machine m_2 working on the m+2nd job in line, followed by additional n jobs that have not yet been identified, and are waiting behind machine m_2.

(n,m_1,0,m_2), where n≥0, with n+2 jobs, machine m_2 is working on the first job, machine m_1 is working on the second job, and there are additional n jobs of unidentified type waiting behind the two machines.

(m,m_1), where m≥0, machine m_1 is working on the first job, followed by m jobs which are all of type b, and machine m_2 is idle.

(0,m_2), machine m_2 is working on a single job in the system, and machine m_1 is idle.

(0), the empty system state.

Arrivals join the queue or activate an available machine. On completion of service by machine m_1, if the queue is not empty, the earliest job goes into service. On completion of service by machine m_2, the machine moves along the line of n unidentified jobs until it finds a type-a job to process, or if none is available, it becomes idle. The jobs which are passed over in this search are now identified as type b and are queued between the two machines. Each of the unidentified n jobs at the end of the line is of type b with a probability γ=λ_b/(λ_a+λ_b ). Thus machine m_2 will identify j-1 jobs of type b and then find a type-a job with probability p_j=(1-γ) γ^(j-1),j=1,…,n, or identify all n jobs as type b and become idle with probability γ^n.

We can now derive the equilibrium equations for this system. Visschers, Adan and Weiss solved these equations for the stationary probabilities, denoted by π(x) for state (x), by assuming that the solution is of product-form, and found the unique choice of η for which the product-form solution exists. The results are summarized in the following theorem [1]:

Theorem 1 The system is ergodic if and only if λ_b/μ_1 <1, (λ_a+λ_b)/(μ_1+μ_2 )<1. The system has a product-form solution if and only if η=λ_a/(2λ_a+λ_b ). The stationary distribution in that case is

π(n,m_2,m,m_1 )=B(λ_b/μ_1 )^m ('〖'λ_a+λ'〗'_b/'〖'μ_1+μ'〗'_2 )^n,m,n≥0,

π(n,m_1,0,m_2 )=B μ_1/μ_2 '〖' ('〖'λ_a+λ'〗'_b/'〖'μ_1+μ'〗'_2 )'〗'^n,n≥0,

π(m,m_1 )=B (μ_1+μ_2)/λ_a (λ_b/μ_1 )^m,m≥0,

π(0,m_2 )=B μ_1/μ_2 (μ_1+μ_2)/(λ_a+λ_b ),

π(0)=B (2λ_a+λ_b)/(λ_a+λ_b ) μ_1/λ_a (μ_1+μ_2)/(λ_a+λ_b ),

where the normalizing constant is given by

B=(μ_2 λ_a (λ_a+λ_b )^2 (μ_1-λ_b )(μ_1+μ_2-'〖'λ_a-λ'〗'_b))/(μ_1 (μ_1+μ_2)(μ_2^2 λ_a^2+μ_1 λ_a (μ_1-λ_b )(λ_a+λ_b )+μ_1 μ_2 (2λ_a+λ_b)(μ_1+μ_2-λ_b))) .

The product-form solution satisfies partial balance: flux into a state due to arrivals equals flux out of that state due to departures. If one expands the state space to include the identity of the job in service by machine m_1, then there is no product-form solution.

Waiting time distributions and approximation formulas

We will now try to develop an intuition for what the approximation formulas should be. Consider the system is in steady state. Jobs of both types a and b arrive as a Poisson stream, and hence they see the queue in steady state, and find machines m_1,m_2 busy with probability π('⋅',m_1,'⋅',m_2). If at least one of the machines is idle, jobs of type a will go into service immediately, and the waiting time will be zero. If machine m_1 is idle, regardless the state machine m_2 is in, jobs of type b will also go into service immediately. Else, they will have to wait a sum of exponential waiting times.

Type-a jobs wait only when both machines are busy, which occurs with probability ∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) +∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) . The waiting time is exponentially distributed with rate μ_1+μ_2-λ_a-λ_b, which is the waiting time in an M/M/1 queue with arrival rate λ_a+λ_b and service rate μ_1+μ_2. Then the Laplace-Stieltjes transform (LST) of the steady-state waiting time W_a of a job of type a, as derived by Visschers, Adan and Weiss, is equal to

E(e^(-s W_a ) )=(π(0)+π(0,m_2 )+∑_(m=0)^∞'▒'π(m,m_1 ) )'⋅'1

+(∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) +∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) )'⋅'( μ_1+μ_2-λ_a-λ_b)/( μ_1+μ_2-λ_a-λ_b+s) ,

(2)

where π('⋅',m_1,'⋅',m_2) is given by Theorem 1.

Things become more difficult for jobs of type b. Note that type-b jobs will have to wait if machine m_1 is busy, regardless the state machine m_2 is in. We will treat all three cases individually. When only machine m_1 is busy, and machine m_2 is idle, which occurs with probability ∑_(m=0)^∞'▒'π(m,m_1 ) , the waiting time is exponentially distributed with rate μ_1-λ_b, which again is the waiting time in an M/M/1 queue, but now with arrival rate λ_b and service rate μ_1. When the system is in the state ('⋅',m_1,0,m_2 ), with probability ∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) , the arriving job will have to wait for the n jobs of unidentified type waiting behind the two machines. The waiting time in that case is again exponentially distributed with rate μ_1+μ_2-λ_a-λ_b. Finally, when the system is in the state ('⋅',m_2,'⋅',m_1 ), with probability ∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) , the arriving job will have to wait for the m jobs of type b waiting between machine m_1 and machine m_2 and for the n jobs of unidentified type waiting behind the two machines, hence it will have to wait a sum of exponential waiting times, with rates μ_1-λ_b and μ_1+μ_2-λ_a-λ_b.

The LST of the steady-state waiting time W_b of a job of type b is then equal to

E(e^(-s W_b ) )=(π(0)+π(0,m_2 ))'⋅'1

+∑_(m=0)^∞'▒'π(m,m_1 ) '⋅'(μ_1-λ_b)/(μ_1-λ_b+s)

+∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) '⋅'( μ_1+μ_2-λ_a-λ_b)/( μ_1+μ_2-λ_a-λ_b+s)

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) '⋅'(μ_1-λ_b)/(μ_1-λ_b+s) ( μ_1+μ_2-λ_a-λ_b)/( μ_1+μ_2-λ_a-λ_b+s) ,

(3)

where π('⋅',m_1,'⋅',m_2) is given by Theorem 1.

Approximation based on weighted average residual service

We will now use the same ideas behind the LSTs to derive the first approximation formula for the mean waiting time. To this end we sum conditional waiting time terms for all states in which the arriving job has to wait multiplied by the steady-state probabilities given by Theorem 1 to uncondition. To illustrate, first, the mean waiting time can be found by using the mean value technique. For type-a jobs, if not both servers are busy on arrival the waiting time is zero. If both servers are busy and there are zero or more jobs waiting, then an arriving job first has to wait until the first departure and then continues to wait for as many departures as there were jobs waiting upon arrival. An inter-departure time is the minimum of two exponential (residual) service times with means 1/μ_1 and 1/μ_2, and thus it is exponential with mean 1/(μ_1+μ_2). So we obtain

E(W_a )=Π_(W_a )'⋅'1/(μ_1+μ_2 )+E(L^q )'⋅'1/(μ_1+μ_2 ) ,

(4)

where'〖' Π'〗'_(W_a ) is the probability that a type-a job has to wait. From Little’s law we get

E(W_a )=Π_(W_a )'⋅'1/(1-ρ)'⋅'1/(μ_1+μ_2 ),"where" ρ=(λ_a+λ_b)/(μ_1+μ_2 ) .

(5)

We will now rewrite the 1/(μ_1+μ_2) term in terms of the mean residual service times. With probability μ_1/(μ_1+μ_2) machine m_1 will have completed its current task first, and hence will be able to serve the arriving job. In that case the waiting time is one residual service time with mean 1/μ_1. With probability μ_2/(μ_1+μ_2) machine m_2 will have completed its current task first. The arriving job then waits a residual service time with mean 1/μ_2. To find the mean, we just multiply by ½. We get

1/(μ_1+μ_2 )=1/2 ( μ_1/(μ_1+μ_2 ) 1/μ_1 +μ_2/(μ_1+μ_2 ) 1/μ_2 ) .

(6)

Hence the desired approximation formula for the mean waiting time for type-a jobs for exponentially distributed service requirements is

E(W_a )=Π_(W_a )'⋅'1/(1-ρ)'⋅'1/2 ( μ_1/(μ_1+μ_2 ) 1/μ_1 +μ_2/(μ_1+μ_2 ) 1/μ_2 ) ,

(7)

where'〖' Π'〗'_(W_a )=∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) +∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) .

Note that from (2) we get exactly the same without using the mean value technique. For type-b jobs things are again somewhat more difficult. When only machine m_1 is busy, and machine m_2 is idle, the arriving job waits only for all type-b jobs already in line. By using the mean value technique and applying Little’s law, we get a mean waiting time of 1/(1-ρ_1)'⋅'1/μ_1 where ρ_1=λ_b/μ_1 . When the system is in the state ('⋅',m_1,0,m_2 ), it waits for all jobs of unidentified type waiting behind the two machines, hence the mean waiting time is 1/(1-ρ)'⋅'1/(μ_1+μ_2). Finally, when the system is in the state ('⋅',m_2,'⋅',m_1 ), it again waits a sum of exponential waiting times, hence we get a mean waiting time of 1/(1-ρ_1 )'⋅'1/μ_1 + 1/(1-ρ)'⋅'1/(μ_1+μ_2). In summary, for type-b jobs the approximation formula for the mean waiting time for exponentially distributed service requirements is

E(W_b )=∑_(m=0)^∞'▒'π(m,m_1 ) '⋅'1/(1-ρ_1 )'⋅'1/μ_1

+∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) '⋅'1/(1-ρ)'⋅'1/2 ( μ_1/(μ_1+μ_2 ) 1/μ_1 +μ_2/(μ_1+μ_2 ) 1/μ_2 )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) '⋅'(1/(1-ρ_1 )'⋅'1/μ_1 +1/(1-ρ)'⋅'1/2 ( μ_1/(μ_1+μ_2 ) 1/μ_1 +μ_2/(μ_1+μ_2 ) 1/μ_2 )) ,

(8)

where ρ=(λ_a+λ_b)/(μ_1+μ_2 ) and ρ_1=λ_b/μ_1 . Note again that (3) gives exactly the same formula without using the mean value technique.

For general service requirements we now substitute the residual waiting times 1/μ_i with the general expressions E(R_i),i=1,2. For μ_i we write 1/E(B_i) ,i=1,2, the reciprocal of the mean service time. Note that this also changes the stationary distribution given by Theorem 1. Hence we simply get

E(W_a )=(∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) +∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) )

⋅'1/(1-ρ)'⋅'1/2 ( μ_1/(μ_1+μ_2 ) E(R_1 )+μ_2/(μ_1+μ_2 ) E(R_2 )),

(9)

and

E(W_b )=∑_(m=0)^∞'▒'π(m,m_1 ) '⋅'1/(1-ρ_1 )'⋅'E(R_1)

+∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) '⋅'1/(1-ρ)'⋅'1/2 ( μ_1/(μ_1+μ_2 ) E(R_1)+μ_2/(μ_1+μ_2 ) E(R_2) )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) '⋅'(1/(1-ρ_1 )'⋅'E(R_1 )+1/(1-ρ)'⋅'1/2 ( μ_1/(μ_1+μ_2 ) E(R_1 )+μ_2/(μ_1+μ_2 ) E(R_2 ))) ,

(10)

where ρ=(λ_a+λ_b)/(μ_1+μ_2 ) , ρ_1=λ_b/μ_1 , μ_i=1/E(B_i) "for" i=1,2 and the stationary distribution π( '⋅' ) is given by Theorem 1.

Note that this approximation formula assumes that the stationary distribution given by Theorem 1 is insensitive, in that the stationary distribution depends on the service time distributions only through their means. This result however was proven by Adan, Hurkens and Weiss for a closely related system [7].

Approximation based on light traffic-heavy traffic interpolation

Next we will derive the second approximation formula for the mean waiting time which is based on interpolation between light traffic and heavy traffic behaviour. For this we use the work done by Reiman and Simon [8], who proposed the use of a Taylor expansion using both light traffic and heavy traffic information to obtain an approximation for all traffic intensities.

In heavy traffic we take the limit of the normalized waiting time process as the arrival rate, λ, approaches its critical value, ζ. In the case of System A we have λ=λ_a+λ_b and ζ=μ_1+μ_2. We will show that in the limit the system behaves as though it is an M/G/1 queue, hence there is just one limit for both job types.

In light traffic the unnormalized limit is trivial, since in the limit there are no customers in the system. We can obtain more sensitive information on light traffic behaviour by calculating derivatives with respect to λ at λ=0. In theory, with n derivatives, this approach leads to an nth degree polynomial in λ. Here we only consider second degree polynomials.

Using light traffic information alone leads to an approximation that may be quite poor in moderate and heavy traffic. For example, the average waiting time function, W(λ), increases to infinity as λ→ζ, so no polynomial will approximate it well. Heavy traffic theory supplies valuable information by showing that lim'┬'(λ→ζ)''〖' (ζ-λ)'⋅'W(λ)'〗' has a finite nonzero limit, and offering a simple formula for the limit. Therefore, if we “normalize” W(λ) by (ζ-λ), and consider the function W ̂(λ)=(ζ-λ) W(λ), a polynomial will be a much more reasonable approximation. In the remainder of this section we will obtain the polynomial approximation for System A utilizing heavy-traffic limits, light-traffic limits, and light-traffic derivatives. In Sect. 3.2 we will give the general result.

We will first consider the behaviour of System A in heavy traffic, that is, as the arrival rate, λ, converges to its critical value, ζ. Intuitively in that case the system will behave as though it is an M/G/1 queue and hence the waiting time is just the mean residual service time.

This also follows from the third form of the Pollaczek-Khinchine formula. Let f[s] be the Laplace-Stieltjes transform of the service time distribution of both machines combined (for exponential service times this is just (μ_1+μ_2)/(μ_1+μ_2+s) ). Note that E(B)=-1/(f^' [0] ) and ρ=-λ f^' [0]. Then

〖'W ̂(ζ)=lim'┬'(λ→ζ)'〗''〖' (ζ-λ)'⋅'W(λ)'〗'=lim'┬'(ρ→1)''((1-ρ) lim'┬'(s→0)''(-d/ds (((1-ρ) s)/(λ f[s]+s-λ))) )

〖'=lim'┬'(ρ→1)'〗''((1-ρ) lim'┬'(s→0)''( (1-ρ)(λ-λ f[s]+s λ f^' [s])/(λ f[s]+s-λ)^2 ) )

〖'=lim'┬'(ρ→1)'〗''((1-ρ) lim'┬'(s→0)''( (d^2/(ds^2 ) ((1-ρ)(λ-λ f[s]+s λ f^' [s])))/(d^2/(ds^2 ) ((λ f[s]+s-λ)^2 ) ) ) )

〖'=lim'┬'(ρ→1)'〗''((1-ρ)'⋅'(λ (1-ρ) f^'' [0])/(2 (1+λ f^' [0])^2 ))

〖'=lim'┬'(λ→-1/(f^' [0] ))'〗''(1/2 λ f^'' [0])

=-(f^'' [0])/(2 f^' [0] )

=E(B^2 )/(2 E(B) ) ,

(11)

and hence W ̂(ζ)=lim'┬'(λ→ζ)''〖' (ζ-λ)'⋅'W(λ)'〗'=E(R).

Calculation of f[s] for more complex systems can be quite challenging, however the same limit can be obtained from the first approximation formula. For type-a jobs we have

〖'(W_a ) ̂(ζ)=lim'┬'(ρ→1)'〗''(1-ρ)E(W_a )=lim'┬'(ρ→1)''(1-ρ) (Π_(W_a )'⋅'E(R)+E(L^q )'⋅'E(R))

=lim'┬'(ρ→1)''(1-ρ) (Π_(W_a )'⋅'1/(1-ρ)'⋅'E(R))

=lim'┬'(ρ→1)''〖'Π_(W_a ) '〗'⋅'E(R)

=E(R).

(12)

Now we rewrite E(R) the same way we did previously, hence the heavy-traffic limit is

〖'(W_a ) ̂(ζ)=lim'┬'(ρ→1)'〗''(1-ρ)E(W_a )=E(R)

=1/2 ( μ_1/(μ_1+μ_2 ) E(R_1 )+μ_2/(μ_1+μ_2 ) E(R_2 ))

=1/2 ( μ_1/(μ_1+μ_2 ) E('〖'B_1'〗'^2 )/(2 E(B_1 ) )+μ_2/(μ_1+μ_2 ) E('〖'B_2'〗'^2 )/(2 E(B_2 ) ) ) ,

(13)

where μ_1=1∕E(B_1 ) and μ_2=1∕E(B_2 ) . It can also be proven that (1-ρ)E(W_b ) has the same limit, hence there is just one heavy-traffic limit for both job types.

Now we examine the behaviour of System A in light traffic, that is for small values of λ. The unnormalized limit is trivial, since in the limit there are no customers in the system. To obtain more sensitive information, we calculate derivatives of W ̂(λ)=(ζ-λ) W(λ) with respect to λ at λ=0. More specifically we calculate the first nonzero derivative. For W(λ) we choose to use the first approximation formula. Note that now we have to make a difference between the two job types, hence we will use E(W_a ) and E(W_b ) instead of W(λ). Also note that the conditional waiting times E[W_i |'〖' Π'〗'_(W_i )>0] do not contain any λ. Only the steady-state probabilities Π_(W_i ) contain λ, and, as Adan, Hurkens and Weiss showed, these are insensitive to the service time distributions, which justifies our choice. Hence in light traffic we have

(W_i ) ̂^((m)) (0)='├' d^m/(dλ^m ) ((ζ-λ) W_i (λ))'┤'|_(λ=0),i=a,b,

(14)

where W_i (λ)=E(W_i ) are given by (9) and (10) and m the smallest integer such that (W_i ) ̂^((m)) (0)>0.

We now combine these results to obtain the second approximation formula using Taylor expansion. Assume that we have found W ̂(0) as well as the first nonzero derivatives of (W_i ) ̂(λ). In addition, assume that we have determined W ̂(ζ). There is then a unique second degree polynomial, F(λ),i=a,b, such that

F(0)=W ̂(0),

F^' (0)=(W_i ) ̂^((m) ) (0),

"and"

F(ζ)=(W_i ) ̂(ζ).

(15)

Let F(λ) be of the form A+B ρ+C ρ^2, then

A=W ̂(0),

B=(W_i ) ̂^((m) ) (0),

A+B+C=(W_i ) ̂(ζ),

(16)

and hence

F_i (λ)=W ̂(0)+(W_i ) ̂^((m) ) (0) ρ+((W_i ) ̂(ζ)-W ̂(0)-(W_i ) ̂^((m) ) (0)) ρ^2,i=a,b.

(17)

The mean waiting times for type¬-a and type-b jobs then are

E(W_a )=(W ̂(0)+(W_a ) ̂^((m) ) (0) ρ+((W_a ) ̂(ζ)-W ̂(0)-(W_a ) ̂^((m) ) (0)) ρ^2)/(1-ρ),

(18)

and

E(W_b )=(W ̂(0)+(W_b ) ̂^((m) ) (0) ρ+((W_b ) ̂(ζ)-W ̂(0)-(W_b ) ̂^((m) ) (0)) ρ^2)/(1-ρ).

(19)

The Multidimensional Model

In this section we consider the general K machine system. To demonstrate we will first analyse System B of Fig. 1. We will then give the general result in Section 3.2. Finally, we assess the quality of the approximation formulas using simulation.

We introduce the following notation:

λ_X '≔' ∑_(c'∈'X)'▒'λ_c , where X'⊂'C.

μ_Y '≔' ∑_(m_i'∈'Y)'▒'μ_(m_i ) , where Y'⊂'M.

C(Y) '≔' total set of job types that can be handled by the machines in Y'⊂'M, which is equal to '⋃'_(m_i'∈'Y)'▒'〖'C(m_i)'〗'.

U(Y) '≔' set of job types unique to the machines in Y'⊂'M, thus the set of job types that cannot be handled by machines outside Y. We have U(Y)=(C(Y ̅)) ̅.

A system with three servers and three job types

We will now analyse System B of Fig. 1, in which there are three job types a,b,c and three machines with C(m_1 )={a,b},C(m_2 )={a,c},C(m_3 )={a}. Arrivals are Poisson, with rates λ_i,i=a,b,c, and service times are generally distributed with machine-dependent rates μ_(m_i ),i=1,2,3. Service is on an FCFS basis. For ergodicity it is necessary and sufficient that

(λ_a+λ_b)/μ_(m_1 ) <1, (λ_a+λ_c)/μ_(m_2 ) ,"and" (λ_a+λ_b+λ_c)/(μ_(m_1 )+μ_(m_2 )+μ_(m_3 ) )<1

(20)

should hold. There are several control parameters:

η_(m_j ) (S) is the rate at which server m_j'∈'S becomes busy when the set of idle servers is S.

We can obtain these rates from the formula derived by Adan, Hurkens and Weiss [7]

η_(m_j ) (S)=η(S)'⋅'(1+∑_(m_j'∈'S\\\\{m_k})'▒'(η_(m_j ) (S\\\\{m_k}))/(η_(m_k ) (S\\\\{m_j})))^(-1),"where " η(S)=∑_(m_j'∈'S)'▒'〖'η_(m_j ) (S) '〗'=∑_(c'∈'C(S))'▒'λ_c .

(21)

In our Markovian description the system can be in the following states:

π(0),π(n,m_1 ),π(m,m_2 ),π(0,m_3 ),

π(n,m_1,m,m_2 ),π(n,m_1,0,m_3 ),π(n,m_2,0,m_3 ),

π(n,m_2,m,m_1 ),π(n,m_3,m,m_1 ),π(n,m_3,m,m_2 ),

π(n,m_1,m,m_2,0,m_3 ),π(n,m_2,m,m_1,0,m_3 ),π(n,m_3,m,m_1,p,m_2 ),

π(n,m_1,0,m_3,m,m_2 ),π(n,m_2,0,m_3,m,m_1 ),π(n,m_3,m,m_2,p,m_1 ).

In state (n,m_3,m,m_2,p,m_1 ), where m,n,p≥0, with m+n+p+3 jobs, machine m_1 is working on the first job followed by p jobs waiting between machine m_1 and machine m_2, all of which must be of type b, followed by machine m_2 working on the p+2nd job in line, followed by m jobs waiting between machine m_2 and machine m_3, all of which must be of type b or c, followed by machine m_3 working on the m+p+3rd job in line, followed by additional n jobs that have not yet been identified, and are waiting behind machine m_3.

We will now derive the first approximation formula for the mean waiting time for System B. Let M^n denote the last n servers (in order) in state (x), thus M^1={m_1 },M^2={m_1,m_2 } "and" M^3={m_1,m_2,m_3 } when the system is in the state described above. Let'〖' M'〗'^n (i) denote the'〖' i'〗'^th element of the set M^n, thus M^3 (1)=m_1. We introduce the following notation:

W_X'≔'μ_('〖' M'〗'^j )/(μ_('〖' M'〗'^j )-λ_X )'⋅'1/j'⋅'∑_(i=1)^j'▒'μ_('〖' M'〗'^j (i) )/μ_('〖' M'〗'^j ) E(R_('〖' M'〗'^j (i) ) ),

where j is the number of elements in X. Recall that for System A we had

E(W_a )=(∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) +∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) )

⋅'1/(1-ρ)'⋅'1/2 ( μ_1/(μ_1+μ_2 ) E(R_1 )+μ_2/(μ_1+μ_2 ) E(R_2 )),

and

E(W_b )=∑_(m=0)^∞'▒'π(m,m_1 ) '⋅'1/(1-ρ_1 )'⋅'E(R_1)

+∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) '⋅'1/(1-ρ)'⋅'1/2 ( μ_1/(μ_1+μ_2 ) E(R_1)+μ_2/(μ_1+μ_2 ) E(R_2) )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) '⋅'(1/(1-ρ_1 )'⋅'E(R_1 )+1/(1-ρ)'⋅'1/2 ( μ_1/(μ_1+μ_2 ) E(R_1 )+μ_2/(μ_1+μ_2 ) E(R_2 ))) ,

where ρ=(λ_a+λ_b)/(μ_1+μ_2 ) , ρ_1=λ_b/μ_1 and μ_i=1/E(B_i) "for" i=1,2. We can now write

E(W_a )=∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) '⋅'μ_('〖' M'〗'^2 )/(μ_('〖' M'〗'^2 )-λ_{a,b} )'⋅'1/2'⋅' ( μ_('〖' M'〗'^2 (1) )/μ_('〖' M'〗'^2 ) E(R_('〖' M'〗'^2 (1) ) )+μ_('〖' M'〗'^2 (2) )/μ_('〖' M'〗'^2 ) E(R_('〖' M'〗'^2 (2) ) ))

+∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) '⋅'μ_('〖' M'〗'^2 )/(μ_('〖' M'〗'^2 )-λ_{a,b} )'⋅'1/2'⋅' ( μ_('〖' M'〗'^2 (1) )/μ_('〖' M'〗'^2 ) E(R_('〖' M'〗'^2 (1) ) )+μ_('〖' M'〗'^2 (2) )/μ_('〖' M'〗'^2 ) E(R_('〖' M'〗'^2 (2) ) )),

where'〖' M'〗'^2 are the last two servers in states π(n,m_2,m,m_1 ) and π(n,m_1,0,m_2 ) respectively. Since for type-a jobs the order of machines m_1 and m_2 does not matter, we can write

E(W_a )=(∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) +∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) )'⋅'W_{a,b} .

For type-b jobs we write

E(W_b )=∑_(m=0)^∞'▒'π(m,m_1 ) W_{b} +∑_(n=0)^∞'▒'π(n,m_1,0,m_2 ) W_{a,b} +∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) (W_{b} +W_{a,b} ).

Hence for System B we have

E(W_a )=(∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_1,m,m_2,0,m_3 ) +∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_1,0,m_3,m,m_2 ) '┤

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1,0,m_3 ) +∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,0,m_3,m,m_1 )

├' +∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'∑_(p=0)^∞'▒'π(n,m_3,m,m_1,p,m_2 ) +∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'∑_(p=0)^∞'▒'π(n,m_3,m,m_2,p,m_1 ) )'⋅'W_{a,b,c} , (22)

E(W_b )=∑_(n=0)^∞'▒'π(n,m_1 ) '⋅'W_{b} +∑_(n=0)^∞'▒'π(n,m_1,0,m_3 ) '⋅'W_{a,b} +∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_1,m,m_2 ) '⋅'W_{b,c}

+(∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_1,m,m_2,0,m_3 ) +∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_1,0,m_3,m,m_2 ) )'⋅'W_{a,b,c}

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_3,m,m_1 ) '⋅'('〖'W_{b} +W'〗'_{a,b} )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) '⋅'(W_{b} +W_{b,c} )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1,0,m_3 ) '⋅'(W_{b} +W_{a,b,c} )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,0,m_3,m,m_1 ) '⋅'(W_{b} +W_{a,b,c} )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'∑_(p=0)^∞'▒'π(n,m_3,m,m_1,p,m_2 ) '⋅'(W_{b,c} +W_{a,b,c} )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'∑_(p=0)^∞'▒'π(n,m_3,m,m_2,p,m_1 ) '⋅'(W_{b} +W_{b,c} +W_{a,b,c} ),

(23)

E(W_c )=∑_(n=0)^∞'▒'π(n,m_2 ) '⋅'W_{c} +∑_(n=0)^∞'▒'π(n,m_2,0,m_3 ) '⋅'W_{a,c} +∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1 ) '⋅'W_{b,c}

+(∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,m,m_1,0,m_3 ) +∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_2,0,m_3,m,m_1 ) )'⋅'W_{a,b,c}

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_3,m,m_2 ) '⋅'('〖'W_{c} +W'〗'_{a,c} )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_1,m,m_2 ) '⋅'(W_{c} +W_{b,c} )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_1,m,m_2,0,m_3 ) '⋅'(W_{c} +W_{a,b,c} )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'π(n,m_1,0,m_3,m,m_2 ) '⋅'(W_{c} +W_{a,b,c} )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'∑_(p=0)^∞'▒'π(n,m_3,m,m_2,p,m_1 ) '⋅'(W_{b,c} +W_{a,b,c} )

+∑_(m=0)^∞'▒'∑_(n=0)^∞'▒'∑_(p=0)^∞'▒'π(n,m_3,m,m_1,p,m_2 ) '⋅'(W_{c} +W_{b,c} +W_{a,b,c} ). (24)

The second approximation formula for System B can now be obtained using the first formula as described in the previous section.

Main result

Here we will summarize the results obtained in the previous sections for the general K machine system. Let M={m_1,…,m_K} be the set of machines and C={a,b,c,…} be the set of job types. Arrivals are Poisson, with rates λ_i,i'∈'C. Let C(m_i )'⊂'C, the set of job types that can be handled by machine m_i, be given for all machines m_i'∈'M. Let μ_(m_i )=1∕E(B_(m_i ) ) , the service rate of machine m_i, be known for all machines m_i'∈'M for machine-dependent service time distributions B_(m_i ). Service is on an FCFS basis.

For ergodicity it is both necessary and sufficient that [1]

λ_C(m_i ) /μ_(m_i ) <1 "for all" m_i'∈'M,"and" λ_C/μ_M <1.

(25)

For models with K machines we use the following state description:

(n_i,m_i,n_(i-1),m_(i-1),…,n_1,m_1): States in which there are i machines busy. These machines are denoted by m_1,…,m_i, where {m_1,…,m_i }'⊂'M. The number of jobs between machines m_j and m_(j+1) is denoted by n_j (≥0), with j=1,…,i-1. The number of waiting jobs at the end of the queue, behind machine m_i, is denoted by n_i.

The state space is denoted by S. Let s denote an arbitrary state (n_i,m_i,…,n_1,m_1 )'∈'S. Figure 4 shows a system in such a state.

Fig. 3 General system in state s.

Note that the waiting jobs between machines m_j and m_(j+1) can only be handled by the machines m_1,…,m_j, and not by any of the machines m_(j+1),…,m_i or any of the idle machines, due to the First-Come-First-Served processing order. Thus waiting jobs between machines m_j and m_(j+1) can only be of type c'∈'U({m_1,…,m_j}). The n_i waiting jobs in the back of the queue cannot be handled by any of the idle machines and have to be of type c'∈'U({m_1,…,m_i}).

Let π(s) denote the stationary probability for state s. Let S'⊂'M be the set of idle servers at a given moment. Note that S contains exactly those servers that are not in s. The control parameters for the system are η_(m_j ) (S), the rate at which server m_j'∈'S becomes busy when the set of idle servers is S. These rates can be obtained from [7]

η_(m_j ) (S)=η(S)'⋅'(1+∑_(m_j'∈'S\\\\{m_k})'▒'(η_(m_j ) (S\\\\{m_k}))/(η_(m_k ) (S\\\\{m_j})))^(-1),"where " η(S)=∑_(m_j'∈'S)'▒'〖'η_(m_j ) (S) '〗'=∑_(c'∈'C(S))'▒'λ_c .

(26)

After normalization the stationary probability π(s) for state s=(n_i,m_i,…,n_1,m_1 )'∈'S can be obtained from [1]

π(s)=α_i^(n_i )…α_1^(n_1 ) (η_(m_1 ) (S)…η_(m_i ) (S\\\\{m_1,…,m_(i-1)}))/(μ_(m_1 )…μ_('〖'{m'〗'_1,…,m_i}) ) π(0),

(27)

where

α_j=λ_U({m_1,…,m_j }) /μ_('〖'{m'〗'_1,…,m_j}) ,j=1,2,…,i .

(28)

Let M_s^n denote the last n servers (in order) in state s, thus M_s^1={m_1 },…,M_s^i={m_1,…,m_i }. Let'〖' M'〗'_s^n (i) denote the'〖' i'〗'^th element of the set'〖' M'〗'_s^n, thus'〖' M'〗'_s^i (i)=m_i.

The first approximation formula for the mean waiting time for a job of type c'∈'C is

E(W_c )=∑_(s'∈'S)'▒'〖'π(s) ∑_(c'∈'C(s))'▒'∑_'█'([email protected]'∈'M_s^j )^i'▒'〖'μ_(M_s^j )/(μ_(M_s^j )-λ_U(M_s^j ) )'⋅'1/j'⋅'∑_(k=1)^j'▒'μ_(M_s^j (k) )/μ_(M_s^j ) E(R_(M_s^j (k) ) ) '〗' '〗'.

(29)

Let W_c (λ) denote the mean waiting time for a job of type c'∈'C obtained from (29). Let ζ denote the critical value of the arrival rate λ_C. Consider the function (W_c ) ̂(λ)=(ζ-λ) W_c (λ). Then let

(W_c ) ̂(0)=lim'┬'(λ→0)''〖'(W_c ) ̂(λ)'〗'=lim'┬'(λ→0)''〖' (ζ-λ) W_c (λ)'〗'=0,

(W_c ) ̂^((m) ) (0)=lim'┬'(λ→0)''〖'(W_c ) ̂^((m) ) (λ)'〗'=lim'┬'(λ→0)''〖'd^m/(dλ^m ) ((ζ-λ) W_c (λ))'〗',

(W_c ) ̂(ζ)=lim'┬'(λ→ζ)''〖'(W_c ) ̂(λ)'〗'=lim'┬'(λ→ζ)''〖' (ζ-λ) W_c (λ)'〗'.

(30)

The second approximation formula for the mean waiting time for a job of type c'∈'C is

E(W_c )=(W ̂(0)+(W_c ) ̂^((m) ) (0) ρ+((W_c ) ̂(ζ)-W ̂(0)-(W_c ) ̂^((m) ) (0)) ρ^2)/(1-ρ),

(31)

where ρ=λ_C/μ_M.

Quality assessment

In this last section we will assess the quality of the approximation formulas given by (29) and (31) in Section 3.2 using simulation. We will consider only System A and System B discussed in the report, since the computation of the control parameters and the stationary probabilities needs to be done manually and increases in difficulty with the number of machines.

### References

[1] Visschers, J., Adan, I., Weiss, G.: A product from solution to a system with multi-type jobs and multi-type servers. Queueing Syst, 70, 269-298 (2012)

[2] Adan, I., Wessels, J., Zijm, W.: Queueing analysis in a flexible assembly system with a job-dependent parallel structure. In: Operations Research Proceedings, 1988, pp. 551-558. Springer, Berlin (1989)

[3] Aerts, J., Korst, J., Verhaegh, W.: Load balancing for redundant storage strategies: Multiprocessor scheduling with machine eligibility. J. Sched. 4, 245-257 (2007)

[4] Aksin, Z., Armony, M., Mehrotra, V.: The modern call-center, a multi-disciplinary perspective on operations management research. Prod. Oper. Manag. 16, 665-688 (2007)

[5] Garnett, O., Mandelbaum, A.: An introduction to skill-based routing and its operational complexities (2000). http://iew3.technion.ac.il/serveng/Lectures/SBR.pdf

[6] Adan, I., Foley, R., McDonald, D.: Exact asymptotics of the stationary distribution of a Markov chain: a production model. Queueing Syst. 62, 311-344 (2009)

[7] Adan, I., Hurkens, C., Weiss, G.: A reversible multi-class multi-server loss system. Probability in the Engineering and Informational Sciences. 24, 535-548 (2010)

[8] Reiman, M., Simon, B.: An interpolation approximation for queueing systems with poisson input. Operations Research. 36 (3), 454-469 (1988)

**...(download the rest of the essay above)**