A HEURISTIC ALGORITHM FOR THE OPTIMIZATION OF M/M/s QUEUE WITH MULTIPLE WORKING VACATIONS

(1)

MANAGEMENT OPTIMIZATION

Volume 8, Number 1, February 2012 pp. 1–17

Chia-Huang Wu

Department of Industrial Engineering and Management National Chiao Tung University

Hsingchu 30010, Taiwan

Kuo-Hsiung Wang

Department of Business Administration Asia University

Wufeng, Taichung 41354, Taiwan

Jau-Chuan Ke

Department of Applied Statistics

National Taichung University of Science and Technology Taichung 404, Taiwan

Jyh-Bin Ke

Department of Applied Mathematics National Chung-Hsing University

Taichung 402, Taiwan

(Communicated by Wuyi Yue)

Abstract. This paper focuses on an M/M/s queue with multiple working vacations such that the server works with different service rates rather than no service during the vacation period. We show that this is a generalization of an M/M/1 queue with working vacations in the literature. Service times during vacation period, or during service period and vacation times are all exponentially distributed. We obtain the useful formula for the rate matrix R through matrix-geometric method. A cost function is formulated to determine the optimal number of servers subject to the stability conditions. We apply the direct search algorithm and Newton-Quasi algorithm to heuristically find an approximate solution to the constrained optimization problem. Numerical results are provided to illustrate the effectiveness of the computational algorithm.

1. Introduction. We analyze an M/M/s queue with multiple working vacations such that the server works with variable service rates rather than completely termi- nates service during a vacation period. Such a vacation is called a working vacation (WV) [14]. The server starts a working vacation when the system is empty. When a working vacation ends and there are no customers in the system, the server begins another working vacation. If the server returns from a working vacation and find the not-empty system, he switches to another service rate. The time interval between two successive vacations is called a service period [14, p. 107].

2000 Mathematics Subject Classification. Primary: 60K10, 60K25; Secondary: 90B22, 90B25.

Key words and phrases. Newton-Quasi algorithm; optimization; rate matrix; sensitivity analysis; working vacations.

1

(2)

It is assumed that customers arrive according to a Poisson process with rate λ. The service times during a service period follow exponential distribution with mean 1/µB. The service times during a vacation period follow another exponential distribution with mean 1/µv. When there are no waiting customers in the system, the server begins a working vacation of mean duration 1/η, where vacation times are also exponentially distributed. We assume that arriving customers form a single waiting line based on the order of their arrivals; that is, the first-come, first-served discipline is followed. Suppose that one server can serve only one customer at a time. Customers entering into the service facility and finding that the server is busy have to wait in the queue until the server is available. We also assume that the customer is always assigned or shifted to the server with higher service rate if it is available.

Servi and Finn [14] introduced an Internet Protocol (IP) access network example where each gateway router which connected to the optical network could be formulated as an M/M/1/WV queue. The example illustrated that the use of wavelength division multiplexing (WDM) access network can improve performance in reconfiguration of the wavelengths. A reconfigurable WDM optical access network as Figure 1 is considered. Each access routers with multiple ports are connected to the global IP network through a gateway router. Each port can transfer data over a single wavelength. The wavelength reconfiguration problem is to determine which access router ports should be connected to the gateway router ports according to its wavelengths. A strategy is to reconfigure the additional roving wavelengths to the next router when one access router has empty queue. It is assumed that there are s roving wavelengths (s servers) with the capacity of mean service rate u^∗. For keep- ing the connection, each access router serves permanently at a nominal average rate sµv and all additional roving wavelength are initially configured to access router 1 with a total average service rate of sµv+ su^∗ = sµB. At a service completion, one additional roving wavelength would be assigned to the next router if no more requirement (customer) waiting in the queue. This cycle continually repeats itself.

Figure 1. An optical access Network

(3)

For example, there are s = 5 roving wavelengths which is tunable in the access network. Initially, access router 1 and 2 serve with service rates 5µB and 5µv, respectively. When the number of customers in router 1 changes from 5 to 4 (one idle sever), one roving wavelength would be assigned to router 2. This results in access router 1 instantaneously dropping its average service rate from 5µB to 5µ_v+ 4u^∗ = µ_v+ 4µ_B and access router 2 is operating at an average service rate of 5µ_v+ u^∗= 4µ_v+ µ_B. We consider an optical network with s WDM to configure different roving wavelengths (s servers) which can be allocated to several access routers. From the viewpoint of access router i, it follows an M/M/s queue with working vacation policy. Next, the investigation of queueing systems with server vacations or working vacations are reviewed.

During the last two decades, the queueing systems with server vacations or multiple working vacations have been investigated by many researchers. Past work may be divided into two categories: (i) the category of server vacations, and (ii) the category of multiple working vacations. In the first category the readers may refer to the survey paper by Doshi [6]. The GI/M/1 queue with server vacations have been studied by several authors such as Karaesmen and Gupta [8], Chatterjee and Mukherjee [4], Tian et al. [16], etc. They considered that the vacation time follows exponential distribution in [8], [16] and general distribution in [4], respectively.

