Analysis of an infinite multi-server queue with an optional service

Analysis of an infinite multi-server queue with an optional service

Jau-Chuan Ke

a,⇑

, Chia-Huang Wu

b

, Wen Lea Pearn

b

a

Department of Applied Statistics, National Taichung University of Science and Technology, Taichung, Taiwan, ROC b

Department of Industrial Engineering and Management, National Chiao Tung University, Taiwan, ROC

a r t i c l e

i n f o

Article history: Received 19 June 2012

Received in revised form 18 February 2013 Accepted 22 February 2013

Available online 13 March 2013 Keywords:

First essential channel Linear progressing algorithm Particle Swarm Optimization Quasi-Newton method Second optional channel

a b s t r a c t

This paper deals with an infinite-capacity multi-server queueing system with a second optional service (SOS) channel. The inter-arrival times of arriving customers, the service times of the first essent ial service (FES) and the SOS channel are all exponentially distributed. A customer may leave the system after the FES channel with a probability (1  h), or the completion of the FES may immediately require a SOS with a probability h (0 6 h 6 1). The formulae for computing the rate matrix and stationary probabilities are derived by means of a matrix analytical approach. A cost model is developed to simultaneously deter- mine the optimal values of the number of servers and the two service rates at the minimal total expected cost per unit time. Quasi-Newton method and Particle Swarm Optimization (PSO) method are employed to deal with the optimization problem. Under opt imal operating condit ions, numerical results are provided from which several system performance measure s are calcu lated based on the assumed numerical values of the system parameters.

Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction

In day to day life, one encounter s numerous examples of queue- ing models through which all arriving customer s need an essential service but only some require an additional optional service. In this paper, the quasi-birth–death (QBD) process and matrix analytic methods are used to analyze an M/M/ R queue with a second op- tional service channel. An algorithm is developed to calculate the rate matrix and the stationar y probabilities of the QBD process. A cost model is also constructed to search the optimal parameters at a minimum cost.

A possible application of our model is in a manufacturing sys-

tem for a pump based on the work of Yang, Wang, and Kuo

(2011). Consider a pump manufactur ing industry that manufac- tures different kinds of pumps which require shafts of various dimensions . The arrival of shafts from the turning center to the CNC (computer numerical control) copy turning center follows a random process, in which the center owns multiple CNC machines. The mechanics set up the template in these CNC machines to per- form the copy turning process for shafts, i.e., the first essential ser- vice. The good quality items are kept in storage and are sold. Some of the served (processed) shafts are defective. The defective ones need to be reworked (re-served) to meet the required

specifica-tions. In this scenario, the mechanics (including the CNC-ma- chines) and the re-served action of defective items correspond to the channels and the second optional service, respectively.

It is assumed that customers arrive according to a Poisson pro- cess with parameter k. Customers arriving at the system form a single waiting line and are served in the order of their arrival, that is, first-come-first-served. There are R channels (servers) that pro- vide the first essential service (FES) as well as the second optional service (SOS) to arriving customers. The FES is needed by all arriv- ing customers. The service times of FES and the SOS are exponen- tially distributed with means 1/

l

1and 1/

l

2, respectively . As soon

as the FES of a customer is complete d, they may leave the system with a probability (1  h) or, opt for the SOS provided by the same server with a probability h (0 6 h 6 1). Each channel can serve only one customer at a time and it also provides only either FES or SOS at one time. Customers who, upon entry into the channel facility, find that all channels are busy are required to wait in the queue un- til a channel becomes available. The various stochastic (arrival or service) processes involved in the system are independent of each other.

Analytic steady-state solutions of an M/M/ R queue with a sec- ond optional service channel have not been found. A pioneering

work in this queueing situation was proposed byMadan (2000),

who first introduced the concept of the SOS.Madan (2000)studied

an M/G/1 queue with SOS using the supplementar y variable tech- nique in which he considered a general service time distribution for the FES and an exponential service time distribution for the SOS. He also cited some important applications of this model in

many real-life situation s. Medhi (2002) derived the transient

0360-8352/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved.

http://dx.doi.org/10.1016/j.cie.2013.02.017

⇑Corresponding author. Address: Department of Applied Statistics, National Taichung University of Science and Technology, No. 129, Sec. 3, Sanmin Rd., Taichung 404, Taiwan, ROC. Tel.: +886 422196077; fax: +886 422196331.

E-mail addresses:[email protected](J.-C. Ke),[email protected]

(C.-H. Wu),[email protected](W.L. Pearn).

Contents lists available atSciVerse ScienceDi rect

Computers & Industrial Engin eering

solution and steady-state solution for the ordinary M/G/1 queue with the SOS using the same technique. Medhi’s M/G/1 model

was also investigated byAl-Jararh a and Madan (2003), in which

they developed the time-dependent probability generating func- tions involved in Laplace transforms and further obtained the cor-

responding steady-state results.Choudhury and Madan (2005) and

Choudhury and Paul (2006)studied the queue size distribution at a

random epoch as well as at a departure epoch for an M[x]_/G/1

queueing system with a second optional channel and different con- siderations under N-policy. They also derived a simple procedure to obtain the optimal stationary policy under a suitable linear cost

structure.Tadj, Choudhury, and Tadj (2006a, 2006b)investigated

some bulk service queueing systems under N policy. The reliability

measures were examined byWang (2004)regarding the ordinary

M/G/1 queue with channel breakdow ns and SOS. Recently, Ke

(2008)investigated a batch arrival M[x]/G/1 queueing system with J optional services. He derived the steady-state results, including the system size distribution at a random epoch and at a departure epoch, the distribution s of idle and busy periods, and the waiting

