Contents lists available atScienceDirect
Journal of Computational and Applied
Mathematics
journal homepage:www.elsevier.com/locate/cam
Computational algorithm and parameter optimization for a multi-server
system with unreliable servers and impatient customers
Chia-Huang Wu
a, Jau-Chuan Ke
b,∗aDepartment of Industrial Engineering and Management, National Chiao Tung University, Taiwan, ROC
bDepartment of Applied Statistics, National Taichung Institute of Technology, No. 129, Sec. 3, Sanmin Road, Taichung, 404, Taiwan, ROC
a r t i c l e i n f o
Article history:
Received 20 October 2008
Received in revised form 27 April 2010 Keywords:
Cost Balk
Quasi-Newton method Renege
Matrix analytic approach
a b s t r a c t
We consider an infinite capacity M/M/c queueing system with c unreliable servers, in which the customers may balk (do not enter) and renege (leave the queue after entering). The system is analyzed as a quasi-birth-and-death (QBD) process and the necessary and sufficient condition of system equilibrium is obtained. System performance measures are explicitly derived in terms of computable forms. The useful formulae for computing the rate matrix and stationary probabilities are derived by means of a matrix analytical approach. A cost model is derived to determine the optimal values of the number of servers, service rate and repair rate simultaneously at the minimal total expected cost per unit time. The parameter optimization is illustrated numerically by the Quasi-Newton method.
© 2010 Elsevier B.V. All rights reserved.
1. Introduction
In this paper, the quasi-birth-and-death (QBD) process and the matrix analytic method are used to analyze an infinite capacity multi-server M/M/c queue with unreliable servers, balking and reneging customers. The computational algorithm of stationary probability vectors and optimization of parameters are developed.
In most studies on queueing systems, the customers always wait in the system until his service is completed. In many practical systems, such as telephone switchboard customers, hospital emergency rooms’ handling of critical patients, and perishable goods storage inventory systems, the customers may become impatient and leave (i.e., balk or renege) the system without getting services when the waiting time is intolerable. For example, for a call-in customer who cannot get service immediately by the server, he/she is told how long he/she needs to wait. The customer might hang up (balk) or hold on (non-balk and waiting). This is a balking behavior of the customer when the queue length or waiting time is too long. In addition, a waiting customer might hang up (renege) if he/she becomes impatient.
Someone who wants to buy a train ticket (or meal ticket) might decide not entering the system (balk) if the waiting line is too long. On the other hand, as a customer waiting in the queue, he/she might leave the queue (renege) and choose an automat (or instant food). Interesting examples of the occurrence of balking and reneging in air defence systems can be found in [1,2]. In addition, balking and reneging are also common phenomena in telecommunication networks and machine repair problems (see [3,4]).
Queueing models with balking, or reneging, or both have attracted much attention from numerous researchers since Haight [5,6]. The extensions of their basic model can be found in [7–11,3,12]. Later, Wang and Chang [13] examined
∗Corresponding author. Tel.: +886 4 22196077; fax: +886 4 22196331. E-mail address:[email protected](J.-C. Ke).
0377-0427/$ – see front matter©2010 Elsevier B.V. All rights reserved.
a finite capacity M/M/R queueing system with balking, reneging, and server breakdowns. Al-Seedy [14] proposed a transient solution of the non-truncated queue M/M/2 involving balking, and an additional server for longer queues. The maximum likelihood estimates and confidence intervals of an M/M/R/N queue with balking and heterogeneous servers were investigated in [15]. Recently, Ke [16] gave the operating characteristics analysis of the M[x]
/
G/
1 system with a variant vacation policy and balking using supplementary variable technique. Yue and Yue [17] analyzed a finite capacity M/M/c/N queueing system with balking, reneging and synchronous vacations. In [17] they derived the steady-state probability vector expressed as the inverse of two matrices with the blocked matrix forms.Alternatively, many queueing systems were studied with assumption that the server would never fail. In practice, we often encounter cases where the server may fail and can be repaired. Recently, several researches devoted efforts to investigate the impact of unreliable servers (server breakdowns), in which the server subject to unpredictable breakdowns and can be repaired. Past work of unreliable-server queueing models may be divided into two categories, according to whether the system was studied from finite capacity or infinite capacity. Plenty studies in the first category focused on machine repair problems, interested readers may refer to Ke and Wang [3], references there in. Recently, Wang et al. [4] studied the models in [3] and then optimized the number of servers, balking and reneging rates using direct search and steep decent methods. The second category of authors dealing with papers treating the problem of infinite capacity. A pioneering work in this queueing situation is investigated in [18], who first introduced the concept of service interruptions (server breakdowns). Gaver [18] studied an M/G/1 queue with service interruptions by using the embedded Markovian chain. Tang [19] studied some queueing problems of the system and some reliability problems of the broken server for a single unreliable-server M/G/1 queueing system. The reliability measures, were examined in [20], for the ordinary M/G/1 queue with channel breakdowns and vacations. Wang [21] investigated the reliability behavior of the ordinary M/G/1 queue with server breakdowns and second optional service. Wang et al. [22] used maximum entropy principle to approximate the steady-state probability distributions of M[x]
/
M/
1 with server breakdowns and vacations. Recently, Choudhury and Tadj [23] generalized this type of model by introducing the concept of a server breakdown and a delay-repair-period. An M/G/1 retrial queueing system with two phases of service subject to the server breakdown and repair was investigated in [24], who derived the queue size distribution at a random epoch and departure epoch using supplementary variable technique, various system performance measures are also presented. As related works with control policy, the readers can refer to the excellent survey in [25]. Recently, Ke [26] derived the system characteristics and examined the optimal NT policies for M/G/1 system with server breakdown and startup using stochastic decomposition property. Various system performance measures and sensitivity investigations based on the optimal threshold N at a minimum cost, were studied in [27], for the M/G/1 queueing system under N policy with server startup and breakdowns.Existing unreliable-server queues, including those above, mainly focused on finite capacity multi-server system or infinite
capacity single-server system. Because more complicated structure of the stochastic processes required describes system
behaviors, the infinite capacity multi-server queue with unreliable servers is known to be analytically intractable. Analytic steady-state solutions of an infinite capacity M/M/c queue with unreliable servers have not been found. This motivates us to investigate an infinite capacity multi-server queue M/M/c type with unreliable servers, balking and reneging customers. In this paper, we consider an infinite capacity M/M/c queueing system with unreliable servers and impatient customers. Customers arrive according to a Poisson process with parameter
λ
and their service times are provided by c unreliable servers, in which the service times are assumed to be exponentially distributed with mean 1/µ
. A customer on arrival findsn customers and i breakdown (unavailable) servers in the system, either decides to enter the queue with probability bi,nor
balks (do not enter) with probability 1
−
bi,n. If the service is unoccupied or is not interrupted by a breakdown, an arrivingcustomer immediately starts getting the service. But if a customer enters the queue, it may get impatient and leave the queue without getting service. After entering the queue, each customer will wait a certain time t for service to begin, may leave it without being served. This time t is a random variable which is assumed to be distributed according to the exponential distribution with parameter r. We assume that customers arrive at the server form a single waiting line and are served in the order of their arrivals; that is, the first-come, first served discipline. Whenever the server is working, it is assumed that the server can break down at any time with a Poisson breakdown rate
α
. Whenever the server fails, it is immediately repaired at a repair rateβ
, where the repair times are assumed to be exponentially distributed. The server can break down even if no customers are in the system. Each server can serve only one customer at a time, and that the service is independent of the arrival of the customers. Service is additive and allowed to be interrupted if the server breaks down, and the server is immediately repaired. Once the broken down server is repaired, it immediately returns to serve the customer. Although no service occurs during the repair period of all broken down servers, customers continue to arrive according to a Poisson process.The paper is organized as follows; In Section2, mathematical model and the quasi-birth-and-death (QBD) model of an infinite capacity M/M/c queue with unreliable servers and impatient customers are set up. The matrix-geometric property (matrix analytic method) is used to calculate rate matrix in Section3. Section4we derive an efficient algorithm to stationary probabilities by matrix-geometric method. Some system performance measures are derived in Section5. In Section6, a cost model is developed to determine the optimal values of number of channels, service rate and repair rate, simultaneously, in order to minimize the total expected cost per unit time. We use direct search method and Quasi-Newton method to implement the optimization tasks. Some numerical examples are provided to illustrate the two optimization methods. Section 7 concludes.
2. Mathematical model
We consider an infinite capacity M/M/c queueing system with balking, reneging and server breakdowns. The state of the system can be described by the pair
(
i,
n),
i=
0,
1,
2, . . . ,
c,
n=
0,
1,
2, . . .
, where i denotes the number of broken down servers, and n denotes the number of customers in the system.According to system assumptions, a customer who on arrival finds n customers in the system, either decides to enter the queue with probability bi,nor balk with probability 1
−
bi,n, where bi,nare defined as followsbi,n
=
1
,
0≤
n≤
c−
i−
1ρ(
1−
ρ)
n−(c−i)+1,
0≤
c−
i≤
n.
If n is less than or equal to
(
c−
i)
, the phenomenon of reneging will not occur and the customers which upon the server will get services immediately. Otherwise, there are[
n−
(
c−
i)]
customers have to wait in the queue. Furthermore, since the waiting time is assumed to be distributed according to an exponential distribution with mean 1/γ
, the average reneging rate in this state is given byγ
n= [
n−
(
c−
i)]γ
.In steady-state, we define
Pi
(
n) ≡
probability that there are n customers in the system when there are iservers broken down, where i
=
0,
1,
2, . . . ,
c and n=
0,
1,
2, . . .
andΠ
= [
Π0,
Π1, . . . ,
Πn,
Πn+1, . . .]
withΠn= [
P0(
n),
P1(
n), . . . ,
Pc(
n)],
n=
0,
1,
2, . . .
denotes the steady-stateprobability vector.
Referring to the state-transition-rate diagram as shown inFig. 1, it is noted that changes in the state of the system may occur due to one of the following events: (1) a non-balking arrival, (2) a customer reneges, (3) a server breakdown, (4) a repair completion, or (5) a service completion. The first two events depend upon the number of customers n in the system (also see the preceding definition of bi,nand
γ
n). As n increases, the non-balking rate bi,ndecreases progressively. Whenn is sufficiently large, most new arriving customers would balk and do not enter the system. That is, the likelihood of a
non-balking (successful) arrival is small and converges to a very small value when n
≥
N, where N is a sufficiently large.Therefore, the balking rate is reasonably modifies as
bi,n
=
1,
0≤
n≤
c−
i−
1ρ(
1−
ρ)
n−(c−i)+1,
c−
i≤
n≤
Nρ(
1−
ρ)
N−(c−i)+1,
n>
N.
Based on aforementioned observations, the balking rate is larger when n
≥
N, the increasing of the number of customers inthe system is limited and consequently the reneging rate is restricted as well. For analysis, the reneging rate
γ
i,nis redefinedas
γ
n=
0
,
n<
c−
imin
{[
n−
(
c−
i)], [
N−
(
c−
i)]}γ ,
n≥
c−
i.
From [28], the infinitesimal generator Q of the QBD describing the M/M/c queueing system with balking, reneging and c non-reliable servers as Q
=
A0 B0 C1 A1 B1 C2 A2 B2... ... ...
CN−1 AN−1 BN−1 CN AN BN CN AN BN... ... ...
.
(1)The sub-matrices Bk
(
0≤
k≤
N),
Ak(
0≤
k≤
N)
, and Ck(
1≤
k≤
N)
are diagonal matrices of order(
c+
1)
which definedby Bk
=
diag(λ, λ, . . . , λ
|
{z
}
#=c−k, λ
bc−k,k, λ
bc−k+1,k, . . . , λ
bc,k|
{z
}
#=k+1),
0≤
k≤
c diag(λ
b0,k, λ
b1,k, . . . , λ
bc,k),
c+
1≤
k≤
N Ck=
diag(
kµ,
kµ, . . . ,
kµ,
|
{z
}
#=c−k kµ, (
k−
1)µ + γ , . . . ,
kγ
|
{z
}
#=k+1),
0≤
k≤
c diag(
cµ + (
k−
c)γ , (
c−
1)µ + (
k−
c+
1)γ , . . . ,
kγ ),
c+
1≤
k≤
NFig. 1. State-transition-rate diagram for an infinite capacity M/M/c queueing system with unreliable servers and impatient customers.
where ‘‘diag
(
V)
’’ denotes a diagonal matrix with diagonal elements equal to V .Ak
=
ak,0 cα
β
ak,0(
c−
1)α
2β
ak,2(
c−
2)α
... ...
...
...
...
...
(
c−
1)β
ak,c−1α
cβ
ak,c
,
k=
0,
1,
2, . . . ,
Nwith diagonal elements
ak,s
=
−[
λ + (
c−
s)α +
sβ +
kµ],
0≤
s≤
c−
k−
1−[
λ
bs,k+
(
c−
s)α +
sβ + (
c−
s)µ + (
s+
k−
c)γ ],
c−
k≤
s≤
cif 0
≤
k≤
c ak,s= −[
λ
bs,k+
(
c−
s)α +
sβ + (
c−
s)µ + (
s+
k−
c)γ ],
0≤
s≤
c if c+
1≤
k≤
N.
The steady-state probability vectorΠbe the unique solution toΠQ
=
0 andP
∞n=0Πne
=
1, where e is a column vectorfollowing properties
ΠN+k
=
ΠNRk,
for k≥
1.
(2)The matrix R is the unique non-negative solution with spectral radius less than one of the equation
BN
+
RAN+
R2CN=
0.
(3)From [28,29], we know that R is given by limn→∞Rn, where the sequence
{
Rn}
is defined byR0
=
0,
and Rn+1= −
BNA−N1−
R2
nCNA−N1
,
for n≥
0.
(4)The sequence
{
Rn}
is monotone so that R could be evaluated from(4)by successive substitutions.Algorithm: Quasi Progression Algorithm
INPUT BN
,
AN,
CN, e is a (c+
1) column vector with all elements equal to one, and toleranceε
OUTPUT approximate solution R Step 1 Set R
=
0Step 2 while
|
eT(
BN+
RAN+
R2CN)
e| ≥
ε
do step 3–4Step 3 Set Rnew
= −
BNA−N1−
R2CNA−N1Step 4 Set R
=
RnewStep 5 OUTPUT
It is also well known (Theorem 3.1.1 of [28]) that the steady-state probability vector exists if and only if
xBNe
<
xCNe,
(5)where x is the invariant probability of the matrix F
=
CN+
AN+
BN. x satisfies xF=
0 and xe=
1. First we solve xF=
0,where x
= [
x0,
x1, . . . ,
xc]
. We can write following(
c+
1)
equations−
cx0α +
x1β =
0,
(6a)(
c−
i+
1)
xi−1α −
xi[
(
c−
i)α +
iβ] + (
i+
1)
xi+1β =
0,
1≤
i≤
c−
1,
(6b)xc−1
α −
cxcβ =
0.
(6c)Eq.(6a)implies that x1
=
cβαx0, and solving(6b)–(6c)recursively, we getxi+1
=
(
c−
i)α
(
i+
1)β
xi,
i=
1, . . . ,
c−
1.
(7) Finally, we have xi+1=
(
c−
i)α
(
i+
1)β
xi=
(
c−
i)(
c−
i−
1)
(
i+
1)
iα
β
2 xi−1= · · · =
c i+
1α
β
i+1 x0,
i=
1, . . . ,
c−
1.
(8)Also using the normalization condition x0
+
x1+ · · · +
xc−1+
xc=
1, we can determine x0asx0
=
"
cX
i=0 c iα
β
i#
−1=
1+
α
β
−c.
(9)Substituting BNand CNinto Eq.(5)and doing some routine manipulations, then we have
λ
X
c i=0 c iα
β
i bi,N<
cX
i=0{
(
c−
i)µ + [
N−
(
c−
i)]γ }
c iα
β
i,
(10) or equivalentλ
Ei[
bi,N]
<
E[
H]
µ + (
N−
E[
H]
)γ ,
(11) where Ei[
bi,N] =
cX
i=0 c iα
β
i bi,N 1+
α
β
−c,
E[
H] =
cX
i=0(
c−
i)
c iα
β
i 1+
α
β
−cmean the average of entering probability (no balking) and the average number of normal (no breakdown) servers, respectively. As
γ =
0 and bi,N=
1 for all i, Eq.(11)can be reduced to stability condition of the ordinary M/M/c queueing3. Algorithm solution
Under the stability condition, by solving the equation PQ
=
0 with the normalization condition, we obtainΠ0A0
+
Π1C1=
0,
(12a) Πn−1Bn−1+
ΠnAn+
Πn+1Cn+1=
0,
1≤
n≤
N−
1,
(12b) ΠN−1BN−1+
ΠNAN+
ΠNRCN=
0,
(12c) ΠNRn−N−1BN+
ΠNRn−NAN+
ΠNRn−N+1CN=
0,
N+
1≤
n,
(12d) ∞X
n=0 Πne=
1.
(13)After doing routine substitutions to(12a)–(12c), we have
Π0
=
Π1C1(−
A0)
−1=
Π 1φ
1,
Πn−1=
ΠnCn[−
(φ
n−1Bn−2+
An−1)]
−1=
Πnφ
n,
2≤
n≤
N,
(14) and ΠNφ
NBN−1+
ΠNAN+
ΠNRCN=
0.
(15)Consequently,Πn
(
0≤
n≤
N−
1)
in Eq.(14)can be written in terms ofΠNasΠn=
ΠNΠin=+N1φ
i,
n=
0,
1,
2, . . . ,
N−
1and the rest steady-state vectorsΠN
,
ΠN+1, . . .
can be calculated recursively asΠn=
ΠNRn−N, for n≥
N. OnceΠNisdetermined, the steady-state solutionsΠ
= [
Π0,
Π1, . . . ,
ΠN,
ΠN+1, . . .]
are obtained. The vectorΠN is given by solvingEq.(15)with the following normalization condition. ∞
X
n=0 Πne= [
Π0+
Π1+ · · · +
ΠN−1+
ΠN+
ΠN+1+
ΠN+2+ · · · ]
e= [
ΠNΠi1=Nφ
i+
ΠNΠi2=Nφ
i+ · · · +
ΠNΠiN=Nφ
i+
ΠN+
ΠNR+
ΠNR2+ · · · ]
e=
ΠN"
NX
n=1 Πn i=Nφ
i+
(
I−
R)
−1#
e=
1.
(16)Solving Eqs. (15) and (16) in accordance with Cramer’s rule, we obtain ΠN. Then the prior state probabilities
[
Π0,
Π1,
Π2, . . . ,
ΠN−1]
are computed from (14) and[
ΠN+1,
ΠN+2,
ΠN+3, . . .]
are gained by the formula Πn=
ΠNRn−N
,
n≥
N+
1. The solution procedure of the steady-state probabilities is summarized as below:Algorithm: Recursive Solver
INPUT c
,
N,
B0,
B1, . . . ,
BN, A0,
A1, . . . ,
AN, C1,
C2, . . . ,
CN,
R.OUTPUT approximate solutionΠ0
,
Π1,
Π2, . . .
Step 1 set
φ
1=
C1(−
A0)
−1 Step 2 for i=
2 to N Step 3 setφ
i=
Ci[−
(φ
i−1Bi−2+
Ai−1)]
−1 Step 4 end Step 5 for k=
1 to N Step 6 setΦk=
Πik=Nφ
i Step 7 endStep 8 SolvingΠN
φ
NB+
ΠNAN+
ΠNRCN=
0, andΠN[
P
Nk=1Φk
+
(
I−
R)
−1]
e=
1Step 9 for i
=
0 to N−
1 Step 10 setΠi=
ΠNΦi+1Step 11 end
Step 12 for i
=
N+
1 to infinity Step 13 setΠi+1=
ΠiRStep 14 end Step 15 OUTPUT
4. System performance measures
In this section, we derive some system performance measures of the system such as the expected number of customers in the system (denoted by Ls), the expected number of customers in the queue (denoted by Lq), the expected number of busy,
are given by Ls
=
∞X
n=1 nΠne=
N−1X
n=1 nΠne+
NΠNe+
(
N+
1)
ΠNRe+ · · ·
=
N−1X
n=1 nΠNΦn+1e+
NΠN(
I−
R)
−1e+
ΠNR(
I−
R)
−2e=
ΠN"
N−1X
n=1 nΦn+1+
N(
I−
R)
−1+
R(
I−
R)
−2#
e (17) Lq=
Π1
0...
0 1
+
Π2
0...
1 2
+ · · · +
ΠN
N−
c...
N−
1 N
+
ΠNR
N−
c+
1...
N N+
1
+ · · ·
=
N−1X
n=1 ΠNΦn+1un+
ΠN(
I−
R)
−1uN+
ΠNR(
I−
R)
−1e=
ΠN"
N−1X
n=1 Φn+1un+
(
I−
R)
−1uN+
R(
I−
R)
−1e#
(18) E[
D] =
∞X
j=0 Πj
0 1...
c
=
ΠN"
NX
j=1 Φj+
(
I−
R)
−1#
0 1...
c
(19) E[
I] =
Π0
c c−
1...
0
+
Π1
c−
1 c−
2...
0
+
Π2
c−
2 c−
3...
0
+ · · · +
Πc−1
1 0...
0
=
ΠNΦ1v
0+
ΠNΦ2v
1+ · · · +
ΠNΦcv
c−1=
ΠN cX
j=1 Φjv
j−1 (20) E[
B] =
c−
E[
D] −
E[
I]
(21)where un is a
(
c+
1)
dimensional column vector with the kth element equal to max{
0,
n−
(
c−
k+
1)}
andv
j=
[
c−
j,
c−
j−
1, . . . ,
1,
0|
{z
}
#=c−j+1,
0, . . . ,
0|
{z
}
#=j]
T.We also give the steady-state availability and using the concept of Ancker and Gafarian [1,2], the expected balking rate, the expected reneging rate and the expected rate of customer loss are obtained as follows:
1. The steady-state availability
AV
=
ΠN"
NX
i=1 Φi+
(
I−
R)
−1#
ec+1 (22)where ec+1is a column vector with dimension
(
c+
1)
and all elements equal to 1 except that the(
c+
1)
th element equalto zero.
2. The expected balking rate
B
.
R. = λ
Π0ρ
0+
λ
Π1ρ
1+ · · · +
λ
ΠNρ
N+
λ
ΠNRρ
N+
λ
ΠNR2ρ
N· · ·
=
λ
ΠN[
Φ1ρ
0+
Φ2ρ
1+ · · · +
ΦNρ
N−1+
(
I−
R)
−1ρ
N]
(23) whereρ
jis defined byρ
j=
[0, . . . ,0,0 | {z } #=c−j ,1−bc−j,j, . . . ,1−bc,j | {z } #=j+1 ]T, j≤c [1−b0,j, . . . ,1−bc,j]T, c+1≤j≤N+1.Fig. 2. The effect of N on the expected number of customers in the system.
Fig. 3. The effect of N on the steady-state availability.
3. The expected reneging rate
R
.
R. =
Π1Γ1+
Π2Γ2+ · · · +
ΠN−1ΓN−1+
ΠNΓN+
ΠNRΓN+
ΠNR2ΓN· · ·
=
ΠN[
Φ2Γ1+
Φ3Γ2+ · · · +
ΦNΓN−1+
(
I−
R)
−1ΓN]
(24)whereΓjis a column vector with dimension
(
c+
1)
and the kth element equal to max{
0,
j−
(
c−
k+
1)}γ
.4. The expected loss rate L
.
R. =
B.
R. +
R.
R.
To understand how system performance measures listed above vary with N, we now perform some numerical investigation to the measures based on changing the value of N. It should be noted that N initiates with 5 since the number of server which we considered is from 1 to 4 (the readers can refer to system assumptions). As
α =
0.
05,β =
0.
2, γ =
0.
5,λ =
5 andµ =
20, the numerical illustrations for N versus Ls,
AV,
L.
R.
are graphically presented inFigs. 2–4, respectively.FromFigs. 2–4, we observe that the reneging or balking does not result in a commensurate improvement in the system performance measures. Intuitively, N rarely affects the measures. For computation convenience, we adopt N
=
30 in following numerical examples.For an infinite capacity M/M/c queueing system with unreliable servers and impatient customers, the numerical results of Lsare obtained by considering the following two cases with different values of c.
Case 1.
µ =
20, α =
0.
05,β =
0.
2,γ =
0.
5, varyλ
from 1.0 to 5.0. Case 2.λ =
2.
0, α =
0.
05,β =
0.
2,γ =
0.
5, varyµ
from 20 to 25. Results of Lsare depicted inFigs. 5and6for Cases 1–2, respectively.Fig. 4. The effect of N on the expected loss rate.
Fig. 5. The effect ofλon the expected number of customers in the system.
One sees fromFig. 5that Lsdrastically increases as
λ
increases. It reveals fromFig. 6that Lsdrastically decreases asµ
increases. To further understand how system performance measures listed above vary with system parameters, we also perform some numerical investigations to the measures based on changing the value of
λ, µ, α
, andβ
. The following four cases are performed by considering the different values of c.Case 3. Availability versus
α
from 0.05 to 0.45 whenλ =
2.
5, µ =
3.
0,γ =
0.
5, andβ =
0.
5. Case 4. Availability versusβ
from 0.1 to 0.5 whenλ =
2.
5, µ =
3.
0,γ =
0.
5, andα =
0.
05.Case 5. The expected loss rate versus
λ
from 2.5 to 4.5 whenµ =
5.
0, γ =
0.
5,α =
0.
05, andβ =
0.
2. Case 6. The expected loss rate versusµ
from 3.0 to 5.0 whenλ =
2.
5, γ =
0.
5,α =
0.
05, andβ =
0.
2.The numerical illustrations of the steady-state availability and the expected loss rate are graphically presented in
Figs. 7–10for Cases 3–6, respectively. We observe fromFigs. 7and8that the steady-state availability drastically decreases for as
α
increases orβ
decreases for c=
1, while it is not very insensitive toα
orβ
for c≥
2. FromFigs. 9and10, it seems that the expected loss rate is more sensitive to the change ofλ
than the change ofµ
.5. Optimization analysis
In this section, we construct the total expected cost function per unit time based on the system performance measures for such a system in which the number of servers
(
c)
is a discrete decision variable, and then the service rate(µ)
and the repairFig. 6. The effect ofµon the expected number of customers in the system.
Fig. 7. The effect ofαon the steady-state availability.
rate
(β)
are continuous decision variables. Our main objective is to find the optimum number of servers c∗, and the optimum values of service rate and repair rate(µ
∗, β
∗)
to minimum the cost function. Let us define the following cost elements:Ch
≡
holding cost per unit time per customer present in the system;
Cs
≡
cost per unit time of providing an service rateµ;
Cd
≡
cost per unit time when one server is broken down;
Cr
≡
cost per unit time of providing a repair rateβ;
Cl
≡
loss cost per unit time when one customer balks or reneges;
Cp
≡
fixed cost for purchasing one server.
Using the definition of the cost parameters listed above, the total expected cost function per unit time is given by:
F
(
c, µ, β) =
ChLs+
Csµ +
CdE[
D] +
Crβ +
ClL.
R. +
Cpc (25)where Ls
,
E[
D]
, and L.
R.
are defined previously.The analytic study of the optimization behavior of the expected cost function would have been an arduous task to undertake since the decision variables appear in an expression which is a highly nonlinear and complex. Based on the preceding formulation, we use a direct search method to compute/find the optimal value of the number of servers, say
Fig. 8. The effect ofβon the steady-state availability.
Fig. 9. The effect ofλon the expected loss rate.
c∗, when
µ
andβ
are fixed. We then fix c∗and use the Quasi-Newton method to search/adjust the optimal value ofµ
andβ
, sayµ
∗andβ
∗.5.1. Direct search method
The optimum value c∗ can be determined by the following inequality with satisfying Eq. (10)(stability condition) satisfied
F
(
c∗−
1|
µ, β) >
F(
c∗|
µ, β) <
F(
c∗+
1|
µ, β).
(26) Some examples are performed to illustrate the existence of solution. We setµ =
20, β =
0.
8 and consider the following cost parameters asCh
=
$150/
customer/
unit time,
Cd=
$30/
server/
unit time,
Cr=
$45/
unit time,
Cs
=
$15/
unit time,
Cl=
$120/
unit time,
Cp=
$60/
server.
Under other parameters are given, we observe fromTable 1that (i) the optimal number of servers c∗and its corresponding minimum cost increase as
λ
orα
increases; and (ii) the optimum number of servers, c∗does not affect at all whenγ
changes from 0.2 to 0.8. This seems too insensitive to changes inγ
.Fig. 10. The effect ofµon the expected loss rate.
Table 1
The cost function associated with a different number of servers and values ofλ, α, γ.
(λ, µ, α, β, γ ) c=1 c=2 c=3 c=4 c=5 c=6 (5, 20, 0.05, 0.8, 0.2) 560.562 519.969 562.455 621.086 682.394 744.097 (10, 20, 0.05, 0.8, 0.2) 803.648 633.410 619.825 662.907 720.719 781.740 (15, 20, 0.05, 0.8, 0.2) 1133.07 821.187 711.592 716.238 762.168 820.114 (20, 20, 0.05, 0.8, 0.2) N/A 1141.36 870.786 797.168 812.906 861.167 (5, 20, 0.05, 0.8, 0.5) 560.921 520.018 562.462 621.086 682.394 744.097 (10, 20, 0.05, 0.8, 0.5) 808.851 633.706 619.871 662.915 720.720 781.740 (15, 20, 0.05, 0.8, 0.5) 1142.56 821.620 711.672 716.257 762.172 820.115 (20, 20, 0.05, 0.8, 0.5) N/A 1141.36 870.792 797.170 812.907 861.167 (5, 20, 0.05, 0.8, 0.8) 561.137 520.064 562.469 621.087 682.394 744.097 (10, 20, 0.05, 0.8, 0.8) 812.638 633.970 619.915 662.922 720.721 781.740 (15, 20, 0.05, 0.8, 0.8) 1149.01 822.028 711.770 716.276 762.176 820.116 (20, 20, 0.05, 0.8, 0.8) N/A 1141.36 870.797 797.172 812.908 861.167 (λ, µ, α, β, γ ) c=1 c=2 c=3 c=4 c=5 c=6 (5, 20, 0.1, 0.8, 0.5) 589.183 535.802 571.222 628.410 690.473 753.557 (10, 20, 0.1, 0.8, 0.5) 847.825 668.532 639.509 674.618 730.219 791.597 (15, 20, 0.1, 0.8, 0.5) 1190.21 879.252 749.803 738.112 775.973 831.496 (20, 20, 0.1, 0.8, 0.5) N/A 1220.18 934.448 837.518 836.582 876.824 (5, 20, 0.3, 0.8, 0.5) 676.654 598.853 611.031 658.310 718.775 784.148 (10, 20, 0.3, 0.8, 0.5) 968.223 797.030 726.312 731.023 771.576 828.080 (15, 20, 0.3, 0.8, 0.5) 1337.30 1079.67 903.843 842.111 845.522 882.825 (20, 20, 0.3, 0.8, 0.5) N/A 1481.37 1171.85 1015.11 957.632 959.272 (5, 20, 0.5, 0.8, 0.5) 737.368 655.188 653.066 690.071 745.381 809.312 (10, 20, 0.5, 0.8, 0.5) 1051.49 904.517 813.597 795.664 819.693 866.317 (15, 20, 0.5, 0.8, 0.5) 1438.88 1238.18 1045.98 954.811 930.227 946.340 (20, 20, 0.5, 0.8, 0.5) N/A 1677.72 1373.52 1191.21 1097.83 1065.77
* Denotes system is unstable (i.e., the stable condition does not hold).
5.2. Quasi-Newton method
After we find c∗, we will use Quasi-Newton method to search/adjust
(µ, β)
until the minimum value of F(
c∗, µ, β)
is achieved, say F(
c∗, µ
∗, β
∗)
. The cost minimization problem can be illustrated mathematically asF
(
c∗, µ
∗, β
∗) =
Minimizeµ,βs.t. Eq.(10)F
(
c∗
, µ, β).
(27) The finding of the joint optimal value
(µ
∗, β
∗)
for a given c∗is difficult to implement. We note that the derivative of the cost function F with respect to(µ, β)
indicates the direction which cost function increases. It means that, the optimal value(µ
∗, β
∗)
can be found along this opposite direction of the gradient. (See [30]). An effective procedure that makes it possible to calculate the optimal value(µ
∗, β
∗)
is presented as follows:Table 2 The illustration of the implement process of the Quasi-Newton method. Case (i): (λ, α, γ ) = ( 15 , 0 . 1 , 0 . 5 ) with initial value (µ, β) = ( 20 , 0 . 8 ) Iterations 0 1 2 3 4 5 F ( c ∗,µ, β) 738.1120 730.1445 727.7439 727.6845 727.6845 727.6845 c ∗ 4 4 4 4 4 4 µ 20 15.55359 16.61114 16.81859 16.82467 16 .82467 β 0.8 0.968384 0.985683 0.983808 0.983693 0 .983693 ∂ F ∂µ 5.294032 − 4.03993 − 0.566575 − 0.01552 − 0.000013 − 6 . 8 × 10 − 12 ∂ F ∂β − 12.3544 − 6.33340 − 0.525819 − 0.00908 − 0.000007 − 1 . 8 × 10 − 11 Ls 0.735673 0.931753 0.878269 0.868294 0.868004 0.868003 E [ D ] 0.444444 0.374397 0.368432 0.369069 0.369108 0.369108 L . R . 0.320231 0.518905 0.428566 0.415152 0.414787 0.414787 Hessian 1 . 3099 3 . 1502 3 . 1502 156 . 56 3 . 7520 4 . 1618 4 . 1618 111 . 69 2 . 7621 3 . 4215 3 . 4215 98 . 124 2 . 6135 3 . 3275 3 . 3275 97 . 351 2 . 6093 3 . 3252 3 . 3252 97 . 347 2 . 6093 3 . 3251 3 . 3251 97 . 347 Case (ii): (λ, α, γ ) = ( 10 , 0 . 3 , 0 . 5 ) with initial value (µ, β) = ( 20 , 0 . 8 ) Iterations 0 1 2 3 4 5 6 F ( c ∗,µ, β) 726.3119 693.9217 673.7195 671.5941 671.5440 671.5439 671.5439 c ∗ 3 3 3 3 3 3 3 µ 20 12.2204 13.99813 14.58737 14.62913 14.62891 14 .62891 β 0.8 1.13777 1.432544 1.600997 1.639305 1.640787 1 .640789 ∂ F ∂µ 6.46736 − 10.0017 − 2.16940 − 0.178865 − 0.00212 − 1 . 7 × 10 − 6 − 1 . 8 × 10 − 12 ∂ F ∂β − 123.276 − 56.5653 − 15.5898 − 2.36610 − 0.08194 − 0.00011 − 1 . 0 × 10 − 10 Ls 0.480083 0.733697 0.666151 0.645710 0.644638 0.644667 0.644667 E [ D ] 0.818182 0.625970 0.519467 0.473436 0.464084 0.463729 0.463729 L . R . 0.946575 1.254851 0.864804 0.747327 0.711000 0.730524 0.730523 Hessian 0 . 9484 2 . 6956 2 . 6956 427 . 06 5 . 0396 3 . 5377 3 . 5377 170 . 56 3 . 0123 2 . 3415 2 . 3415 84 . 356 2 . 5564 1 . 8817 1 . 8817 59 . 713 2 . 5221 1 . 8093 1 . 8093 55 . 579 2 . 5219 1 . 8069 1 . 8069 55 . 430 2 . 5219 1 . 8069 1 . 8069 55 . 429
Table 3 The optimal value ( c ∗,µ ∗,β ∗) and its minimum expected value for various values of λ, α ,and γ ,while c ∗is obtained at initial value (µ, β) = ( 20 , 0 . 8 ) . (λ, α, γ ) (5, 0.1, 0.5) (10, 0.1, 0.5) (15, 0.1, 0.5) (5, 0.5, 0.5) (10, 0.5, 0.5) (15, 0.5, 0.5) c ∗ 2 3 4 3 4 5 (µ ∗,β ∗) [9.313793, 0.818691] [14.02913, 0.951822] [16.82467, 0.983693] [8.896326,1.510552] [13.41282, 1.833732] [16.42465, 1.968097] F ( c ∗,µ ∗,β ∗) 449.8870 608.4378 727.6845 540.4650 715.4502 849.7029 Ls 0.493536 0.679731 0.868003 0.530139 0.716304 0.889476 E [ D ] 0.217701 0.285219 0.369108 0.746064 0.856996 1.012926 L . R . 0.606480 0.538771 0.414787 0.476188 0.488205 0.424664 (λ, α, γ ) (15, 0.1, 0.2) (15, 0.3, 0.2) (15, 0.5, 0.2) (15, 0.1, 0.8) (15, 0.3, 0.8) (15, 0.5, 0.8) c ∗ 4 4 5 4 4 5 (µ ∗,β ∗) [16.82713, 0.983389] [17.72648, 1.720770] [16.42762, 1.967314] [16.82225, 0.983991] [17.72113, 1.722408] [16.42172, 1. 968856] F ( c ∗,µ ∗,β ∗) 727.6578 793.3554 849.6615 727.7107 793.4822 849.7433 Ls 0.867909 0.816563 0.883965 0.868097 0.816661 0.889586 E [ D ] 0.369212 0.593833 1.013248 0.369007 0.593352 1.012615 L . R . 0.414464 0.581035 0.424300 0.415105 0.582143 0.425021 (λ, α, γ ) (10, 0.1, 0.2) (10, 0.1, 0.5) (10, 0.1, 0.8) (10, 0.5, 0.2) (10, 0.5, 0.5) (10, 0.5, 0.8) c ∗ 3 3 3 4 4 4 (µ ∗,β ∗) [14.02822, 0.951169] [14.02913, 0.951822] [14.03003, 0.952454] [13.41353, 1.831961] [13.41282, 1.833732] [13.41214, 1. 835407] F ( c ∗,µ ∗,β ∗) 608.3219 608.4378 608.5506 715.3165 715.4502 715.5793 Ls 0.679970 0.679731 0.679507 0.716492 0.716304 0.716125 E [ D ] 0.285397 0.285219 0.285048 0.857647 0.856996 0.856382 L . R . 0.537821 0.538771 0.539685 0.487267 0.488205 0.489114
Algorithm: Quasi-Newton Method
INPUT Cost function F
(
c∗, µ, β),
R, λ, µ, α
,β
,γ
, initial value−
→
θ
(0)= [
µ
(0), β
(0)]
T, and the toleranceε
. OUTPUT approximation solution[
µ
∗, β
∗]
T.Step 1 Set the initial trial solution for
−
→
θ
(0), and compute F(
c∗, µ
(0), β
(0))
.Step 2 While
|
∂
F/∂µ| > ε
or|
∂
F/∂β| > ε
do Steps 3–4Step 3 Compute the cost gradient
∇
E
F(
−
→
θ ) = [∂
F/∂µ, ∂
F/∂β]
Tand the cost Hessian matrixH
(
−
→
θ ) =
∂
2F/∂µ
2∂
2F/∂µ∂β
∂
2F/∂β∂µ
∂
2F/∂β
2 at point−
→
θ
(i).Step 4 Find the new trial solution
−
→
θ
(i+1)=
−
→
θ
(i)− [
H(
→
−
θ
(i))]
−1∇
E
F(
−
→
θ
(i))
.Step 5 OUTPUT
We present two examples to illustrate the optimization procedure shown inTable 2. FromTable 2, we can see that the minimum expected cost per day of 727.6845 is achieved at
(µ
∗, β
∗) = (
16.
82467,
0.
983693)
by using 5 iterations, which is c∗=
4 based on Case (i) with initial value(µ, β) = (
20,
0.
8)
. Based on Case (ii), c∗is 3 and the minimum expected cost per day of 671.5439 is achieved at(µ
∗, β
∗) = (
14.
62891,
1.
640789)
by using 6 iterations.Finally, we perform a sensitivity investigation to the optimal values
(
c∗, µ, β)
. For various values ofλ, α
, andγ
by considering the initial value(µ, β)
of (20, 0.8), the minimum expected cost F(
c∗, µ, β)
and the system performance measures Ls,
E[
D]
and L.
R.
at the optimum values(
c∗, µ, β)
are shown inTable 3.FromTable 3, we find that (i) c∗increases as
λ
orα
increases and is insensitive to the change ofγ
; (ii)µ
∗(β
∗)
increases asλ(α)
increases. (iii)µ
∗andβ
∗slightly changes whenγ
changes from 0.2 to 0.8. Intuitively, this seems too insensitive to changes inγ
.6. Concluding remarks
An infinite capacity M/M/c queueing system with balking, reneging and server breakdowns is studied. This system is formulated as a QBD process and the necessary and sufficient condition for the stability of the system is discussed. It generalizes the model studied by Wang and Chang [13]. We have not only obtained numerically the steady-state probability and the system performance measures using matrix approach but also presented one efficient method to find the optimal number of servers, the optimal service rate and repair rate, simultaneously, so as to reach the minimum cost. We also have performed a sensitivity analysis between the joint optimal values
(
c∗, µ
∗, β
∗)
and specific values ofλ, α
, andγ
.Acknowledgements
The authors would like to express their sincere thanks to the referees and the editor for their helpful suggestions that greatly improved the quality of this paper.
References
[1] C.J. Ancker Jr., A.V. Gafarian, Some queueing problems with balking and reneging: I, Oper. Res. 11 (1963) 88–100. [2] C.J. Ancker Jr., A.V. Gafarian, Some queueing problems with balking and reneging: II, Oper. Res. 11 (1963) 928–937.
[3] J.-C. Ke, K.-H. Wang, Cost analysis of the M/M/R machine repair problem with balking, reneging, and server breakdowns, J. Oper. Res. Soc. 50 (1999) 275–282.
[4] K.-H. Wang, J.-B. Ke, J.-C. Ke, Profit analysis of the M/M/R machine repair problem with balking, reneging, and standby switching failures, Comput. Oper. Res. 34 (2) (2007) 835–847.
[5] F.A. Haight, Queueing with balking, Biometrika 44 (1957) 360–369. [6] F.A. Haight, Queueing with reneging, Metrika 2 (1959) 186–197.
[7] M.O. Abou-El-Ata, A.I. Shawky, The single-server Markovian overflow queue with balking, reneging and an additional server for longer queues, Microelectron. Reliab. 32 (1992) 1389–1394.
[8] M.O. Abou-El-Ata, A.M.A. Hariri, The M/M/c/N queue with balking and reneging, Comput. Oper. Res. 19 (1992) 713–716.
[9] R.O. Al-Seedy, The truncated queue: M/M/2/m/m+Y with balking, spares, machine interference and an additional server for longer queues (Krishnamoorthi discipline), Microelectron. Reliab. 35 (1995) 1423–1427.
[10] R.O. Al-Seedy, Analytical solution of the state-dependent Erlangian queue: M/Ej/1/N with balking, Microelectron. Reliab. 36 (1996) 203–206. [11] A.I. Shawky, The single-server machine interference model with balking reneging and an additional server for longer queues, Microelectron. Reliab.
37 (1997) 355–357.
[12] K.-H. Wang, J.-C. Ke, Probabilistic analysis of a repairable system with warm standbys plus balking and reneging, Appl. Math. Model. 27 (2003) 327–336. [13] K.-H. Wang, Y.-C. Chang, Cost analysis of a finite M/M/R queueing system with balking, reneging, and server breakdowns, Math. Methods Oper. Res.
56 (2002) 169–180.
[14] R.O. Al-Seedy, A transient solution of the non-truncated queue M/M/2 with balking and an additional server for longer queues (Krishnamoorthi discipline), Appl. Math. Comput. 156 (3) (2004) 763–769.
[15] K.-H. Wang, S.-C. Chen, J.-C. Ke, Maximum likelihood estimates and confidence intervals of an M/M/R/N queue with balking and heterogeneous servers, RAIRO Oper. Res. 38 (2004) 227–242.
[16] J.-C. Ke, Operating characteristic analysis on the M[x]/
G/1 system with a variant vacation policy and balking, Appl. Math. Model. 31 (7) (2007) 1321–1337.
[17] D. Yue, W. Yue, Analysis of M/M/c/N queueing systems with balking, reneging and synchronous vacations, in: The 2nd Asia-Pacific Symposium on Queueing Theory and Network Applications, QTNA2007, pp. 53–62.
[18] D.P. Gaver, A waiting line with interrupted service, including priorities, J. R. Stat. Soc. Ser. B 24 (1962) 73–96.
[19] Y.H. Tang, A single-server M/G/1 queueing system subject to breakdowns—some reliability and queueing problem, Microelectron. Reliab. 37 (2) (1997) 315–321.
[20] W. Li, D. Shi, X. Chao, Reliability analysis of M/G/1 queueing systems with server breakdowns and vacations, J. Appl. Probab. 34 (1997) 546–555. [21] J. Wang, An M/G/1 queue with second optional service and server breakdowns, Comput. Math. Appl. 47 (2004) 1713–1723.
[22] K.-H. Wang, M.-C. Chan, J.-C. Ke, Maximum entropy analysis of the M[x]/
M/1 queueing system with multiple vacations and server breakdowns, Comput. Ind. Eng. 52 (2) (2007) 192–202.
[23] G. Choudhury, L. Tadj, An M/G/1 queue with two phases of service subject to the server breakdown and delayed repair, Appl. Math. Model. 33 (2009) 2699–2709.
[24] G. Choudhury, K. Deka, An M/G/1 retrial queueing system with two phases of service subject to the server breakdown and repair, Perform. Eval. 65 (10) (2008) 714–724.
[25] L. Tadj, G. Choudhury, Optimal design and control of queues, TOP 13 (2005) 359–414.
[26] J.-C. Ke, Optimal NT policies for M/G/1 system with a startup and unreliable server, Comput. Ind. Eng. 50 (3) (2006) 248–262.
[27] K.-H. Wang, T.-Y. Wang, W.L. Pearn, Optimal control of the N policy M/G/1 queueing system with server breakdowns and general startup times, Appl. Math. Model. 31 (10) (2007) 2199–2212.
[28] M.F. Neuts, Matrix Geometric Solutions in Stochastic Models: An Algorithmic Approach, The John Hopkins University Press, Baltimore, 1981. [29] G. Latouche, V. Ramaswami, Introduction to Matrix Analytic Methods in Stochastic Modeling, in: ASA–SIAM Series on Statistics and Applied Probability,
1999.