Karaesmen and Gupta [8] analyzed a finite capacity GI/M/1 queue with server vacations. They developed the queue length distribution at arrival and random epochs for the multiple vacations case. The model was previously investigated by Chatter- jee and Mukherjee [4] for infinite capacity queue considering server vacations. They applied the embedded Markov chain technique to obtain the steady-state probability distributions of system size at pre-arrival and at random epochs, separately.

Fuhrmann [7] studied a single server queue of M/G/1 type with server vacations.

The finite capacity M/G/1 queue with server vacations was investigated by Lee [9].

In the second category Servi and Finn [14] first examined an M/M/1 queue with working vacations, where inter-arrival times, service times during service period, service times during vacation period, and vacation times are all exponentially distributed. The queueing model is denoted by M/M/1/WV. Recently, Liu et al. [12]

used stochastic decomposition method to analyze M/M/1/WV queue which is different from Servi and Finn’s method. Discrete time GI/Geo/1 queue with multiple working vacations was first studied by Li et al. [10]. Later Wu and Takagi [18]

extended Servi and Finn’s M/M/1/WV queue to an M/G/1/WV queue. They assumed that service times during service period, service times during vacation period as well as vacation times are all generally distributed. Baba [1] extended Servi and Finn’s M/M/1/WV queue to a GI/M/1/WV queue and derived the steady-state system length distributions at arrival and arbitrary epochs. Banik et al. [2] studied a finite capacity GI/M/1 queue with multiple working vacations. They developed some important system performance measures such as, the probability of blocking and the expected waiting time in the system, etc. Recently, Lin and Ke [11] studied a Markovian queueing system with a single working vacation. They provided a closed-form method for obtaining the rate matrix. In the study of working vacations, existing literature focuses mainly on a single server queue with multiple vacations. This would motivate us to investigate the M/M/s/WV queue.

In Section 2, the steady-state equations are solved via matrix-geometric method.

The computation for rate matrix R follows in Section 3. This queueing model is a generalization of an M/M/1 queue with working vacations shown in Section 4.

(4)

In Section 5, the explicit formulae for system performance measures are developed.

Cost analysis and sensitivity investigation are presented in Section 6. Finally, Sec- tion 7 summarizes the paper with some concluding remarks.

2. Steady-state results. The states of the system are described by the pair (k, n), k = 0, 1, 2, . . . , s − 1, s and n = k, k + 1, k + 2, . . ., where k represents the number of working servers and n represents the number of customers in the system. In steady-state, we introduce the following notations:

P0(n) ≡ probability that there are n customers in the system when all servers are on a working vacation, where n = 0, 1, 2, . . .;

P1(n) ≡ probability that there are n customers in the system when there are (s − 1) servers on a working vacation, where n = 1, 2, . . .;

Pk(n) ≡ probability that there are n customers in the system when there are (s−k) servers on a working vacation, where n = k, k + 1, k + 2, . . ., k = 2, 3, 4, . . . , s.

The steady-state equations for P_k(n)(k = 0, 1, 2, . . . , s) relating to Figure 2 are given by:

(i) k = 0

λP0(0) = µBP1(1) + µvP0(1), (1) (λ + nµ_v+ sη)P₀(n) = λP₀(n − 1) + (n + 1)µ_vP₀(n + 1), 1 ≤ n ≤ s − 1(2) (λ + sµv+ sη)P0(n) = λP0(n − 1) + sµvP0(n + 1), n ≥ s (3) (ii) 1 ≤ k ≤ s − 1

(λ + kµ_B)P_k(k) = (k + 1)µ_BP_k+1(k + 1) + (kµ_B+ µ_v)P_k(k + 1) + (s − k + 1)ηP_k−1(k), (4) [λ + kµB+ (n − k)µv+ (s − k)η]Pk(n) = [kµB+ (n − k + 1)µv]Pk(n + 1)

+λP_k(n − 1) + (s − k + 1)ηP_k−1(n), k + 1 ≤ n ≤ s − 1 (5)

[λ + kµB+ (s − k)µv+ (s − k)η]Pk(n) = [kµB+ (s − k)µv]Pk(n + 1)

+λP_k(n − 1) + (s − k + 1)ηP_k−1(n), n ≥ s (6) (iii) k = s

(λ + sµB)Ps(s) = sµBPs(s + 1) + ηPs−1(s), (7) (λ + sµB)Ps(n) = λPs(n − 1) + sµBPs(n + 1) + ηPs−1(n), n ≥ s + 1. (8) In the next section, we apply the matrix-geometric method to solve (1) to (8) for Pk(n), where k = 0, 1, 2, . . . , s and n = k, k + 1, k + 2, . . .. A matrix geometric method was introduced by Neuts [13] for analyzing many complex queueing systems.