time distribution in the queue.Choudhury and Tadj (2009)

gener-alized this type of model by introducing the concept of a server

breakdown and a delay-repair -period. Recently, Choudhury and

Tadj (2011)studied the optimal management of an Mx_{/G/1 unreli-}

able server queue with optional service under a Bernoulli vacation schedule. Existing work with optional service mentioned earlier, mainly focused on single-server queue . The main reason for this is

such that the steady-state probabili ty vector of a multi-server

queue with SOS is not easily derived. This motivates us to investi- gate an infinite capacity M/M/ R queueing system with a second op- tional service channel which includes parameter optimizati on at a minimum cost.

The paper is organized as follows. In Section 2, the steady-st ate equations and the quasi-birth–death (QBD) model of an infinite capacity M/M/ R queue with SOS channel are set up. The matrix- geometric property (matrix analytic method) is used to calculate the rate matrix in Section 3. In Section 4, we develop an efficient algorithm to obtain the stationary probabilities using a matrix-

geometric and a recursive method. Some system performanc e

measures are derived in Section 5. In Section 6, a cost model is developed to determine the optimal values of the number of chan- nels and the two different service rates used to minimize the total expected cost per unit time. We use the Quasi-Newton method and the direct search method to impleme nt the optimization tasks. The Particle Swarm Optimization (PSO) method is compared with the Quasi-Newton method in establishing the heuristic solution. Some numerical examples are provided to illustrate these two optimiza- tion methods. In Section 7, we offer our conclusions .

2. Markov chain model

For an infinite capacity M/M/ R queueing system with second op- tional service (SOS) channel, the states of the system are described by the pair (i, j), i = 0, 1, 2, . . . and j = 0, 1, 2, . . . , R, where i denotes the

number of customer s in the FES channel (including customers

waiting in the queue) and j denotes the number of customers in the SOS channel. If (i + j) 6 R, i.e., there are available servers, the customers upon arrival to the server will get service immediatel y. If (i + j) > R, i.e., all servers are busy, the newly arriving customer must wait in the queue until a server becomes available. We define the following notations in steady-state:

Pi,j probability that there are i customers in the FES channel

and there are j customers in the SOS channel, where i = 0, 1, 2, . . . and j = 0, 1, 2,. . . , R.

Referring to the steady transition-r ate diagram shown inFig. 1

and using the birth-and-d eath process, the steady-state equations governing the M/M/ R queueing system are:

(i) j = 0

kP0;0¼ ð1  hÞ

l

1P1;0þ

l

2P0;1; ð1Þ

ðk þ i

l

1ÞPi;0¼ kPi1;0þ ði þ 1Þð1  hÞ

l

1Piþ1;0þ

l

2Pi;1; 1 6 i 6 R  1; ð2Þ

ðk þ R

l

1ÞPi;0¼ kPi1;0þ Rð1  hÞ

l

1Piþ1;0þ

l

2Pi;1; R 6 i: ð3Þ

(ii) 1 6 j 6 R  1

ðk þ j

l

2ÞP0; j¼ h

l

1P1; j1þ ð1  hÞ

l

1P1; jþ ðj þ 1Þ

l

2P0; jþ1; ð4Þ

ðk þ i

l

1þ j

l

2ÞPi; j¼ kPi1; jþ ði þ 1Þh

l

1Piþ1; j1þ ði þ 1Þ

ð1  hÞ

l

1Piþ1; jþ ðj þ 1Þ

l

2Pi; jþ1; 1 6 i 6 R  j  1; ð5Þ

½k þ ðR  jÞ

l

₁þ j

l

₂Pi; j¼ kPi1; jþ ðR þ 1  jÞh

l

1Piþ1; j1

þ ðR  jÞð1  hÞ

l

₁Piþ1; jþ ðj þ 1Þ

l

2Pi; jþ1; R  j 6 i: ð6Þ

(iii) j = R

ðk þ R

l

₂ÞP0;R¼ h

l

1P1;R1; ð7Þ

ðk þ R

l

₂ÞPi;R¼ kPi1;Rþ h

l

1Piþ1;R1; 1 6 i: ð8Þ

There is no way of solving Eq.(1)–(8)in a recursive manner to

develop the explicit expressions for the steady-st ate probabilities Pi,j, where i = 0, 1, 2, . . . and j = 0, 1, 2, . . . , R. Alternatively, the

infini-tesimal generator Q of the quasi birth-and-dea th (QBD) process

describin g the M/M/ R queueing system with SOS channel is of

the block-tri diagonal form (seeNeuts, 1981 ):

ð9Þ Each entry of the matrix Q is a square matrix of order R + 1 listed as follows: B ¼ kI; ð10Þ Ai¼ ai;0

l

2 ai;1 2

l

2 ai;2 . . . . . . R

l

₂ ai;R 2 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 5 ; i ¼ 0; . . . ; R ð11Þ Ci¼ ci;0 di;0 ci;1 di;1 ci;2 di;2 . . . . . . ci;R1 di 0 2 6 6 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 7 7 5 ; i ¼ 1; . . . ; R ð12Þ

ai;j¼ ðk þ i

l

1þ j

l

2Þ; 1 6 i þ j 6 R; ½k þ ðR  jÞ

l

1þ j

l

2; i þ j > R;  ð13Þ ci;j¼ ið1  hÞ

l

₁; 1 6 i þ j 6 R; ðR  jÞð1  hÞ

l

1; i þ j > R;  ð14Þ di;j¼ ih

l

₁; 1 6 i þ j 6 R; ðR  jÞh

l

₁; i þ j > R:  ð15Þ Let F ¼ CRþ ARþ B; ð16Þ

It is clear that F is an irreducible generato r (seeNeuts, 1981 ). Let X = [x0, x1, . . . , xR], a 1  (R + 1) vector, is the invariant vector of F.

Then, F satisfies the linear equations:

XF ¼ 0 and Xe ¼ 1; ð17Þ

where 0 and e are column vectors with dimens ions R + 1 and all ele- ments are equal to zero and one, respective ly. Solving the linear Eqs.

(16) and (17), X could be obtained easily. Next, the stability condi-

tion could be establish ed by Theorem 3.1.1 of Neuts (1981), the

standard drift conditio n is:

XBe < XCRe; ð18Þ

which is the necessary and sufficient conditio n for stability of the QBD-Q process. First, we solve xF = 0 x = [x0, x1, . . . , xR] and write fol-

lowing (R + 1) equations:

Rh

l

1x0¼ x1

l

2; ð19aÞ

ðR  i þ 1Þh

l

1xi1þ ½ðR  iÞh

l

1þ i

l

2xi ði þ 1Þ

l

2xiþ1

¼ 0; 1 6 i 6 R  1; ð19bÞ

h

l

1xR1¼ R

l

2xR: ð19cÞ

Eq.(19a)implies that x1¼chll21x0. Solving Eqs.(19b) and (19c) recur-sively, we get: xiþ1¼ ðR  iÞh

l

₁ ði þ 1Þ

l

2 xi; i ¼ 1; . . . ; R  1: ð19dÞ Finally , we have: xiþ1¼ ðc  iÞh

l

1 ði þ 1Þ

l

2 xi¼ ðc  iÞðc  i  1Þ ði þ 1Þi h

l

1

l

2  2 xi1¼  ¼ c i þ 1   _h

_l

1

l

2  iþ1 x0; i ¼ 1;...;R  1: ð19eÞ

Using the normaliz ation conditio n x0+ x1+    + xR1+ xR= 1, x0can

determine x0as: x0¼ XR i¼0 R i   h

l

1

l

₂  i " #1 ¼ 1 þh

l

1

l

₂  R : ð19fÞ

Substitu ting B and CRinto Eq.(18)and using (19f), we have:

l

1ðR  L2Þ > k; ð20aÞ

which is equivalent to: k

l

1ðR  L2Þ

<1; ð20bÞ

where

L2¼ x1þ 2x2þ    þ RxR¼ XR i¼1 ixi¼ XR i¼1 i R i   h

l

1

l

2  i x0 ¼Rh

l

1

l

₂ x0¼ Rh

l

1

l

₂ 1 þ h

l

1

l

₂  R ð21Þ denotes the expected number of customer s in the SOS channe l. It is noted that if h = 0 or

l

2?1 (i.e., L2= 0), Eq.(20)could be reduced

to the stabilit y conditio n for the ordinar y M/M/ R queuein g system without the SOS channe l.

3. The matrix geometric property

Because the infinitesimal generator given in(9)is a special case of Eq. (5.2.1) ofNeuts (1981), we know that the QBD is periodic and a positive recurrent. Denote by P the stationary probability vector of Q. This implies that the unique solution of the system PQ = 0 under stability condition s. We partition the vector P as P = [P0, P1,

P2, . . . , PR1, PR, PR+1, . . .], where Pi= [Pi,0, Pi,1, . . . , Pi,R] is a row vector

with a dimension of R + 1. Our aim is to obtain the steady-state vector P by means of the matrix analytic method and normaliza- tion. By applying the matrix-geom etric method, the steady-state probabilitie s [PR+1, PR+2, PR+3, . . .] can be obtained as Pi= PRTiR, for

i P R + 1, where T is the minimal nonnegative solution to the ma- trix quadratic equation:

T2_C

Rþ TARþ B ¼ 0: ð22Þ

The matrix T is a very important matrix needed in the evalua tion of the performan ce measures of a QBD process. It is known as the rate matrix of the Markov chain Q. Develop ing a closed- form solution for the rate matrix by taking the nonlinea r Eq.(22)is very difficult be- cause the matrix structu re of AR, B and CRis not consiste nt. We de-

velop some matrix analytic properti es to approximat e the rate

matrix T.

Let us decompo se the level space into two groups as

‘(J) = {‘(0), ‘(1), . . . , ‘(R)} and ‘(K) = {‘(R + 1), ‘(R + 2), . . .}. The QBD model of this paper is partially level-dep endent up to a certain le- vel (group ‘(J)) and thereafter becomes a infinite level-independ ent (group ‘(K)). An infinite level-indepe ndent QBD has a matrix-g eo- metric form which can be solved from the matrix quadratic equa-

tion (Latouche & Ramaswa mi, 1999 ). The level-independ ent

structure in our paper can be solved by Cramer’s rule. Thus, we

can use the finite level-depend ent algorithm first and then the

algorithm of infinite level-indepe ndent QBDs to derive the state probabilitie s.

Note from the matrix (9)that starting from level ‘(R) the matri- ces CR1and AR1change to CRand AR, respectively . This implies

that the process holds an infinite level-independ ent QBD with

group ‘(K). First, we reduce the QBD- Q into a finite level-depen- dent QBD- Q⁄_as:

ð23Þ

FromNeuts (1981), we know that the QBD- Q⁄_{is a periodic and irre-}

ducible infinitesimal generato r with finite dimens ions. The matrix H in(23)represe nts the transiti ons between the states belongi ng to the imagina ry level group ‘(K). The boundary steady-stat e prob- ability vector PR+1based on ‘(R + 1) is given by solving the following

equations:

PRB þ PRþ1H ¼ PRþ1; ðfrom QBD- QÞ ð24aÞ