Using this method, we derive the steady-state probabilities P_k(n)(k = 0, 1, 2, . . . , s).

(5)

Figure 2. State-transition-rate diagram for an M/M/s/WV queue

The infinitesimal generator for the quasi birth and death process has the block- tridiagonal form (see Neuts [13]):

Q =







B₀ C A₁ B₁ C

A₂ B₂ C

. .. . .. . .. As−1 Bs−1 C

As Bs C As Bs C

As Bs C . .. . .. . ..







. (9)

Each entry of the matrix Q is a square matrix of order s + 1 listed as follows:

C = −λI, (10)

B0=





 λ

−µB λ + µB

−2µB λ + 2µB

. .. . ..

−(s − 1)µB λ + (s − 1)µ_B

−sµB λ + sµ_B





 ,

(11)

(6)

B_i=





 b_i,0

b_i,1 b_i,2

. ..

bi,s−1

bi,s







, i = 1, 2, . . . , s (12)

A_i=







a_i,0 −sη

a_i,1 −(s − 1)η

a_i,2 −(s − 2)η . .. . ..

a_i,s−1 −η ai,s







, i = 1, 2, . . . , s (13)

where I is the identity matrix of order s + 1, and ai,j=

−iµv+ jµB, 1 ≤ i + j ≤ s

−[(s − j)µv+ jµB], i + j > s (14) b_i,j =

λ + iµ_v+ jµ_B+ (s − j)η, 1 ≤ i + j ≤ s

λ + (s − j)µv+ jµB+ (s − j)η, i + j > s. (15) Let P denote the corresponding steady-state probability vector of Q. By parti- tioning the vector P as P = [P0, P1, P2, . . . , Ps−1, Ps, Ps+1, . . .], where Pn = [P0(n), P1(n + 1), P2(n + 2), . . . , Ps(n + s)] is a row vector with dimension (s + 1), we can write the steady-state equations PQ = 0 and normalization condition as

P₀B₀+ P₁A₁ = 0 (16)

Pn−1C + PnBn+ Pn+1An+1 = 0, n = 1, 2, . . . , s − 1 (17) P_n−1C + P_nB_s+ P_n+1A_s = 0, n = s, s + 1, s + 2, . . . (18)

s−1

X

n=1

Pne + Ps(I − R)⁻¹e = 1 (19)

where R is the minimal non-negative solution to the matrix quadratic equation (20) and e is a column vector with dimension (s + 1) whose transpose is [1, 1, . . . , 1].

R²A_s+ RB_s+ C = 0. (20)

Using (16)-(17) and (19), the steady-state vectors [P₀, P₁, P₂, . . . , P_s] can be solved by the block Gauss-Seidel iteration method. The rest steady-state vector [Ps+1, Ps+2, Ps+3, . . .] can be determined recursively by Pn = PsR^n−s, for n ≥ s. Once the rate matrix R is determined, the steady-state solutions [P0, P1, P2, . . . , Ps−1, Ps, Ps+1, . . .] are obtained.

3. Computation for rate matrix R. To obtain [P₀, P₁, P₂, . . . , P_s−1, P_s, P_s+1, . . .], it is necessary to solve the minimal solution of (20). Based on the upper triangular structures of matrices A_s, B_s, and C, the rate matrix solution R is also upper triangular. For the multiple-server queue, we assume that there exists a matrix R satisfying the following equation:

R²As+ RBs+ C = 0,

where A_s, B_s, and C are given in (13), (12), (10), respectively.

(7)

Using the Maple computer program and matrix algorithm, we develop the explicit formula for rate matrix R as follows:

R =







t1,1 t1,2 t1,3 · · · t1,s+1

0 t2,2 t2,3 · · · t2,s+1

. .. ... ... . .. ... ...

t_s,s t_s,s+1 0 t_s+1,s+1







, (21)

where

ti,i = λ + θ_i+ (s − i + 1)η −p[λ + θ_i+ (s − i + 1)η]²− 4λθ_i 2θi

, for 1 ≤ i ≤ s + 1

ti,j =

Pj−1

m=i+1ti,mtm,jθj+ (s − j + 2)ηPj−1

m=iti,mtm,j−1

λ + θ_j+ (s − j + 1)η − (t_i,i+ t_j,j)θ_j , for 1 ≤ i ≤ j ≤ s + 1 ti,j = 0, for i > j, and

θ_j = (s − j + 1)µ_v+ (j − 1)µ_B.

Note that t_i,jis the corresponding eigenvalue of the rate matrix R. Once the rate matrix R is determined, the steady-state solutions can be evaluated.

4. Generalization of Servi and Finn’s result. The queueing model considered in this paper is more general than the work of Servi and Finn [14] who first introduced an M/M/1/WV queue. In this section, we will show that our queueing model generalizes the M/M/1/WV queue. For a singer server queue (that is, s = 1), we obtain from (20) that