PRB þ PRþ1ARþ PRþ2CR¼ PRþ1: ðfrom QBD- Q Þ ð24bÞ

Solving Eqs.(24a) and (24b), we obtain:

H ¼ ARþ TCR: ð25Þ

Substitutin g(25)into Eq.(23)yields the following system of linear equations : KQ¼ ½P0;P1;P2;...;PRþ1 A0 B C1 A1 B C2 A2 . . . . . . . . . AR B CR ARþ TCR 2 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 5 ¼ 0; ð26Þ

where Pi= [Pi,0, Pi,1, Pi,2, . . . , Pi,R], i = 0, 1, 2, . . . , R + 1.

By the arguments ofLatouche and Ramaswa mi (1999), there ex-

ists an infinitesimal generator U of the transient continuous-ti me Markov chain that is restricted to level ‘(R + 2) before it reaches ‘(R + 1) from group level ‘(J). It is given by:

U ¼ ARþ BðUÞ1CR¼ ARþ BG ¼ ARþ TCR;

where

T ¼ BðUÞ1;

G ¼ ðUÞ1CR;

H ¼ U:

Based on the analysis above, we summarize an algorithm to ob-

tain the approximat ion for the rate matrix T (see Latouche &

Ramaswa mi, 1999 ).

Algorithm 1. Linear Progressi on Algorithm

INPUT B, AR, CR, e = [1, . . . , 1]T, I is the identity matrix, and the

tolerance d.

OUTPUT approximat e solution T Step 1 G = (AR)1CR

Step 2 while ke  Gek P d do Step 3–4 Step 3 set H = AR+ BG Step 4 and G = (H)1_C R Step 5 set T = B(H)1 Step 6 OUTPUT 4. Probability computation

By solving Eq.(24)recursivel y, we obtain:

P0¼ P1C1ðA0Þ1¼ P1/1; ð27aÞ P1¼ P2C2½ð/1B þ A1Þ1¼ P2/2; ð27bÞ P2¼ P3C3½ð/2B þ A2Þ1¼ P3/3; ð27cÞ .. . PR1¼ PRCR½ð/R1B þ AR1Þ1¼ PR/R; ð27dÞ PR¼ PRþ1CR½ð/RB þ ARÞ1¼ PRþ1/Rþ1; ð27eÞ PRT½/Rþ1B þ H ¼ 0: ð27fÞ

where /1= C1(A0)1, /2= C2[(/1B + A1)]1, . . . , /i= Ci[(/i1B +

Ai1)]1, and /R+1= CR[(/RB + AR)]1. Conseque ntly, the Pi

(0 6 i 6 R  1) steady-st ate probabilitie s Pi(0 6 i 6 R  1) in Eqs.

P1¼ PRQ2i¼R/i; . . . ;PR1¼ PRQRi¼R/i, and the rest of the steady-stat e

vector [PR, PR+1, PR+2, . . .] can be determine d recursively using

Pi= PRTiR, for i P R. Once the level probability PRis obtained , the

steady- state solution s [P0, P1, P2, . . . , PR1, PR, PR+1, . . .] can be deter-

mined. The steady-stat e probability PRcan be solved by using the fol-

lowing normaliz ation equation:

X1 n¼0 Pne ¼ PR Y1 i¼R /iþ PR Y2 i¼R /iþ ... þ PR YR i¼R /iþ PRþ PRT þ PRT2þ ... " # e ¼ PR XR k¼1 Yk i¼R /iþ I þ TðI  TÞ 1 " # e ¼ 1: ð28Þ It is clear that we need O(R + 2) equations to obtain the steady-stat e probability PR. Solving Eqs.(27f) and (28)in accordanc e with Cra-

mer’s rule, we obtain PR. Next, computin g the prior state probabil-

ities [P0, P1, P2, . . . , PR1] from (27)we obtain [PR+1, PR+2, . . .] by the

formula, Pi= PRTiR, i P R + 1. We summarize the procedu re below:

Algorithm 2. Recursive Solver

INPUT R, B, A0, A1, . . . , AR, C1, C2, . . . , CR, T, and H

OUTPUT approximat e solution P0, P1, P2, . . .

Step 1 set /1= C1(A0)1

Step 2 for i = 2 to R

Step 3 set /i= Ci[(/i1B + Ai1)]1

Step 4 end

Step 5 set /R+1= CR[(/RB + AR)]1

Step 6 for k = 0 to R  1 Step 7 set Uk¼Qkþ1i¼R/i

Step 8 end Step 9 Solving PRT½/Rþ1B þ H ¼ 0; PR PRk¼1 Qk i¼R/iþ I þ TðI  TÞ1 h i e ¼ 1 Step 10 for i = 0 to R  1 Step 11 set Pi= PRUi Step 12 end Step 13 for i = R + 1 to . . . Step 14 set Pi+1= PiT

Step 15 end Step 16 OUTPUT

5. System performance measures

The system performance measure s, such as the expected num-

ber of customers in the FES channel (denoted by L1), the expected

number of customers in the SOS channel (denoted by L2), the ex-

pected number of idle servers (denoted by E[I]) and the expected number of busy servers in the system (denoted by E[B]), can be evaluated from the steady-state probabilities Pi= [Pi,0, Pi,1, Pi,2,

-. -. -., Pi,R]. The expressions for L1, L2, E[I] and E[B] are given by:

L1¼ X1 i¼1 iPie ¼ XR1 i¼1 iPiþ RPRþ ðR þ 1ÞPRT þ    " # e ¼ X R1 i¼1 iPiþ RPRðI  TÞ1þ PRTðI  TÞ2 " # e; ð29Þ L2¼ X1 i¼0 PiJ ¼ XR1 i¼1 Piþ PRþ PRT þ    " # J ¼ X R1 i¼1 Piþ PRðI  TÞ1 " # J; ð30Þ Ls¼ L1þ L2 ¼X R1 i¼1

Piðie þ JÞ þ PRðI  TÞ1ðRe þ JÞ þ PRTðI  TÞ2e; ð31Þ

E½I ¼X

R1 i¼0

Pi

v

i; ð32Þ

E½B ¼ R  E½I: ð33Þ

where J and e are column vectors with dimension R + 1 as [0, 1, 2, . . ., R]T _{and [1, . . . , 1]}T_{, respective ly. For each 0 6 i 6 R  1,}

the jth element of vector viis max (0, R  i  j + 1), j = 1, 2, . . . , R. That

is,

v

i¼ maxð0; R  i  1 þ 1Þ maxð0; R  i  2 þ 1Þ maxð0; R  i  3 þ 1Þ .. . maxð0; R  i  R þ 1Þ 2 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 5 ¼ maxð0; R  iÞ maxð0; R  i  1Þ maxð0; R  i  2Þ .. . maxð0; i þ 1Þ 2 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 5 ð34Þ

For an infinite capacity M/M/ R queuein g system with the SOS chan- nel, the numerica l results of Lsare obtained by considering the fol-

lowing three cases with different values of R:

Case 1:

l

1= 15,

l

2= 5, h = 0.05, vary the values of k from

0.5 to 10.

Case 2: k = 10,

l

1= 15, h = 0.05, vary the values of

l

2 from

2.5 to 10.

Case 3: k = 10,

l

2= 5, h = 0.05, vary the values of

l

1 from

15 to 25.

Results for LS are depicted inFigs. 2–4for Cases 1–3, respec-

tively. One sees from Fig. 2 that Lsdrastically increases as k

in-creases for R = 1, while Ls slightly increases as k increases for

R P 2. From Figs. 3 and 4we can see that Lsdrastically decrease s

as

l

1or

l

2 increases for R = 1, while LS is not sensitive to

l

1or

l

2for R P 2.

6. Optimization analysis

We construct the total expected cost function per customer per unit time based on the system performanc e measures presented in the previous section. Our main objective is to determine the opti- mum number of servers R, say R⁄_{, and the optimal value of the ser-}

vice rate

l

= (

l

1,

l

2), say

l

¼

l

1;

l

2

 

, simultaneously , so that the expected cost function is minimized. To do this, we define the fol- lowing cost elements:

Ch cost per unit time per customer present in the system,

C1 cost per unit time when one server is busy,

C2 cost per unit time of providing a service rate

l

1,

C3 cost per unit time of providing a service rate

l

2,

C4 fixed cost for purchase of one server.

Using these cost elements listed above, the expected cost func- tion F(R,

l

1,

l

2) per customer per unit time is given by:

FðR;

l

₁;

l

₂Þ ¼ ChLsþ C1E½B þ C2

l

1þ C3

l

2þ C4R: ð35Þ

The cost function in(30)is assumed to be linear in the mean num-

ber of indicate d quantity and it would have been a difficult task to develop analytic results for the optimum value R_;

_l



1;

l

2

 

because

the expected cost function is highly complex and non-linea r in

terms of (R,

l

1,

l

2). In the next sectio n, we first use the QuasiNew

Fig. 2. The expected number of customers in the system versus k.

Fig. 3. The expected number of customers in the system versus l2.

(

l

1,

l

2), say

l

1;

l

2

 

, and then use the direct search method to search for the optimal value of discrete variable R, say R⁄

. 6.1. Quasi-Newton method

In practice, the number of servers is bounded by a positive inte-

ger RUP1 because of the purchase budget. We want to find the

joint optimal value

l



1;

l

2

 

for each given R in the feasible set {1, 2, . . . , RU}. The cost minimizatio n problem can be illustrated

mathematical ly as: F R;

l

 1;

l

2   ¼ min and s:t:ð20ÞfFðR;

l

1;

l

2ÞjRg; R ¼ 1; 2; . . . ; RU ð36Þ

For the optimization problem in(31), it is difficult to show the convexity of F(R,

l

1,

l

2) in (

l

1,

l

2).We note that the derivative of

the cost function F with respect to (

l

1,

l

2) indicates the direction

in which the cost function increases. It means that the optimal va- lue

l



1;

l

2

 

can be found along this opposite direction of the gradi- ent (see Chong & Zak (2001)). That is, for a fixed R, the Quasi-

Newton method is employed to search (

l

1,

l

2) until the minimum

value of F(R,

l

1,

l

2) is achieved, say F R;

l

1;

l

2

 

. An effective proce- dure that makes it possible to calculate the optimal value

R;

l

 1;

l

2

 

is presented as follows: Algorithm 3. Quasi-Newton Method

INPUT Cost function F(R,

l

1,

l

2), R, k, h initial value

l

ð0Þ_¼

_l

ð0Þ 1 ;

l

ð0Þ 2

h iT

, and the tolerance

e

.

OUTPUT approximat ion solution

l



1;

l

2

 T

.

Step 1 Set the initial trial solution for

l

(0)_{, and compute F(}

_l

(0)_).

Step 2 While j@F/@

l

1j >

e

or j@F/@

l

2j >

e

do Step 3–4

Step 3 Compute the cost gradient ~rFð

l

Þ ¼ ½@F=@

l

1; @F=@

l

2 T

and the cost Hessian matrix

Hð

l

Þ ¼ @ 2_F=@

_l

2 1 @ 2_F=@

_l

1@

l

2 @2F=@

l

2@

l

1 @ 2_F=@

_l

2 2 " # at point ~

l