R²A1+ RB1+ C = 0, (22)

It yields from (10), (12)-(13) that A1=

−µv −η 0 −µB

, B1=

λ + µv+ η 0

0 λ + µB

, C =

−λ 0

0 −λ

.

Since A₁, B₁, C are 2×2 upper triangle matrices. Therefore we set R =

a b 0 d

. From (22), we get

a b 0 d

−µv −η 0 −µB

+

a b 0 d

λ + µv+ η 0

0 λ + µB

+

−λ 0

0 −λ

= 0. (23) Solving (23) yields

a =λ + µ_v+ η −p(λ + µ_v+ η)²− 4λµ_v 2µv

, b = ηa²

µB(1 − a), d = λ µB

.

Thus, a rate matrix R is obtained. Let P = [P₀, P₁] and P₀= [P₀(0), P₁(1)]. Since PQ = 0 and P₁= P₀R, it implies from (16) that

P₀B₀+ P₁A₁= P₀B₀+ P₀RA₁= P₀(B₀+ RA₁) = 0. (24) We can get P₀= [P₀(0), P₁(1)] using (24) and the normalization condition of

P₀(I − R)⁻¹e = 1, (25)

(8)

where I is the identity matrix of order 2, and e is a column vector whose transpose is [1, 1]. We find that

B0=

λ 0

−µB λ + µB

, (I − R)⁻¹= 1 (1 − a)(1 − d)

1 − d −b 0 1 − a

.

Solving (24) and (25), we obtain

(λ − µva)P0(0) − µBP1(1) = 0, (26) (−ηa − bµ_B)P₀(0) + µ_BP₁(1) = 0, (27) (1 + b − d)P0(0) + (1 − a)P1(1) = (1 − a)(1 − d). (28) Solving (26)-(28) again, after the algebraic manipulation, we get

P0(0) = (1 − a)(µ_B− λ) µB− aµv

. (29)

where a = ^λ+µ^v^+η−

√

(λ+µv+η)²−4λµv

2µv =^λz_µ^∗

v . Note that z^∗is the smaller root of the equation

A(z) = λz²− (λ + η + µ_v)z + µ_v = λ(z − z^∗)(z − ˆz) = 0, for |ˆz| ≥ 1 (30) which is (A.4) of [14, p. 49]. Following Servi and Finn’s result, we have λz^∗z = µˆ _v. It implies that a = λz^∗/µ_v= 1/ˆz. Substituting a = 1/ˆz into (29) yields

P0(0) = (µB− λ)(1 − 1/ˆz)

µ_B− µ_v/ˆz = (µB− λ)(ˆz − 1)

µ_Bz − µˆ _v . (31) Thus from Servi and Finn’s equation (1.2) (see [14, p. 42]), we obtain

P r(N = 0) = r(1 − λ/µ_B) + (1 − r)(1 − ˆz⁻¹)

= z − 1ˆ ˆ

z −(λˆz − µB)(ˆz − 1)(λˆz − µv)µB

µ_Bz(µˆ _Bz − µˆ _v)(λˆz − µ_B)

= z − 1ˆ ˆ z

1 − λˆz − µv

µBˆz − µv

= (ˆz − 1)(µB− λ) µBz − µˆ v

which is in agreement with expression (31).

5. System performance measures and some numerical examples. The system performance measures, such as the average number of customers in the system (Ls), the average number of servers during service period (E[N W ]), and the average number of servers during vacation period (E[W V ]), can be obtained from the steady-state probabilities P_n = [P₀(n), P₁(n + 1), P₂(n + 2), . . . , P_s(n + s)] as

(9)

follows:

Ls =

∞

X

n=0

Pn





 n n + 1 n + 2

... n + s







=

∞

X

n=0

Pn(v + ne)

= P0v + P1(v + e) + · · · + Ps−1[v + (s − 1)e] + Ps(v + se) +P_sR[v + (s + 1)e] + · · ·

=

s−1

X

n=0

Pn(v + ne) + Ps(v + se) + PsR[v + (s + 1)e]

+PsR²[v + (s + 2)e] + · · ·

=

s−1

X

n=0

Pn(v + ne) + Ps(I + R + R²+ · · · )(v + se) +Ps(R + 2R²+ 3R³+ · · · )e

=

s−1

X

n=0

Pn(v + ne) + Ps(I − R)⁻¹(v + se) + PsR(I − R)⁻²e, (32)

E[N W ] =

∞

X

n=0

Pn





 0 1 2 ... s







=

∞

X

n=0

Pnv =

s−1

X

n=0

Pnv + Ps(I + R + R²+ · · · )v

=

s−1

X

n=0

Pnv + Ps(I − R)⁻¹v, (33)

E[W V ] = s − E[N W ] (34)

where v and e are column vectors with dimension (s + 1) whose transposes are denoted by [0, 1, 2, . . . , s] and [1, 1, . . . , 1], respectively. Moreover, the expected customers served by normal working servers (E[LsN W ]), the expected customers served by working vacation servers (E[LsW V ]), and the expected idle working vacation servers that do not serve customers currently (E[IDW V ]) are also obtained as follows:

E[L_sN W ] =

∞

X

n=0

P_n





 0 1 2 ... s







=

s−1

X

n=0

P_nv + P_s(I − R)⁻¹v = E[N W ], (35)

(10)

E[LsW V ]

= P1





 1 1 ... 1 1 0





 + P2





 2 2 ... 2 1 0







+ · · · + Ps





 s s − 1

... 2 1 0







+ Ps+1





 s s − 1

... 2 1 0





 + · · ·

=

s−1

X

n=1

Pn







max{min{n, s}, 0}

max{min{n, s − 1}, 0}

...

max{min{n, s − (s − 2)}, 0}

max{min{n, s − (s − 1)}, 0}

max{min{n, s − s}, 0}







+ Ps(I − R)⁻¹





 s s − 1

... 2 1 0





 , (36)

E[IDW V ] = P₀





 s s − 1

... 2 1 0





 + P₂





 s − 1 s − 2

... 1 0 0







+ · · · + P_s−1





 1 0 ... 0 0 0







=

s−1

X

n=0

Pn







max{0, s − n}

max{0, s − n − 1}

...

max{0, s − n − (s − 2)}

max{0, s − n − (s − 1)}

max{0, s − n − s}







. (37)

It could be verify that E[LsN W ] + E[LsW V ] + E[IDW V ] = s. (Therefore, E[LsW V ] + E[IDW V ] = E[W V ]).

Consider an Internet Protocol (IP) access network with s roving wavelengths, some numerical results about the number of customers in the access router (Ls) are presented graphically on the basis of the following three cases:

Case 1. µB= 3.5, µv = 2.0, η = 0.05, vary the values of λ from 1.0 to 3.0.

Case 2. λ = 1.5, µ_B= 3.5, η = 0.05, vary the values of µ_v from 0.5 to 3.0.

Case 3. λ = 1.5, µ_v= 2.0, η = 0.05, vary the values of µ_B from 2.5 to 4.5.

Results of Lsare depicted in Figures 3-5 for cases 1-3, respectively. Figure 3 shows that Lsdrastically increases as λ increases for s = 1, while Lsslightly increases as λ increases for s ≥ 2. Figure 4 illustrates that Lsdrastically decreases as µv increases for s = 1, while Ls slightly decreases as µv increases for s ≥ 2. From Figure 5, we find that L_sslightly decreases as µ_B increases for any s. It is noted that the more roving wavelengths, the more flexible the optical access network is. Also, diving and transferring the additional wavelength partly would be more elastic than shifting all the additional wavelength at a time.

6. Cost analysis and sensitivity investigation. A steady-state expected cost function per unit time is developed. Suppose the average arrival rate is λ, Our

(11)

Figure 3. The mean number of customers in the system versus λ

object is to determine the value of the number of tunable wavelengths s, nominal average service rate µ_v, and fully average service rate µ_B under stability condition ρ = λ/sµ_B< 1 and a natural (rational) constraint µ_B = µ_v+ u^∗≥ µ_v. The discrete variable s is required to be a positive integer, and the continuous variables µv and µB are positive numbers. The following cost parameters are considered:

C_h ≡ holding cost per unit time per customer present in the system;

C_B ≡ cost per unit time when one server is normal working;

C_V ≡ cost per unit time per customer served by one working vacation server;

C_I ≡ cost per unit when one server is idle on working vacation state;

Cs ≡ cost per unit time of providing a specific server rate.

The cost minimization problem can be written as

s,µmin_v,µ_BF (s, µv, µB) (38) subject to

sµB> λ and µB≥ µv, (39)

where

F (s, µv, µB) = ChLs+ CBE[LsN W ] + CVE[LsW V ] + CIE[IDW V ] + Cs(µB+ µv).

(40)

(12)

Figure 4. The mean number of customers in the system versus µ^v

(see White et al. [17], p. 218). It is extremely difficult to prove the convex property and then develop the analytical results for the optimum values (s^∗, µ^∗_v, µ^∗_B) because the cost function is highly complex.

Due to the fact s is a discrete quantity, µ_v and µ_B are continuous quanti- ties, we attempt to obtain an approximate (local) minimum solution (ˆs, ˆµ_v, ˆµ_B) of (s^∗, µ^∗_v, µ^∗_B) by Newton-Quasi method for various values of s. Newton’s method is one of the most powerful and well-known numerical methods for solving the root- finding problems. Consequently, Newton’s method is employed to solve the root of G(x) = ∇F (s, µv, µB) = 0 which is the first order necessary condition of local minimum solution of function G. From Burden and Douglas [3], the sequence