ðiÞ_.

Step 4 Find the new trial solution

l

ðiþ1Þ_¼

_l

ðiÞ_{ ½Hð}

_l

ðiÞ_Þ1_r~

Fð

l

ðiÞ_Þ.

Step 5 OUTPUT

To demonst rate the validity and the process of the optimization

method, some examples are given inTable 1which consider the

following cost parameters

Ch¼ $250=customer=unit-time; C1¼ $180=server=unit-time;

C2¼ $15=unit-time; C3¼ $30=unit-time; and C4¼ $60=server:

Under other given parameters, one can find from Table 1that

the minimum expected cost per unit time of 1682.21 is achieved

at

l



1;

l

2

 

¼ ð27:3756; 14:0267Þ using six iterations, which is R = 3 based on Case (i) with the initial value (

l

1,

l

2) = (20, 10).

Based on Case (ii) with R = 2 and initial value (

l

1,

l

2) = (20, 20),

the minimum expected cost per unit time of 1737.30 is achieved at

l



1;

l

2

 

¼ ð28:8310; 18:7206Þ using six iteration s. 6.1.1. Direct search method

After we obtain the joint optimal value

l



1;

l

2

 

of the contin- uous variable (

l

1,

l

2), we will use the direct search method to

obtain the optimal R such that the expected cost function F R;

l



1;

l

2

 

attains a minimum, say F R_;

_l



1;

l

2

 

. Therefore,

the cost minimizatio n problem can be illustrate d mathemati cally as: F R ;

l

1;

l

2   ¼ min R2f1;2;...;RUg F R;

l

 1;

l

2   ð37Þ The procedure to find the optimal solution is described below. A

numerica l example is shown in Table 2 and is based on: (i)

(k, h) = (15, 0.5) and (ii) (k, h) = (20, 0.8). Algorithm 4. Direct Search Method

INPUT RU, F⁄= M which M is a sufficiently large number

OUTPUT approximat ion solution S

¼ R;

l

 1;

l

2   and F ¼ F R;

l

 1;

l

2   . Step 1 for R = 1 to RU

Step 2 Set a initial trial solution (

l

1,

l

2)

Step 3 Use Quasi-Newton method to find the optimal

value

l



1;

l

2

 

and the cost function F R;

l

 1;

l

2

 

Step 4 If the algorithm is diverge, back to step 2 end if Step 5 If F R;

l

 1;

l

2   <F Step 6 F¼ F R;

l

 1;

l

2   and S¼ R;

l

 1;

l

2   Step 7 end if Step 8 end

Step 9 OUTPUT S⁄_{and F}⁄

Based on Table 2, it is noted that the optimal value

R

;

l

 1;

l

2

 

¼ ð3; 22:86016; 11:64466Þ and the correspondi ng mini-

mum cost F⁄_{= 1463.830 for Case (i). For Case (ii),}

R_;

_l

 1;

l

2

 

¼ ð4; 25:40649; 16:13801Þ and F⁄_{= 1891.530 are opti-}

mal. Finally, we perform a sensitivity investiga tion to the optimal value R_;

_l



1;

l

2

 

based on changes in specific values of the system

paramete rs. The numerical results are shown inTable 3for various

values of h and k. We find that (i) R⁄

increases as k or h increases; and (ii)

l



1

l

2

 

increases as k (h) increases. Moreove r, the minimum expected cost increases as k or h increases.

6.2. Particle Swarm Optimization

In this section, the PSO method is implemented to deal with the cost optimization problem. A comparison between the Quasi- Newton method and the PSO method are also performed. The

PSO algorithm, introduce d by Kennedy and Eberhart (1995)

and Kennedy, Eberhart, and Shi (2001), works by having a pop- ulation of particles and includes the idea of exploitation and explorati on searches . Each particle having a position and velocity denotes a candidate solution. Each particle’s movement is influ-enced by its best known local and global positions in the search- space.

Algorithm 5. Particle Swarm Optimization

INPUT R, initial solution X, learning paramete r w and the tolerance d

OUTPUT approximat e solution ^

l

¼ ½^

l

1; ^

l

2 and FðR; ^

l

Þ

Step 1 Initialization, do Step 2 to Step 4 Step 2 Initialize partial best solution PB = X Step 3 Initialize global best solution

GB = arg min x{F(R, x); x 2 PB}

Step 5 Generate two random numbers U1 U(0, 1) and

U2 U(0, 1)

Step 6 Update the particle’s velocity and positions as:

V ¼ wV þ U1ðPB  XÞ þ U2ðeGB  XÞ and X ¼ X þ V:

Step 7 Update the partial best solution PB and GB = arg min x{F(R, x); x 2 PB}

Step 8 Repeat Step 5 to Step 7 until max{F(R, x1)  F(R, x2); x1, x22 PB} < d

Step 9 Output approximate solution ^

l

¼ GB and FðR; ^

l

Þ

Since the PSO approach does not include the computation of

the gradient, it is suitable for the non-differen tiable objective

function. On the basis of system parameters h = 0.5 and

R = 2, 3, . . . ,7, the Quasi-Newton method is impleme nted with

k= 10, 15, 20, and initial solution

l

(0)_{= [k, 10]} T_{. Based on same}

condition s, the PSO algorithm is executed with an initial solution with 20 particles (randomly generated) and learning paramete r w = 0.2. Some numerical results, including the approximat e opti- mal solutions, computati on time (in seconds), iterations needed to reach convergence and the minimum cost obtained by the

two approaches are shown inTable 4. All numerical results are

obtained by the mathematical program MAPLE 9, which are per- formed using a computer with a CPU-Pentium i3-2100, RAM 4.00 GB.

6.3. Comparis on

A comparative analysis for the two methods is shown in

Ta-ble 4, with changes in initial values of the decision variables for

given system paramete rs. It is noted that (i) the approximat e

optimal solutions and the correspondi ng minimum cost estab-

lished by Quasi-Ne wton method and PSO algorithm are very

close and (ii) the computation time and convergence iteration s

of Quasi-Newton method are evidently less than those of PSO method. That is, the Quasi-Newton method is significantly more effective than PSO method because the PSO method has numer- ous computations /calculations in updating the candidate solution matrix.

Table 1

The illustration of the implement process of the Quasi-Newton method.

Iterations 0 1 2 3 4 5 6

Case (i): (k, h) = (20, 0.5) with R = 3 and initial value (l1,l2) = (20, 10) F(R⁄ ,l1,l2) 1862.22 1735.76 1689.68 1682.43 1682.21 1682.22 1682.21 l1 20 22.7766 25.5320 27.0701 27.3668 27.3756 27.3756 l2 10 11.4360 12.9115 13.8155 14.0192 14.0267 14.0267 @F @l1 32.3746 12.3033 3.53768 0.51339 0.01504 0.00001 3  10 9 @F @l2 74.8311 28.6228 8.49048 1.33248 0.04521 0.00005 4.7  10 9 Ls 2.88890 2.22232 1.83577 1.67455 1.64478 1.64379 1.64379 E[B] 2.00000 1.75253 1.55783 1.46265 1.44412 1.44350 1.44350

Case (ii): (k, h) = (15, 0.8) with R = 2 and initial value (l1,l2) = (20, 20) F(R⁄ ,l1,l2) 1829.50 1760.25 1739.61 1737.33 1737.30 1737.30 1737.30 l1 20 23.8016 27.0887 28.6094 28.8273 28.8310 28.8310 l2 20 19.2062 18.8294 18.7303 18.7207 18.7206 18.7206 @F @l1 30.3036 10.9797 2.84444 0.32356 0.00538 1.8  10 6 _0. @F @l2 6.92424 3.14200 1.01269 0.13222 0.00232 9.9  10 7 3  10(10) Ls 2.26602 1.72458 1.73603 1.66634 1.65691 1.65674 1.65674 E[B] 1.35000 1.25501 1.19104 1.16498 1.16134 1.16128 1.16128

The bold value means that the optimum value ofli.

Table 2

The optimal value (l1,l2) and the corresponding minimum expected cost. R Initial value Coverage value l

1;l2   Iteration Cost ⁄ (i) (k, h) = (15, 0.5) R = 1 [30, 25] [44.20521, 24.33688] 7 2022.146 R = 2 [20, 20] [27.50290, 14.50211] 6 1527.743 R = 3 [15, 15] [22.86016, 11.64466] 6 1463.830 R = 4 [15, 10] [21.33382, 10.71376] 6 1492.969 R = 5 [15, 10] [20.88151, 10.44900] 5 1545.927 (ii) (k, h) = (20, 0.8) R = 1 [50, 30] [61.14970, 40.31473] 9 2890.717 R = 2 [40, 30] [35.80379, 23.29807] 8 2056.578 R = 3 [30, 25] [28.23610, 18.09640] 8 1896.310 R = 4 [25, 20] [25.40649, 16.13801] 5 1891.530 R = 5 [20, 15] [24.38956, 15.44162] 5 1933.145 Table 3

The optimal value R_;l 1;l2

 

and minimum expected value F R_;l 1;l2

 

for various values of k and h.

(k, h) (5, 0.2) (10, 0.2) (20, 0.2) (5, 0.8) (10, 0.8) (20, 0.8) R⁄ 2 2 3 2 3 4 l 1;l2   [13.0953, 4.35200] [19.9021, 6.80977] [26.3424, 8.64436] [13.7175, 8.80645] [18.2622, 11.6276] [25.4065, 16.1380] F R_;l 1;l2   _{729.6488 1011.985 1391.119 976.8809 1356.801 1897.530} Ls 0.690286 0.983412 1.346797 0.958229 1.326524 1.864544 E[B] 0.611596 0.796155 1.221962 0.818713 1.235596 1.778650 (k, h) (10, 0.2) (10, 0.5) (10, 0.8) (20, 0.2) (20, 0.5) (20, 0.8) R⁄ 2 3 3 3 3 4 l 1;l2   [19.9021, 6.80977] [17.9854, 9.09991] [18.2622, 11.6276] [26.3424, 8.64436] [27.37559, 14.02674] [25.4065, 16.1380] F R_;l 1;l2   1011.985 1215.012 1356.801 1391.119 1682.213 1897.530 Ls 0.983412 1.173003 1.326524 1.346797 1.643788 1.864544 E[B] 0.796155 1.105460 1.235596 1.221962 1.443501 1.778650

7. Concluding remarks

In this paper, we modeled an infinite capacity M/M/ R queue using some practical situations wherein arrivals may require an additional optional service (second optional channel by the server). The stationary probabilitie s were able to be efficiently computed by using the two algorithms, in matrix and using matrix approach with the aid of computer software. We also presented efficient search approaches to determine the optimal number of channels and the optimal service rates simultaneou sly to incur minimum cost and we evaluated various system performance measures un- der the optimal operating conditions. The comparison between the Quasi-Newton method and the PSO method were also under- taken. The results rendered are useful in the contexts of modeling banking service systems, computer job processing, automatic ma- chine quality control services channels and many other related applications .

References

Al-Jararha, J., & Madan, K. C. (2003). An M/G/1 queue with second optional service with general service time distribution. Information and Management Sciences, 14, 47–56.