xn+1= xn− ∇G(xn)⁻¹G(xn) = xn− [H(xn)]⁻¹∇F (x) = 0 (41) could help us to obtain an approximate minimum solution where H(xn) denotes the Hessian matrix of function F . Therefore, equation (41) could be regard as the Newton’s method that aims at the root-finding problem of the first derivation equation. Based on Chong and Zak [5], the steps of Newton-Quasi algorithm can be described as follows for any iteration

Step 1. For each value of s, let µ_n= [µv, µB]^T. Step 2. Set the initial trial solution for µ_n with n = 0.

(13)

Figure 5. The mean number of customers in the system versus µ^B

Step 3. Compute F (s, µ_n), the cost gradient ∇F (s, µ_n) = [∂F/∂µ_v, ∂F/∂µ_B]^T|µ_n and the cost Hessian matrix

H(s, µ_n) =

∂²F/∂µ²_v ∂²F/∂µ_v∂µ_B

∂²F/∂µB∂µv ∂²F/∂µ²_B

µ_n

.

Step 4. Find the new trial solution µ_n+1= µ_n− [H(s, µ_n)]⁻¹· ∇F (s, µ_n).

Step 5. Set n = n + 1 and repeat Steps 3-4 if |∂F/∂µ_v| > ε1 or |∂F/∂µ_B| > ε2 or kµ_n+1− µ_nk > ε3, where ε₁, ε₂, and ε₃ are the tolerances; otherwise, go to Step 6.

Step 6. Find the local minimum value F (s, µ_n) = F (s, µ^∗_v, µ^∗_B).

Because the cost function is complex, the gradient ∇F and the Hessian matrix H may be derived (approximated) numerically while their expressions could not be explicitly obtained. The gradient vector could be computed directly as

∇F (s, µv, µ_B) ≈

F (s, µv+ ∆, µB) − F (s, µv, µB)

∆ ,F (s, µv, µB+ ∆) − F (s, µv, µB)

∆

^T

(42)

(14)

with ∆ is small (For example, 10⁻⁶). The elements of Hessian matrix are approximated as

∂²F (s, µv, µB)

∂µ²_v ≈∂F (s, µv+ ε, µB)/∂µv− ∂F (s, µv, µB)/∂µv

ε

≈

F (s,µ_v+ε+∆,µ_B)−F (s,µ_v+ε,µ_B)