Chong, E. K. P., & Zak, S. H. (2001). An introduction to optimization (2nd ed.). John Wiley and Sons, Inc..

Choudhury, G., & Madan, K. C. (2005). A two-stage batch arrival queueing system with a modified Bernoulli schedule vacation under N-policy. Mathematical and Computer Modelling, 42, 71–85.

Choudhury, G., & Paul, M. (2006). A batch arrival queue with a second optional service channel under N-policy. Stochastic Analysis and Applications, 24, 1–21. Choudhury, G., & Tadj, L. (2009). An M/G/1 queue with two phases of service subject

to the server breakdown and delayed repair. Applied Mathematical Modelling, 33(6), 2699–2709.

Choudhury, G., & Tadj, L. (2011). The optimal control of an Mx

/G/1 unreliable server queue with two phases of service and Bernoulli vacation schedule. Mathematical and Computer Modelling, 54(1-2), 673–688.

Ke, J.-C. (2008). An M[x]

/G/1 system with startup server and J additional options for service. Applied Mathematical Modelling, 32, 443–458.

Kennedy, J., & Eberhart, R. C. (1995). Particle swarm optimization. Proceedings IEEE International Conference on Neural Networks, 4, 1942–1948.

Kennedy, J., Eberhart, R. C., & Shi, Y. (2001). Swarm intelligence . CA: Morgan Kaufmann.

Tadj, L., Choudhury, G., & Tadj, C. (2006a). A quorum queueing system with a random setup time under N-policy and with Bernoulli vacation schedule. Quality Technology and Quantitative Management, 3(2), 145–160.

Tadj, L., Choudhury, G., & Tadj, C. (2006b). A bulk quorum queueing system with a random setup time under N-policy and with Bernoulli vacation schedule. Stochastics: An International Journal of Probability and Stochastics Processes, 78(1), 1–11.

Latouche G. & Ramaswami V. (1999). Introduction to matrix analytic methods in stochastic modeling (ASA-SIAM series on statistics and applied probability). Madan, K. C. (2000). An M/G/1 queue with second optional service. Queueing

Systems, 34, 37–46. Table 4

The comparison between the Quasi-Newton method and the PSO method with various values of k and initial solutio ns.

R 2 3 4 5 6 7 (i) k;lð0Þ 1 ;l ð0Þ 2   ¼ ð10; 10; 10Þ Quasi-Newton method Iterations 9 7 7 7 7 7 CPU time 17.110 20.173 21.015 26.844 35.424 47.046 l 1 20.82313 17.98540 17.18035 16.98173 16.94016 16.93263 l 2 10.88468 9.099914 8.608544 8.493474 8.470412 8.466353 F R;l 1;l2   _{1228.797 1215.012 1259.429 1316.463 1375.963 1435.886}

Particle Swarm Optimization method

Iterations 22 23 18 13 27 15 CPU time 119.765 187.532 159.459 158.098 417.968 327.624 l 1 20.00670 17.96513 17.17655 16.98197 16.73007 16.75662 l 2 9.544196 9.161628 8.615406 8.493446 8.270186 8.534586 F R;l 1;l2   _{1238.469 1215.027 1259.429 1316.464 1376.148 1435.930} (ii) k;lð0Þ 1 ;l ð0Þ 2   ¼ ð15; 15; 10Þ Quasi-Newton method Iterations 10 7 7 6 6 6 CPU time 19.798 19.344 21.265 21.000 28.017 34.047 l 1 27.50290 22.86016 21.33382 20.88151 20.76724 20.74225 l 2 14.50211 11.64466 10.71376 10.44900 10.38488 10.37130 FðR;l 1;l2Þ 1527.743 1463.830 1492.969 1545.927 1604.499 1664.242

Particle Swarm Optimization method

Iterations 25 14 16 21 17 29 CPU time 150.061 115.296 133.749 248.909 256.893 559.404 l 1 23.42697 22.86019 21.33448 20.95719 20.76728 20.42189 l 2 14.68239 11.64474 10.71396 10.31709 10.38477 10.90717 F R;l 1;l2   _{1542.834 1463.831 1492.968 1545.984 1604.499 1665.103} (iii) k;lð0Þ 1 ;l ð0Þ 2   ¼ ð20; 20; 10Þ Quasi-Newton method Iterations 15 7 6 6 6 6 CPU time 27.359 17.998 16.859 21.357 27.64 38.028 l 1 33.86672 27.37559 25.03292 24.24485 24.01695 23.95997 l 2 17.95195 14.02674 12.60472 12.14069 12.01170 11.98050 F R;l 1;l2   1799.006 1682.213 1693.087 1740.360 1797.412 1856.803 Particle Swarm Optimization method

Iterations 15 24 20 18 15 23 CPU time 83.075 189.732 168.342 206.8 233.689 455.797 l 1 33.86676 27.42405 25.05725 24.24761 24.01718 23.95270 l 2 17.95193 13.96756 12.60232 12.12005 12.01238 11.95651 F R;l 1;l2   _{1799.007 1682.223 1693.087 1740.360 1797.414 1856.803}

Medhi, J. (2002). A single server Poisson input queue with a second optional channel. Queueing Systems, 42, 239–242.

Neuts, M. F. (1981). Matrix geometric solutions in stochastic models: An algorithmic approach. Baltimore: The John Hopkins University Press. 1981.

Wang, J. (2004). An M/G/1 queue with second optional service and server breakdowns. Computers and Mathematics with Applications, 47, 1713–1723.

Yang, D.-Y., Wang, K.-H., & Kuo, Y.-T. (2011). Economic application in a finite capacity multi-channel queue with second optional channel. Applied Mathematics and Computation, 217 (18), 7412–7419.