∆ −^{F (s,µ}^v^+∆,µ^B_∆^{)−F (s,µ}^v^,µ^B⁾

ε (43)

∂²F (s, µ_v, µ_B)

∂µv∂µB

≈ ∂F (s, µ_v+ ε, µ_B)/∂µ_B− ∂F (s, µ_v, µ_B)/∂µ_B ε

≈

F (s,µv+ε,µB+∆)−F (s,µv+ε,µB)

∆ −^{F (s,µ}^v^,µ^B^{+∆)−F (s,µ}_∆ ^v^,µ^B⁾

ε (44)

with ε is small. For simplification, the equations (43) and (44) become

∂²F (s, µ_v, µ_B)

∂µ²_v ≈ F (s, µ_v+ 2∆, µ_B) − 2F (s, µ_v+ ∆, µ_B) + F (s, µ_v, µ_B)

∆² (45)

∂²F (s, µ_v, µ_B)

∂µv∂µB

≈ F (s, µ_v+ ∆, µ_B+ ∆) − F (s, µ_v+ ∆, µ_B)

∆²

−F (s, µv, µB+ ∆) − F (s, µv, µB)

∆² (46)

while we set ε = ∆. Similarly, ∂²F/∂µ_B∂µ_v and ∂²F/∂µ²_B could be derived easily by the exchange of symbols µ_B and µ_v.

Note that Newton-Quasi method aims at unconstrained optimization problem while our problem has nonnegative constraint and two functional constraints. First, the constraint sµ^∗_B > λ > 0 often holds because the queueing system would boil up and cause the rapid increase of Ls (F ) if the stability condition is violated. Some adjustments steps are provided to handle the following three situations:

Case 1. µ^∗_B> 0, µ^∗_v< 0, set µ^∗_v= 0, re-apply Newton’s method to find µ^∗_B. Case 2. µ^∗_B≥ µ^∗_v≥ 0, the optimization solution is feasible.

Case 3. µ^∗_v> µ^∗_B> λ/s, set µv= µB= t, re-apply Newton’s method to find t^∗. For Cases 1 and 3, we simplify the original problem to the approximate optimization problem with single variable, (i.e., discuss the optimal solution at binding equation).

As the example of IP access network mentioned in Section 1, for an access router with mean arrival rate λ = 5 customers/minute, mean working vacation time 1/η = 2 minutes, and the following cost parameters (per minute):

Ch= $30/customer, CB= $180/customer, CV = $45/customer, CI = $15/server, Cs= $30/unit.

The Newton-Quasi method is employed to deal with this wavelength reconfiguration problem and some numerical illustrations are provided. Table 1 presents the approximated optimal solutions and the corresponding system performance measures under various numbers of roving wavelengths (servers). From Table 1, it is observed that the minimum cost is 566.347 if the IP access network implements exhaustive schedule (assigns all tunable wavelengths to the next router when the queue becomes empty). The minimum expected cost and the optimal solution are 416.863 and (ˆs, ˆµv, ˆµB) = (3, 3.70716, 3.70716) respectively. That is, for this access router, providing three tunable wavelengths is appropriate for cost saving.

(15)

Table 1. The approximate value (ˆµ^∗_v, ˆµ^∗_B) under various values of servers (λ = 5.0, η = 0.5)

ˆ

s 1 2 3 4 5 6

ˆ

µ^∗_v 6.30991 4.39560 3.70716 3.36074 2.95602 2.61867 ˆ

µ^∗_B 7.27773 4.39560 3.70716 3.68261 3.96625 4.23769 F 566.347 427.706 416.863 433.770 454.700 474.614 Ls 2.75094 1.68139 1.50040 1.47136 1.49100 1.50565 E[LsN W ] 0.28704 0.36595 0.47377 0.55777 0.61370 0.65668 E[W V ] 0.71296 1.63405 2.52623 3.44223 4.38630 5.34332 E[L_sW V ] 0.46136 0.77155 0.87497 0.87658 0.86803 0.84668 E[IDW V ] 0.25160 0.86250 1.65126 2.56565 3.51827 4.49664

The result of Newton-Quasi algorithm in searching the optimal solution is shown in Table 2. In Table 2 (I) (the first stage), the optimal solution found by using Newton-Quasi method is Case 3 (µ^∗_v > µ^∗_B > λ/s). Therefore, some adjustment procedures (set µv = µB = t) are executed as shown in Table 2 (II). It appears from Table 2 that the Newton-Quasi method is working well and converges very fast.

Table 2. Newton-Quasi algorithm in searching the optimal solution (s = 3, λ = 5.0, η = 0.5, µ = [5.0/3 + 1.0, 5.0/3 + 2.0]^T)

(I) Newton-Quasi algorithm Stage I

Iteration 1 2 3 4 5 6

ˆ

µ^∗_v 2.66666 3.55229 3.88561 3.92150 3.92185 3.92185 ˆ

µ^∗_B 3.66666 3.61194 3.55529 3.54783 3.54776 3.54776 F 427.458 417.317 416.598 416.591 416.591 416.591

∂F/∂ ˆµv -20.5784 -4.36434 -0.38011 -0.00367 −3 × 10⁻⁷ −8 × 10⁻¹⁴

∂F/∂ ˆµB -9.85432 -1.71118 -0.12899 -0.00125 −1 × 10⁻⁷ −2 × 10⁻¹⁴ Ls 1.98542 1.57542 1.47034 1.46023 1.46013 1.46013 E[L_sN W ] 0.61936 0.50486 0.47733 0.47479 0.47477 0.47477 E[W V ] 2.38064 2.49514 2.52267 2.52521 2.52523 2.52523 E[L_sW V ] 1.02338 0.89421 0.85005 0.84547 0.84543 0.84543 E[IDW V ] 1.35726 1.60094 1.67262 1.67974 1.67981 1.67981

(II) Newton-Quasi algorithm Stage II

Iteration 1 2 3 4 5

ˆ

µ^∗_v= ˆµ^∗_B= t^∗ 3.6921146 3.7070019 3.7071552 3.7071552 3.7071552 F 417.49514 416.86883 416.86299 416.86299 416.86299

∂F/∂t -8.2040178 -0.7069619 -0.0063172 −5 × 10⁻⁷ 3 × 10⁻¹³ Ls 1.5913926 1.5084666 1.5004011 1.5004011 1.5004011 E[LsN W ] 0.5156125 0.4774836 0.4738062 0.4737730 0.4737729 E[W V ] 2.4843875 2.5225164 2.5261938 2.5262271 2.5262271 E[LsW V ] 0.8937283 0.8767539 0.8749862 0.8749701 0.8749701 E[IDW V ] 1.5906592 1.6457625 1.6512076 1.6512569 1.6512569 Next, we investigate the effect of system parameters on the optimal solution and the system performance measures. For various values of λ and η, some numerical results of the optimal solution and the corresponding system performance measures

(16)

are shown in Table 3. From Table 3, it is observed that (i) F increases when λ or η become larger; because the raise of average arrival rate or average vacation rate would increase the queueing length and the number of customers in the system;

(ii) from the first three columns, s^∗is insensitive to the parameter η; (iii) the optimal value s^∗ increases when the arrival rate becomes larger; that is, more tunable wavelengths are required when the traffic intensity grows up; and (iv) µ^∗_v and µ^∗_B becomes larger as λ increases which is a reasonable result.

7. Concluding remarks. The infinite capacity M/M/s queue with multiple working vacations is studied (denoted by M/M/s/WV). The matrix geometric method works well for computing the steady-state probabilities in this paper. This paper generalizes the M/M/1/WV queue. Moreover, an efficient algorithm (Newton-Quasi algorithm) is developed for searching the approximate optimum values (ˆs^∗, ˆµ^∗_v, ˆµ^∗_B) that minimize the cost function.

Table 3. System performance measures of an M/M/s queue with multiple working vacations under optimal operating conditions.

(λ, η) (5.0, 0.3) (5.0, 0.6) (5.0, 0.9) (2.5, 0.5) (5.0, 0.5) (7.5, 0.5) ˆ

s^∗ 3 3 3 2 3 3

ˆ

µ^∗_v 3.57169 3.73321 3.12192 2.79076 3.70716 4.75075 ˆ

µ^∗_B 3.57169 3.76697 4.30445 2.79076 3.70716 4.75075 F 399.313 423.440 436.139 300.583 416.863 503.624 Ls 1.57691 1.47988 1.53670 1.12063 1.50040 1.87384 E[LsN W ] 0.37561 0.51109 0.60014 0.31588 0.47377 0.51854 E[W V ] 2.62439 2.48891 2.39986 1.68412 2.52623 2.48146 E[LsW V ] 1.02429 0.82362 0.77411 0.57993 0.87497 1.06016 E[IDW V ] 1.60010 1.66529 1.62575 1.10419 1.65126 1.42130

Acknowledgments. The authors thank the referees for their remarks and com- ments, which helped us to improve the clarity of the article.

REFERENCES

[1] Y. Baba,Analysis of a GI/M/1 queue with multiple working vacations, Operations Research Letters, 33 (2005), 201–209.

[2] A. D. Banik, U. C. Gupta and S. S. Pathak,On the GI/M/1/N queue with multiple working vacations-analytic analysis and computation, Applied Mathematical Modelling, 31 (2006), 1701–1710.

[3] R. L. Burden and J. Douglas, “Numerical Analysis,” 7^thEdition, Brooks/Cole, USA, 2001.

[4] U. Chatterjee and S. P. Mukherjee, GI/M/1 queue with server vacations, Journal of the Operational Research Society, 41 (1990), 83–87.

[5] E. K. P. Chong and S. H. Zak, “An Introduction to Optimization,” 2^ndEdition, Wiley, New York, 2001.

[6] B. T. Doshi,Queueing systems with vacations–a survey, Queueing Systems Theory Appl., 1 (1986), 29–66.

[7] S. W. Fuhrmann and R. B. Cooper, Stochastic decompositions in the M/G/1 queue with generalized vacations, Operations Research, 33 (1985), 1117–1129.

[8] F. Karaesmen and S. M. Gupta, The finite capacity GI/M/1 queue with server vacations, Journal of the Operational Research Society, 47 (1996), 817–828.

[9] T. Lee, The M/G/1/N queue with vacation and exhaustive service discipline, Operations Research, 32 (1984), 774–784.

(17)

[10] J.-H. Li, W.-Y. Liu and N.-S. Tian, Discrete time GI/Geo/1 queue with multiple working vacations Queueing Systems, 56 (2007), 53–63.

[11] C.-H. Lin and J.-C. Ke,Multi-server system with single working vacation, Applied Mathe- matical Modelling, 33 (2009), 2967–2977.

[12] W.-Y. Liu, X.-L. Xu and N.-S. Tian, Stochastic decomposition in the M/M/1 queue with working vacations, Operations Research Letters, 35 (2007), 595–600.

[13] M. F. Neuts, “Matrix-Geometric Solutions in Stochastic Models. An Algorithmic Approach,”

Johns Hopkins Series in the Mathematical Sciences, 2, Johns Hopkins University Press, Bal- timore, Md., 1981.

[14] L. D. Servi and S. G. Finn,M/M/1 queues with working vacations (M/M/1/WV), Perfor- mance Evaluation, 50 (2002), 41–52.

[15] H. Takagi, “Queueing Analysis: A Foundation of Performance Evaluation,” Vol. 1, Vacation and Priority Systems, Part 1, North-Holland Publishing Co., Amsterdam, 1991.

[16] N. Tian, D. Zhang and C. Cao, The GI/M/1 queue with exponential vacations, Queueing Systems Theory Appl., 5 (1989), 331–344.

[17] J. A. White, J. W. Schmidt and G. K. Benett, “Analysis of Queueing System,” Operations Re- search and Industrial Engineering, Academic Press [Harcourt Brace Jovanovich, Publishers], New York-London, 1975.

[18] D. A. Wu and H. Takagi,M/G/1 queue with multiple working vacations, Performance Eval- uation, 63 (2006), 654–681.

Received December 2009; 1st revision July 2010; final revision June 2011.

E-mail address: [email protected] E-mail address: [email protected] E-mail address: [email protected] E-mail address: [email protected]