On a batch matching system with impatient servers and boundedly rational customers

Applied Mathematics and Computation 354 (2019) 308–328

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

On a batch matching system with impatient servers and boundedly rational customers Xudong Chai a,∗, Liwei Liu a,∗, Baoxian Chang a, Tao Jiang b, Zhen Wang a a b

School of Science, Nanjing University of Science and Technology, Nanjing 210094, Jiangsu, China College of Economics and Management, Shandong University of Science and Technology, Qingdao 266590, Shandong, China

a r t i c l e

i n f o

Keywords: Batch matching system Impatient servers Boundedly rational customers Provider’s utility PSO

a b s t r a c t In this paper we consider a batch matching system with impatient servers and boundedly rational customers, which is essentially a double-ended batch service system. Each server is responsible for a batch of customers and serves them instantaneously. Especially, we assume that the servers have the characteristic of ”impatience”, which will lead to incomplete batch processing of each service. To derive the stationary distribution of system states, we combine two methods (called probability generating function method and G-matrix method), and some main performance measures are captured. In addition, we assume the customers are boundedly rational, that is, they are unable to estimate their sojourn times and utilities accurately. To deal, a logit model is employed to model this characteristic of the customers, and by which, we obtain the customers’ equilibrium joining fraction which is compared with the equilibrium joining probability of fully rational customers afterwards. Depending on whether the service is visible or invisible, we conduct some numerical analyses to show main parameters’ impact on equilibrium fraction, along with some intuitive explanations. Finally, to make the system viable economically, a function of service provider’s utility is developed and several numerical plots are presented to study the impact of several main parameters on the utility. Through the PSO (Particle Swarm Optimization) algorithm, we numerically obtain the maximum utility point within a certain range of multiple parameters. © 2019 Elsevier Inc. All rights reserved.

1. Introduction The batch matching system we research on is essentially a kind of so-called double-ended (or double-sided) queueing system, in which servers and customers arrive at both sides in the system. There are not a few practical examples which can be characterized by double-ended queueing model. For instance, the passenger-taxi system, a taxi rank at a railway station or airport. These matching system were first proposed by Kendall [1], who introduced the passenger-taxi model as a typical example. In his model, one side is a queue of passengers, and the other side is a queue of taxis, the two queues will never coexist since the matching process is instantaneous. And Kendall further illustrated quite a lot of service systems can be analyzed as double-ended queueing systems. Afterwards, Dobbie [2] considered a kind of this system with nonhomogeneous Poisson arrival of passengers and taxis based on Kendall’s model, and obtained transient probability of queue length. Then Kashyap [3,4] applied the method of supplementary variable to capture the expected queue lengths of passengers and taxis ∗

Corresponding authors. E-mail addresses: [email protected] (X. Chai), [email protected] (L. Liu).

https://doi.org/10.1016/j.amc.2019.02.004 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

309

for the system with finite buffers. Recently, the double-ended queue model was still a hot topic for many scholars. Crescenzo et al. [5] considered a double-ended system with catastrophes occurring at constant rate, and exponentially distributed repair times. They obtained both the transient and steady-state probability, apart from these, they also derived a heavytraffic approximation. Shi and Lian [6] studied a passenger-taxi problem with limited taxi waiting space, and analyzed the strategic behavior of passengers in both observable and unobservable cases. Wang et al. [7] considered a double-ended queueing system with gated policy, they also obtained the strategic behavior of passengers and social optimal strategies. Nowadays, the matching systems exit in a variety of fields, including parallel processing, communication protocols, database concurrency control, inventory management and so on. Take the communication for example, networks with synchronization nodes which can be considered as a matching queueing system have been discussed in various areas [8]. However, aforementioned papers are all discussed on the assumption of one-to-one service, in our model, the service mechanism of the servers is batch service, which implies that each server has a fixed service capacity. So far, the research on the batch service system has aroused the interest of many scholars. Neuts [9] seems to be the first one who proposed the concept of general bulk service rule, i.e., servers serve customers in batches and the batch size varies from a minimum value to a maximum value. Afterwards, some investigators further studied on queueing system with bulk service, such as the research of Chaudhry and Templeton [10]. In recent years, Jiang et al. [11,12] studied on a multi-server polling system with batch service of an unlimited size. Bountali and Economou [13] researched on a system with batch service, and considered customers’ strategic based on the information of system status. Bu et al. [14] did research on a limited clearing queueing model with non-persistent customers and retry mechanism, in fact, the so-called limited clearing queueing system also can be seen as a kind of batch service system. Taking the psychology of individuals into consideration, we now introduce impatience which was usually used in queueing models by plenty of literature. Barrer [15] firstly took customer’s impatience into account in queueing system. He focused on the system with constant impatience times and captured the steady state distribution of queue length. Altman and Yechaili [16] had analyzed a queueing system with impatient customers and server’s vacations, and derived some closedform results of M/M/1 and M/G/1 queueing system respectively. Certainly impatient behaviour is also admitted in the field of communication networks. Brandt and Brandt [17] considered a queueing system with s-servers and the calls may leave the system due to customers’ impatience. So far, scholars have studied the impatience of customers in the queueing system quite thoroughly. Nevertheless, little work has been done to research on the servers impatience, which is exactly a part of our research. Now we focus on customer’s behavior. In traditional queueing literature, the analysis of customer’s behavior was always based on the hypothesis that customers are fully rational. Naor [18] firstly proposed this hypothesis to research incorporate customer’s behavior. He assumed that customers are able to accurately estimate their expected sojourn time, as a result, the expected utility of joining the system. Then subsequent researchers continued his work based on this hypothesis of customer’s full rationality. Yechiali [19] extended Naors model to a GI/M/1 queueing system, he showed that customer joins the system if and only if the queue length is less than a specific number. However, motivated by the reality, it’s difficult for customers to be fully rational. Customers always estimate their sojourn time roughly based on the information they have obtained, and then decide whether to join the system or not. Up to now, researches on boundedly rational customers in queueing systems are not much more than traditional researches on fully rational customer. Su [20] considered the classical newsvendor model with a boundedly rational decision maker. Huang et al. [21] studied a canonical service model, and analyzed customers’ behaviors based on visible queue case and invisible queue case, respectively. In general, the practical point of view, the research on bounded rationality of customers makes sense. In this paper, the innovation points and main contributions are summarized as follows:

• In this batch matching system, we propose a hypothesis of servers’ impatience, which can also be considered as a service mechanism formulated by the service provider to maximize revenue. Under this mechanism, the servers are likely to serve the customers in incomplete batch once the number of customers in the batch exceeds a certain number. And the research on servers impatience in batch matching system hasn’t been done in open literature. • We focus on a kind of systems where the lengths of customers and servers are unobservable. And we divide the discussion into two parts based on the availability of the service (i.e. whether servers exist in the system or not). Then the generating function method is applied to solve the conditional steady-state probabilities of the states when the service is unavailable. The results of this part are pretty neat. To solve remaining conditional steady-state probabilities, we employ the G-matrix method whose calculation is far less than the R-matrix method’s. Finally, using the normalization condition, we capture all the steady-state probabilities of system states. • We take the bounded rationality of customers into account in this batch matching system, that is, customers’ estimates of their own utilities are biased. To solve, we assume the deviation of customer’s estimate of sojourn time follows a logistic distribution. Then we obtain the equilibrium fraction of the customers who will join the system based on a logit model. Subsequently we investigate the impact of customer’s bounded rationality, and meanwhile compare the value of the joining fraction with the joining probability of fully rational customer. • Based on whether the service status is visible or invisible, we study the relationship between the equilibrium joining fraction of customers and main parameters through some specific examples. The results show the local monotonicity of the equilibrium fraction of joining customers with respect to each parameter within given range.

310

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

• We also take the service provider’s utility per unit time into account. Based on whether the service status is visible or invisible, we study the impact of each main parameter on the service provider’s utility. Then the PSO algorithm is employed to search the feasible combination of main parameters which achieves the maximum of the service provider’s utility. It is found that a system with visible service always brings greater benefit to the service provider than a system with invisible service under the same parameters. The rest of the paper is structured as follows. In Section 2, we formulate the batch matching system as a two-dimensional Markov process, and transition rate diagram is given. Then according to the transition rate diagram, we construct an infinitesimal generator matrix. In Section 3, we obtain the necessary and sufficient conditions for the stability condition of the system. In Section 4, steady-state distribution is derived by the generator function method and the G-matrix method. In Section 5, some various measures of system performance are presented. Based on the sojourn time given in Section 5, equilibrium arrival rate of boundedly rational customers is discussed in Section 6. Furthermore, in Section 7, we construct a utility function of the service provider and the PSO algorithm is introduced and applied to the optimization problem of provider’s utility. Finally, Section 8 concludes the paper. 2. Model formulation We consider a batch matching queueing system with impatient servers and boundedly rational customers. In this model, servers and customers arrive at both sides with infinite waiting space. We assume potential primary customers arrive singly according to a Poisson process with rate λ > 0. And the arrival time intervals of servers are supposed to be independent and exponentially distributed with rate θ > 0. Note that whether they are servers or customers, they won’t abandon the system halfway. And the customers and servers are matched in accordance with a first-come-first-served (FCFS) discipline. Each server serves the customers in batches with service capacity N and take his customers out of the system immediately once the matching completes. We are interested in the impatience of servers who provide batch services, which is commonly involved in some queueing systems with batch services in real life. Such as the system of coach departure station, the drivers will not always wait the passengers to fill the coach before leaving. Such mechanism is defined as impatience in a broad sense. Actually, this mechanism of impatience is sometimes planned by the service provider who pursues for considerable revenue. Now we illustrate the mechanism of servers’ impatience in this model specifically. When the service is available (there is at least one server in the system) and the number of customers in incomplete batch exceeds a certain threshold denoted by M − 1 (0 < M < N), the server starts to be impatient. After an exponential time (which is called impatient time in this paper) with mean β1 , he immediately takes his batch away from the system regardless of whether his batch is full or not. If customers fill

the batch during the impatient time, the server definitely serve the complete batch immediately. Furthermore, the servers are homogeneous for granted and the customers realize the mechanism of servers’ impatience which can also be viewed as an obtainable information. Last but not least, it is assumed that the inter-arrival times of servers, the inter-arrival times of customers and the impatient times are mutually independent. In addition, we also research on the strategic behavior of the customers. We know that a customer regards his joining or balking dilemma with respect to the information he obtains. For ease of handling, It is necessarily to assume that customers are all risk neutral and strategic, i.e., their aim is to maximize their estimated profit without considering its variability. In coincident with the reality, we take the customers’ bounded rationality into account, which will lead to the inaccurate estimation of their sojourn time. To model this characteristic of the customers, a logistic distribution is employed to approximate the distribution of customers’ evaluated error. So that a logit model is built to capture the equilibrium joining fraction of customers. The discussion about the bounded rationality of customers is not covered here at length and will be conducted detailedly in Section 6. In our model, we focus on a kind of system where the queue length of customers is invisible to a new arriving customer. We discuss the strategic behavior of boundedly rational customers in two circumstances, invisible service and visible service. Under the former circumstance, new arriving customers are incapable of observing the service state before they join the system. That is to say, they are unable to judge whether there is a server in the system before his joining. Conversely, new arriving customers can observe the service state before joining under the latter circumstance. And compared to the situation where the service is unavailable, customers are more likely to join the system when the service is available. Based on the above statement and for convenience, we denote λ1 and λ2 as the effective arrival rate when the service is unavailable and the effective arrival rate when the service is available, respectively. For the case of invisible service, λ1 = λ2 . And for another case of visible service, λ1 < λ2 . We formulate the process as a special quasi birth-and-death(QBD) process. And define the state space of the queueing system at time t by {X(t), Y(t)}. X(t) denotes the number of batches when the service is not available (X(t) ≥ 0) or the number of servers in the system when the service is available (X(t) < 0) at time t. And Y(t) denotes the number of customers in incomplete batch at time t. Thus {X(t), Y(t)} is a bivariate Markov process, the corresponding state space is

 = {0, ±1, ±2, . . .} × {0, 1, 2, . . . , N − 1}. In order to study the steady-state behavior of the system, let X ≡ limt→∞ X (t ), Y ≡ limt→∞ Y (t ), and the steady-state probability distribution pi, j = P {X = i, Y = j}, i = 0, ±1, ±2, . . . ± ∞, j = 0, 1, 2, . . . , N − 1. The transition rate diagram of this

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

311

Fig. 1. Transition rate diagram of the batch matching system with impatient servers and boundedly rational customers.

system is represented in Fig. 1. To obtain the systems steady-state probabilities, system’s states are arranged in the following order:

{. . . , S2 , S1 , S0 , S−1 , S−2 , . . .},

Si = {(i, 0 ), (i, 1 ), (i, 2 ), . . . , (i, N − 1 )}.

Thus we construct an infinitesimal generator matrix denoted by Q as shown below:

⎛ ⎜ ⎜ ⎜ ⎜ Q =⎜ ⎜ ⎜ ⎜ ⎝

..

.

..

. B

⎞

..

. C B

A0 C B

A0 C A2

A0 A1 A2

A0 A1 .. .

A0 .. .

..

⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎠ .

where the submatrices A0 , A1 , A2 , B and C, each of order N × N, are given by ⎛ ⎞ 0 ⎜ .. ⎟ ⎞ ⎛ ⎜ . ⎟ θ ⎜ ⎟ ⎜ 0 ⎟ θ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ A0 = A 2 = ⎜ β ⎟, .. ⎝ ⎠, ⎜ . ⎟ . ⎜ . ⎟ ⎜ . ⎟ θ

⎝ β ⎠ λ2 + β

(1)

312

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

⎛

−(λ2 + θ )

⎜ ⎜ ⎜ ⎜ ⎜ ⎜ A1 = ⎜ ⎜ ⎜ ⎜ ⎜ ⎝

⎛ ⎞

0 ⎜ .. ⎟ B = ⎜ . ⎟, ⎝ ⎠ 0

⎞

λ2 −(λ2 + θ )

λ2 ..

..

.

. −(λ2 + θ )

λ2

−(λ2 + θ + β )

λ2 ..

⎛

−(λ1 + θ )

⎜ ⎝

..

..

. −(λ2 + θ + β )

⎞

λ1

C=⎜

.

⎟ ⎟. ⎠

..

.

. −(λ1 + θ )

λ1

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎟ ⎠ λ2 −(λ2 + θ + β )

λ1

−(λ1 + θ ) Note that all elements in blank of the above matrices are 0. Given these submatrices, then we consider the condition for this system’s stability in the next section. 3. Analysis of the model Theorem 3.1. The bivariate Markov process {X(t), Y(t), t ≥ 0} is stable if and only if

λ1 N

<θ <

λ2 

M + λβ2 − λβ2

λ2

N−M ,

(2)

λ2 +β

Proof. To facilitate subsequent analysis, let

⎛

−λ2

⎜ ⎜ ⎜ ⎜ ⎜ A ≡ A0 + A1 + A2 = ⎜ ⎜ β ⎜ ⎜ .. ⎜ . ⎝ β λ2 + β ⎛

−λ1

⎜ ⎜ D ≡ B + C + A0 = ⎜ ⎜ ⎝

λ1 −λ1

⎞

λ2 −λ2

λ2 ..

.

..

. −λ2

λ2

−(λ2 + β )

λ2 ..

.

..

. −(λ2 + β )

λ2

⎟ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎟ ⎠

(3)

−(λ2 + β )

⎞ λ1 ..

.

λ1

..

. −λ1

⎟ ⎟ ⎟, ⎟ λ1 ⎠ −λ1

(4)

and eN ≡ (1, 1, . . . , 1 )T , a N-dimensional column vector of ones. According to the statement of this model, since the space capacities of both sides are all infinite, so in order to ensure the stability of the system, we can’t let the state extend indefinitely to either end. Referring to Neuts [22], the stability condition can be presented as

 1 A0 eN < π  1 A2 eN , π

(5)

 2 BeN < π  2 A0 eN . π

(6)

 1 is the invariant probability vector of Inequality (5) ensures the number of servers in the system doesn’t go to infinity, π A, and satisfies the linear equations:

 1 A = 0, π  1 eN = 1. π Solving the linear system of Eq. (7) yields

(7)

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328





313





  2  λ2 λ2 λ2 λ2 λ2 λ2 N−M −1 N−M  1 = M+ π − × 1, 1, . . . , 1 , , ,..., , β β λ2 + β λ2 + β λ2 + β λ2 + β ↑ Mth

 1 into the inequality (5), through some algebraic manipulation, we obtain substituting π

θ<

λ2

(8)

M + λβ2 − λβ2 ( λ λ+2 β )N−M 2

 2 is the Analogously, inequality (6) ensures the number of customers in the system doesn’t grow indefinitely, and π invariant probability vector of D, i.e., the unique solution of the following equations:

 2 D = 0, π  2 eN = 1. π

(9)

2 = Thus, we obtain π

λ1 N

1 N ( 1, 1, . . . , 1 ).

 2 into inequality (6), it follows that Substituting π

< θ.

(10)

Combining with (8) and (10) yields

λ1 N

<θ <

λ2

, λ λ λ2 N−M M+ β2 − β2 ( λ + ) β 2

which is exactly the necessary and sufficient condition

for the system to be stable and ensures the number of either customers or servers will not increase infinitely. Therefore, the steady-state of the system is guaranteed.  λ

λ

λ

Remark 1. In Theorem 3.1, since M + β2 − β2 ( λ +2 β )N−M = M + 2 λ2 λ1 that N−M > N .

N−M i=1

i

( λ λ+2 β ) < M + N − M = N, and λ2 ≥ λ1 , it’s obvious 2

λ λ2 λ M+ β2 − β2 ( λ + ) 2 β

Remark 2. When the value of β tends to infinity, that is, once the server finds the number of customers he matches exceeds M − 1, this batch will be served and taken away from the system immediately. Therefore, in this case, the steady-state λ λ condition becomes N1 < θ < M2 . 4. Steady-state distribution In this section, our discussion will concentrate on the steady-state distribution of the process {X(t), Y(t), t ≥ 0}. To this end, some definitions are provided as below:

pk ≡ ( pk,0 , pk,1 , . . . , pk,N−1 ),

k = 0, ±1, ±2, . . . ,

p ≡ (. . . , p2 , p1 , p0 , p−1 , p−2 , . . . ), Ls (t ) ≡ the number of servers in the system at time t, L(t ) ≡ the number of customers in the system at time t. To consider the steady-state of this system, we may define Ls ≡ limt→∞ Ls (t ) and L ≡ limt→∞ L(t ) at steady-state. Specifically, according to the definitions of X and Y, L can be expressed as



L=

NX + Y,

X ≥ 0,

Y,

X < 0,

which will be utilized in latter computation of conditional expected queue size. For convenience in subsequent calculations, we denote lk = P {L = k, X ≥ 0} as the joint probability that the service is unavailable and there are k customers in the system. And the corresponding (partial) probability generating function is

l (z ) =

∞ 

z n ln .

n=0

Now we conduct our work on the calculating of the steady-state probability distribution pi,j . Due to the special structure of the Q-matrix (1), we split the calculation into two steps. Firstly, we use probability-generating function method to obtain the stationary probabilities of the upper half states of the transition rate diagram shown in Fig. 1. And the results are all expressed in terms of p0,0 (or l0 ). According to the transition rate diagram and utilizing the notation lk , we write the balance equation of the states where service is unavailable as below:

314

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

(λ1 + θ )l0 = θ lN + β

N−2 

p−1, j + (λ2 + β ) p−1,N−1 ,

(11)

j=M

(λ1 + θ )lk = θ lk+N + λ1 lk−1 ,

k = 1, 2, 3, . . . .

(12)

Taking the sum of (11) and (12) with respect to k derives

θ

N−1 

li = β

i=0

N−2 

p−1, j + (λ2 + β ) p−1,N−1 .

(13)

j=M

Multiplying (11) and (12) by zk and summing over k derives

θ

(λ1 + θ )l (z ) =

z

(l (z ) − N

N−1 

li zi ) + λ1 zl (z ) + β

N−2 

i=0

p−1, j + (λ2 + β ) p−1,N−1 ,

(14)

j=M

combined with (13), it follows that



l (z ) =



N−1 N i θ N−1 i=0 li z − θ i=0 li z . (λ1 + θ )zN − λ1 zN+1 − θ

(15)

Obviously, z = 1 is a zero of both the numerator and the denominator in (15). To simplify the expression for l(z), we next prove the denominator of the right-hand side of (15) has only one zero outside the disk |z| ≤ 1. Thus we provide the following Lemma 4.1. Lemma 4.1. (λ1 +θ )zN − λ1 zN+1 − θ = 0 has only one root outside the disk |z| ≤ 1. Proof. Let u = 1z , so that (λ1 +θ )zN − λ1 zN+1 − θ translates into ( 1u )N+1 ((λ1 + θ )u − λ1 − θ uN+1 ). We now denote f (u ) = (λ1 + θ )u and g(u ) = −(λ1 + θ uN+1 ). It’s clearly that f(u) has only one zero u = 0. For any δ in (0,1), we focus on |u| = 1 − δ, we have

|g(u )| ≤ λ1 + θ (1 − δ )N+1 = λ1 + θ (1 − (N + 1 )δ ) + o(δ ) = λ1 + θ − θ (N + 1 )δ + o(δ ), | f (u )| = (λ1 + θ )(1 − δ ) = λ1 + θ − (λ1 + θ )δ > λ1 + θ − θ (N + 1 )δ, thus,

|g(u )| < | f (u )|, f or every |u| = 1. By Rouch’s theorem, f(u) and f (u ) + g(u ) have the same number of zeros in the disk |u| < 1, that is to say, f (u ) + g(u ) = (λ1 + θ )u − λ1 − θ uN+1 has only one zero in the disk |u| < 1, we denote this zero by u0 . In other words, it implies that (λ1 +θ )zN − λ1 zN+1 − θ = 0 has only one root outside the disk |z| ≤ 1, and this zero is u10 .  According to the expression for the numerator and denominator of l(z), it’s not difficult to find that the numerator has N zeros while denominator has N + 1 zeros. Utilizing Lemma 4.1, the denominator has only one zero outside the disk |z| ≤ 1. And when |z| ≤ 1, the partial generating function l(z) converges. Thus l(z) is an analytic function in the disk |z| ≤ 1, which leads to

(λ1 +θ )zN − λ1 zN+1 − θ = 0 z−

1 u0



= C0

θ

N−1 

li z − θ

i=0

N

N−1 

li z

i

,

(16)

i=0

where C0 is a constant. Combining (16) with the expression of l(z) shown in (15), we obtain

l (z ) =

1 C0 (z −

1 u0

)

.

(17)

u

u

Since l (0 ) = − C0 = l0 = p0,0 , the constant C0 can be expressed by p0,0 , that is C0 = − p 0 . Substituting it back into (17) 0,0

0

derives

l (z ) =

∞

∞

j=0

j=0

  1 p0,0 = ( u 0 z ) j p0,0 = ( u 0 j p0,0 ) z j , 1 − u0 z

(18)

which implies that the probabilities of the states where the service is unavailable are all expressed as

lk = uk0 l0 ,

k = 0, 1, 2, . . .

(19)

or +j pi, j = uNi p0,0 , 0

i = 0, 1, 2, . . . ,

j = 0, 1, . . . , N − 1.

(20)

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

315

Reemphasize that u0 is the unique root of (λ1 + θ )u − λ1 − θ uN+1 = 0 inside the disk |u| < 1. Since the expression of the probability set {pi,j |i ≥ 0} in terms of p0,0 has been captured, the probability vector p0 is presented as

p0 = (1, u0 , u20 , . . . , uN−1 ) p0,0 , 0

(21)

which is crucial for the subsequent calculation of the remaining steady-state probabilities. Now we conduct the second step of the calculation to calculate the remaining steady-state probability distribution {pi,j |i < 0}. To this end, we apply the method of matrix geometric analysis (see, e.g.,Neuts [22], Latouche and Ramaswami  [23]). We firstly denote  as the state Sk , which can be viewed as the initial state. Thus we reconstruct the state space k>0

as ∗ = {, S0 , S−1 , S−2 , . . .}, and the corresponding infinitesimal generator matrix is presented as below

⎛



B0 ⎜B2

S0 Q ∗ = S−1 S−2 .. .

⎜ ⎜ ⎜ ⎝

⎞

B1 C A2

A0 A1 A2

A0 A1 .. .

A0 .. .

..

⎟ ⎟ ⎟. ⎟ ⎠

(22)

.

It should be noted that, since we have obtained p0 in terms of p0,0 , the concrete expressions for B0 , B1 and B2 are nonsignificant for the subsequent  analysis. According to matrix geometric method, the conventional way is R-matrix method. The so called rate matrix R = ri, j (of order N × N)) satisfies

A0 + RA1 + R2 A2 = 0.

(23)

If we obtain the the minimal non-negative solution R to the matrix Eq. (23), the remaining steady-state probabilities can be solved, and have the matrix-geometric form

p−k = p0 Rk ,

k = 1, 2, . . . .

(24)

However, all the elements in R are uncertain and hard to obtain. Thus, we offer an alternative method called G-matrix method to calculate {pi,j |i < 0}. The explanation and significance of the matrix G can be referred to [23,24]. Similar to the R-matrix method, we should look for the minimal non-negative solution of

A2 + A1 G + A0 G2 = 0

(25)

Because only the first column of matrix A2 is non-zero, so is the matrix G. Therefore, we can assume that G has the same structure as

⎛ ⎜ ⎝

G=⎜

g0 g1 .. .

⎞

⎟ ⎟. ⎠

gN−1 Substituting G into Eq. (25) and expanding the equation yields

−(λ2 + θ )gi + λ2 gi+1 + θ g0 g1 = 0,

i = 0, 1, . . . , M − 1,

β − (λ2 + θ + β )g j + λ2 g j+1 + θ g0 g j = 0,

(26)

j = M, M + 1, . . . , N − 2,

(27)

λ2 + β − (λ2 + θ + β )gN−1 + θ g0 gN−1 = 0.

(28)

By (26) and (27), we have

gi =

λ2 + θ − θ g0 λ2



gj =

i

λ2 + θ − θ g0 λ2

g0 ,

M

i = 0, 1, . . . , M,

β g0 − θ + β − θ g0

(29)



λ2 + θ + β − θ g0 λ2

 j−M +

β , θ + β − θ g0

j = M, M + 1, . . . , N − 1, (30)

316

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

Substituting gN−1 which derived by (30) into (28), yields



λ2 + θ − θ g0 λ2

M

β g0 − θ + β − θ g0



λ2 + θ + β − θ g0 λ2

N−M−1 +

β λ2 + β = . θ + β − θ g0 λ2 + θ + β − θ g0

(31)

According to the matrix geometric theory, there must exist a g0 ∈ (0, 1) and thus further we obtain all elements of matrix G. The G-matrix method provides another way to obtain the limiting probabilities, according to [25], the relation between G and R can be expressed as below:

R = A0 (−A1 − A0 G )−1 .

(32)

So that it gives us an indirect approach to obtain R explicitly in our model. Since the matrix R has been captured explicitly, combining (21) with (24), we finally got the expression for p−k in terms of p0,0 . Remark 3. Basically, although the G-matrix and R-matrix are equivalent, they have own advantage and limitation. Compared with the G-matrix, the R-matrix method can provide a more succinct representation of the steady-state probability and the sufficient and necessary condition on the ergodicity of Markov chain. However, it should be noted that there is no bound on the elements of R-matrix, while all elements of G-matrix must be in the range of (0, 1), so the computation of G is numerically stable and simple. Based on the above analysis, we have obtained the expressions for pi,j in terms of p0,0 , then use normalization condition to derive p0,0 . We summarize our results in following theorem. Theorem 4.1. At steady state, for a batch matching system with impatient servers and boundedly rational customers, the stationary probabilities pi,j can be expressed as



pi, j =

+j uNi p0,0 , 0

i = 0, 1, 2, . . . ,

p0 R−iI j+1 ,

i = −1, −2, . . . , j = 0, 1, . . . , N − 1,

j = 0, 1, . . . , N − 1,

where Ik = (0, 0, . . . , 1 , 0, . . . , 0 )T , and p0,0 = ↑ kth

1 1−u0

(33)

1 . +(1,u0 ,u20 ,...,uN−1 )R(I−R )−1 eN 0

Proof. The theorem has been proved in the above.



It should be noted that the effective arrival rates λ1 and λ2 are still unknown, and will be discussed in Section 6 under the consideration of the bounded rationality of customers. So λ1 and λ2 can be viewed as two symbols before Section 6. 5. Performance measures For a new arrival, we divide his sojourn time into two periods. The first period is the time the customer should wait his server to match with him. We call this period the matching time. And the second period lasts from the epoch he matches with his server to his departure. Due to the difference in arriving time of customers, the two periods of time could be 0. In this section, we will calculate the expected queueing length of customers and servers first. Then we aim to obtain the expected matching time and the expected sojourn time, among which, some performance measures are analyzed based on whether the service state is visible or invisible. More concretely, the results are presented as follows: The probability that the service is unavailable is

P {X ≥ 0 } =

p0,0 , 1 − u0

(34)

and the probability that the service is available is

P {X < 0} = p0 R(I − R )−1 eN



or = 1 −



p0,0 . 1 − u0

(35)

When the service is unavailable, the conditional expected queueing length of customers is

E ( L|X ≥ 0 ) =

∞  k=0

∞

 lk u0 k= (1 − u0 )uk0 k = . P {X ≥ 0} 1 − u0

(36)

k=0

When the service is available, the conditional expected queueing size of customers is

E ( L|X < 0 ) =

∞





 p0,0 1 −1 p0 RivN = p0 R(I − R ) vN / 1 − , P {X < 0} 1 − u0 i=1

where vN = (0, 1, 2, . . . , N − 1 )T .

(37)

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

And the expected queueing length of servers is

E ( Ls ) =

∞ 

i p−i eN =

i=1

∞ 



i p−1 R

i−1

eN = p−1

i=1

∞ 

317

iR

i−1

eN = p−1 (I − R )−2 eN (38)

i=1

= p0 R(I − R )−2 eN . To discuss the expected sojourn time of a new arriving customer, we first calculate two conditional expected sojourn times of the arrival, and summarize them in following lemmas. Lemma 5.1. When the system is just at the state {(i, k )|i = −1, −2, . . . , k = 0, 1, . . . , N − 1}, for a new arrival, his expected sojourn time denoted by S(k) can be expressed as

S (k ) =

⎧   λ N−M ⎪ 2 ⎪ M−1−k + 1 1− , ⎪ ⎨ λ2 β β + λ2   λ N−1−k ⎪ ⎪ ⎪ 2 ⎩ 1 1− , β β + λ2

k = 0, 1, . . . , M − 1, (39) k = M, M + 1, . . . , N − 1.

Proof. By routine analysis, we know that S(k) satisfies the following linear system of equations

S (0 ) = S (1 ) =

1

λ2 1

λ2

+ S ( 1 ),

(40)

+ S ( 2 ),

(41)

.. . S (M − 2 ) = S (M − 1 ) = S (M ) =

1

+ S ( M − 1 ),

λ2

λ2 S ( M ), β + λ2

(43)

λ2 S ( M + 1 ), β + λ2

(44)

1

+

β + λ2

1

β + λ2

+

(42)

.. . S (N − 3 ) = S (N − 2 ) =

1

β + λ2 1

β + λ2

+

λ2 S ( N − 2 ), β + λ2

(45)

,

(46)

S ( N − 1 ) = 0.

(47)

Substituting (46) into (45) yields S(N − 3 ), then substituting S(N − 3 ) into the previous equation yields S(N − 2 ). Step by step, after some algebraic manipulation we can capture all S(k ), k = 0, 1, . . . , N − 1, and the results are concluded as (39).  Lemma 5.2. When the system is at the state {(i, j )|i = 0, 1, 2, . . . , j = 0, 1, . . . , N − 1}, for a new arrival, his expected matching time is S1 (i, j ) = i+1 θ . And his expected sojourn time denoted by S(i, j) can be expressed as follows

S(i, j ) =

i+1

θ



+

θ

λ1 + θ

i+1 N− j−1

 k=0

i+k k



λ1

λ1 + θ

k

S ( j + k ).

(48)

Proof. The first half of the theorem is apparent, we focus on the second half the theorem. We define two independent random variables W1 (i) and Z (N − j − 1 ), which are Erlang-(i + 1 ) and Erlang-(N − j − 1 ) with rates θ and λ1 , corresponding to the time to wait the server which the customer matches with and the completion time of the current incomplete batch, respectively. Consider, the target customer enter the system when the system is at a state (i, j), we discuss the problem in two cases:

318

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

Case 1: If W1 (i ) ≥ Z (N − j − 1 ), that is to say, the completion time of the current incomplete batch is less than the time the target customer to wait the server, so that the expected sojourn time is S1 (i, j ) = i+1 θ . Case 2: If W1 (i ) < Z (N − j − 1 ), the server which serve the batch the target customer in arrive at the system before the batch completes. However, the server will be impatient when the number of customers in incomplete batch is greater than M − 1, and the number of new arrival customers during Erlang-(i + 1 ) time W1 (i) denoted by N(W1 (i)) is unknown. So that according to the distribution of N(W1 (i)), we can obtain the sojourn time by the total probability formula. For example, N (W1 (i )) = k, k < N − j − 1, the target customer’s sojourn time is equal to S1 (i, j ) + S( j + k ). Based on the above analysis of two cases, we have

S(i, j ) = =

i+1

θ



N− j−2

P {W1 (i ) > Z (N − j − 1 )} +



i+1

θ

k=0

i+1

θ

i+1 

 + S( j + k ) P {N (W1 (i )) = k}

N− j−2

P {N (W1 (i )) ≥ N − j − 1} +



θ

P {N (W1 (i )) = k}

k=0

N− j−2

+

S( j + k )P {N (W1 (i )) = k}

k=0

=

=

i+1

θ

i+1

θ



N− j−2

+

S( j + k )P {N (W1 (i )) = k}

k=0



N− j−1

+

S( j + k )P {N (W1 (i )) = k}.

k=0

The expression for S(j) is valid in (39), and P {N (W1 (i )) = k} can be computed as below:

P {N (W1 (i )) = k} =



∞

0

 =

∞

P {N (t ) = k}

(λ1 t )k

=

k!

0



i+k i

e−λt

λ1

θ i+1t i i!

θ i+1t i i!

e−θ t dt

e−θ t dt

k

λ1 + θ

θ

i+1 .

λ1 + θ

Substituting the expression for P {N (W1 (i )) = k} into S(i, j) derives

S(i, j ) =

i+1

θ



+

θ

i+1 N− j−1



λ1 + θ



i+k k

k=0

λ1

k

λ1 + θ

S ( j + k ).

 Utilizing the above conclusions, we calculate the expected matching time and the expected sojourn time under two conditions based on whether the system is visible or invisible. I. The service state is visible Theorem 5.1. For the batch matching system with impatient servers and boundedly rational customers, the expected matching time of a new arriving customer upon seeing unavailable service is

E (Wm |X ≥ 0 ) =

∞ N−1

  k + 1 kN+ j 1 u P0,0 = P {X ≥ 0} θ 0 k=0 j=0

1

θ (1 − uN0 )

.

(49)

We can also apply Little’s law to capture the result. According to (36), we have

E (Wm |X ≥ 0 ) =

E ( L|X ≥ 0 )

λ1

=

u0

λ1 (1 − u0 )

.

(50)

Combined with (λ1 + θ )u0 − λ1 − θ u0 N+1 = 0, it’s not difficult to find that the result (50) is equivalent to (49). Remark 4. Here is a intuitive explanation for why the expected matching time can be calculated by applying Little’s law as above. As we know that the matching time is the time length for a new arriving customer has to wait for his server. Once he matches the server, he leaves the subsystem where the customers wait to be matched, and corresponding queueing length decreases. Therefore, Little’s law is suitable in this situation.

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

319

Theorem 5.2. For the batch matching system with impatient servers and boundedly rational customers, the expected sojourn time of an arrival upon seeing the service is available is

E (W |X < 0 ) =

∞

 p0,0 1  N = p0 R(I − R )−1 w  N / (1 − p0 Ri w ), P {X < 0} 1 − u0

(51)

i=1

 N = (S(0 ), S(1 ), . . . , S(N − 1 ))T , and the expected sojourn time of an arrival upon seeing the service is unavailable is where w

E (W |X ≥ 0 ) =

 θ u0 (1 − u0 ) N−1 (k + 1 )uk0 S(k ) λ1 k=0

1 + θ (1 − uN0 )

(52)

Proof. The expression of the expected sojourn time of an arrival upon seeing the service is available in (51) is obvious. We focus on the expected sojourn time of an arrival upon seeing the service is unavailable, combining with (48), we have

E (W |X ≥ 0 ) =

∞ N−1

 1 +j S(i, j )uNi p0,0 = E (Wm |X ≥ 0 ) 0 P {X ≥ 0} i=0 j=0





j−1 ∞ N−1   N−  1 i+k + k P {X ≥ 0} i=0 j=0

=

i+1

θ

λ1 + θ



k=0

j−1 ∞ N−1  N−   i+k 1 + ( 1 − u ) 0 N k θ ( 1 − u0 ) j=0 i=0 k=0

θ

k

λ1

λ1 + θ i+1

+j S( j + k )uNi p0,0 0

uNi 0

λ1 + θ

λ1

λ1 + θ

k S( j + k )u0j .

Then let

Tk =

∞ 



i=0

i+k k



i+1

θ

uNi 0 =

λ1 + θ



θ

λ1 + θ



∞  i+k k i=0

θ uN0 λ1 + θ

i ,

combining with the equation (λ1 + θ )u0 − λ1 − θ uN+1 = 0, the sequence {Tk , k = 0, 1, 2, . . .} has the following recursion 0

Tk =

(λ1 + θ )u0 Tk−1 . λ1

Therefore, we obtain



Tk =

(λ1 + θ )u0 λ1

here we note (



(53)



k

T0 =

(λ1 + θ )u0 λ1

 k

θ u0 , λ1

(54)

0 ) = 1. 0

So that j−1 N−1  N−  1 E (W |X ≥ 0 ) = + ( 1 − u ) 0 θ (1 − uN0 ) j=0 k=0

=

=

1

θ (1 − uN0 ) 1

θ (1 − uN0 )

+

+



(λ1 + θ )u0 λ1

θ u0 ( 1 − u0 ) λ1

j−1 N−1  N− 

θ u0 ( 1 − u0 ) λ1

N−1 

j=0

k

k θ u0 λ1 S( j + k )u0j λ1 λ1 + θ

u0j+k S( j + k )

k=0

(k + 1 )uk0 S(k )

k=0

 II. The service state is invisible Theorem 5.3. For the batch matching system with impatient servers and boundedly rational customers, when the service state is invisible, λ1 = λ2 , the expected matching time of an arrival is

E (Wm ) = E (Wm |X ≥ 0 )P {X ≥ 0} + E (Wm |X < 0 )P {X < 0} = And the expected sojourn time of an arrival is

p0,0 . θ (1 − u0 )(1 − uN0 )

(55)

320

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

E (W ) = E (W |X ≥ 0 )P {X ≥ 0} + E (W |X < 0 )P {X < 0} =

p0,0 + θ (1 − u0 )(1 − uN0 )

 θ u0 p0,0 N−1  N, (k + 1 )uk0 S(k ) + p0 R(I − R )−1 w λ1 k=0

(56)

It’s certain that the facilitators hope the servers’ sojourn time is as short as possible for the sake of high profit. So the expected sojourn time of servers is also an important performance indicator, and the expression for which is given in the following theorem. Theorem 5.4. For the batch matching system with impatient servers and boundedly rational customers, the expected sojourn time of servers is

E (Ws ) =

1

θ

p0 R(I − R )−2 eN

(57)

Proof. According to (38) and Little’s law, the result is obtained immediately.



6. Equilibrium arrival rate of boundedly rational customers In many literature, the authors always assume the customer is perfectly rational, which means they have the capability to calculate his expected sojourn time and utility accurately. However, it’s not realistic. In our model, we assume that the customer is boundedly rational, i.e., they are incapable of estimating their sojourn time and utility accurately. To model this characteristic of the customers, we add a error term to the customer’s sojourn time. To this end, we draw on the thoughts of Huang et al. [21] who models customer queueing choices with a logit model to consider a regular service system but with boundedly rational customers. We now apply Huang’s model to solve the equilibrium arrival rate of bounded rational customers. In our model, we define the actual utility function of customers as U = R0 − p − CE (W ), and the estimated utility function of customers as U1 = R0 − p − C (E (W ) + ε ), where R0 > 0 is customer’s reward on completion of service, p is the price of the service, C > 0 is the average sojourn time cost per unit of time, W is the expected sojourn time, and ε is the random error term which represents the noisy sojourn time estimation. And we know the joining behavior of the customer depends on whether the estimated utility U1 is positive or not. And if the customer balks, he will get zero utility. To model the error term of customers’ estimation, we assume the random error term ε follows a logistic distribution

F (x ) =

ex/η , 1 + ex/η

(58)

where η > 0. Referring to [26], the parameter η measures the error level of customer expected sojourn time estimation (as a result, customers’ expected utility estimation). One way of explaining this is that, let η → 0, it’s not difficult to understand, customers tend to be perfectly rational, whether they join or not depends entirely on their values of actual utility. While η → ∞, customers tend to be completely irrational, so that they join or balk with equal fraction. Thus, based on (58), we can obtain the fraction of customers who will join the system as follows

  U φ=P ε≤ = C

eU/Cη 1 + eU/Cη

(59)

It should be emphasized that φ denotes the fraction of customers who will join the system, doesn’t mean customers choose to play mixed strategies. For the boundedly rational customers, there are only two kinds of strategic behaviors, either joining or not. According to (59), due to the bounded rationality of customers, the fraction of customers who will join the system will not be 0. It’s not difficult to find that the function of the customer’s utility shown in (59) is implicit and continuous with respect to φ . Thus, (59) leads to a fixed-point problem. Due to the complexity of the function on the right-hand side of (59), it is unrealistic to calculate theoretically and furthermore, the number of fixed points and the convergence near the fixed point is not clear. Therefore the conventional fixed point iteration algorithm cannot be applied. To obtain the roots of (59), we adopt the method of global search within the interval of (0,1) to search the values of φ . To this end, we use a numerical example to explore the roots of Eq. (59), consider the system with invisible service and N = 10, M = 5, λ = 12, θ = 1, β = 40, R0 = 20, p = 10, C = 5, η = 1. Where the red curve and blue line in Fig. 2 represent the functions φ1 =

eU (φ0 )/C η 1+eU (φ0 )/C η

and φ1 = φ0 , respectively. Intuitively,

for low values of the effective arrival rate, the time to complete a batch is long, so that few customers have an incentive to join. As the effective arrival rate increases, the time to complete the batch decreases, so customers become more willing to join. However, for high values of the effective arrival rate, the system becomes congested. These imply the mixed ATC-FTC nature appears in our model, and due to which, the function φ1 =

eU (φ0 )/C η 1+eU (φ0 )/C η

is always monotonically increasing in φ 0 firstly

and then decreasing shown as the red curve in Fig. 2. Remark that the red curve represents the function under the steadystate, so the solution should satisfy the steady-state condition. In this numerical example with given parameters, there is only one root within the range of the steady-state, that is the fixed point 3 shown in the Fig. 2. Moreover, given the sensible

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

321

1 0.9 0.8 0.7

fixed point 3 =0.645 0

1

0.6 0.5 0.4 0.3 0.2

fixed point 1 =0 0 fixed point 2 =0.052 0

0.1 0 0

0.2

0.4

0.6

0.8

1

0 Fig. 2. An example of searching fixed points with given parameters when service is invisible.

values of the parameters, there generally exit only one root within the range of the steady-state. And it’s rare that there are no roots for (59) within the range of the steady-state, so we can ignore these extreme situations. Certainly, at steady-state, there may be two roots when the value of M is pretty small relative to N, in which case we still select the higher value between the two roots. We may explain it intuitively in this way: we can view φ 1 as the fraction of subsequent arriving customers who will join the system. Virtually φ 1 is a reaction to φ 0 which relates to previous customers. And the relation of the values of φ 1 and φ 0 determines the trend of joining fraction. When φ 1 is greater than φ 0 , the joining fraction tends to increase, otherwise, tends to decrease. Owing to the mixed ATC-FTC nature, the function φ1 =

eU (φ0 )/C η 1+eU (φ0 )/C η

is always concave

within the range of the steady-state. Thus only the larger root is stable, because a slight disturbance will not cause moving away from the equilibrium based on the previous intuitive explanation for the relation of φ 1 and φ 0 . Therefore we select the larger root as the equilibrium joining fraction. Depending on whether the service state is visible or not, we calculate the equilibrium arrival rate, respectively. When the service state is invisible, λ1 = λ2 , we substitute (56) into (59). Here we denote φ (i) as the joining fraction when the service is invisible. And the specific algorithm steps are presented in Algorithm 1. Algorithm 1: Calculate the equilibrium joining fraction when service is invisible.

1 2 3

4 5

Input: R0 , p, C, η, λ, θ , β , ε , M, N Output: φ (i ) Set the step length 0.001, change the values of φ0 from Mλθ to min(1, Nλθ ); U (φ )/C η Substituting the values of φ0 into φ1 = e U (0φ )/Cη yields a series of values of φ1 ; 1+e

0

Sift through the values of φ1 which is greater than corresponding φ0 and among which pick the maximum denoted by φ (i ) ; Check if the value of φ (i ) meets the condition of steady-state, if so, output the result, otherwise, delete the case. end

When the service state is visible, λ1 < λ2 , we denote φ1(v ) and φ2(v ) as the joining fractions when the service is unavailable and available, respectively. Utilizing (51) and (52), we solve the system of binary nonlinear equations shown as follows to obtain the joining fractions:

⎧ ⎪ ⎨φ1(v) = ⎪ ⎩φ2(v) =

U (1 ) (φ

e

(v ) )/C η

1

U (1 ) (φ

1+e

U (2 ) (φ

e

(v ) )/C η

,

1

(v ) )/C η

(v ) )/C η . U (2 ) (φ

1+e

2

2

(60)

322

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

Because of the complexity of the function, we still solve it numerically and the concrete steps are shown in Algorithm 2. Then, using the computer software Matlab, we can acquire some numerical results of φ (i) and φ (v ) and performed sensitivity analyses by taking into account the following five cases: Case 1: N = 10, M = 5, λ = 10, θ = 1, β = 40, R0 = 20, p = 10, C = 5, change η from 0 to 10; Case 2: N = 10, M = 5, λ = 10, β = 40, R0 = 20, p = 10, C = 5, η = 1, change θ from 0.2 to 1.5; Case 3: N = 10, M = 5, λ = 10, β = 40, θ = 1, R0 = 20, C = 5, η = 1, change p from 0 to 15; Case 4: N = 10, M = 5, λ = 10, θ = 1, R0 = 20, C = 5, p = 10, η = 1, change β from 1 to 100; Case 5: λ = 30, θ = 1, β = 10, R0 = 20, p = 10, C = 5, η = 1, change M from 1 to 18 for N = 20, 22, 24, 26, 28, respectively. Recall that φ (i) is the fraction of customers who will join the system where the service is invisible. φ1(v ) and φ2(v ) are the fractions of customers who will join the system when service is unavailable and available, respectively. φa(vve) = φ1(v ) P {X ≥ 0} + φ2(v ) P {X < 0} is the average joining fraction when the service is visible. For case 1, whether the service state is visible or invisible, the fraction of customers who will join the system φ (i) , φ1(v) , φ2(v) are all monotonically decreasing in bounded rationality η. As η → 0, means that customers’ joining behav-

ior converges to full rationality, then φ (i ) → 0.762, φ1(v ) → 0.679, φ2(v ) → 1, which are the equilibrium joining probabilities of fully rational customers. And we also conducted some numerical examples to prove that, when η → ∞, φ (i ) , φ1(v ) , φ2(v ) are all converge to 0.5, which exactly implies that the customers tend to be completely irrational.

1

0.9

(v) (0)=1 2

0.95

(i) (v) ave (v) 1 (v) 2

0.9 0.85

0.8 0.7 0.6

0.8

(i)

(0)=0.762

0.5

0.75 0.4

(i)

0.7 (i)

0.65

(v) 1

(v) ave (v) 1 (v) 2

0.3

(v) 2

0.2

0.6

0.1

0.55 (v) (0)=0.679 1 0.5 1

2

3

4

5

6

7

8

9

10

0 0.2

(a) Case 1

0.4

0.6

0.8

1

1.2

1.4 1.5

(b) Case 2

1

0.9 (i) (v) ave (v) 1 (v) 2

0.9

(v) 2

0.85 (i) (v) ave (v) 1 (v) 2

0.8

0.8 0.75 0.7 (i)

0.7 0.6 0.65 0.5

(v) 1

0.6

0.4 0

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15

0.55

p

(c) Case 3

(d) Case 4

Fig. 3. (a)–(g) The the fraction of customers who will join versus various system parameters.

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

0.6 N=20

N=22

N=24

N=26

323

0.865

N=28

N=20 N=22 N=24 N=26 N=28

0.86

0.58

0.855

(v) 2

0.56

(v) 1

0.54

0.85

0.845

0.52 0.84 0.5 0.835 0.48

0.83

0.46 1

2

3

4

5

6

7

8

0.825

9 10 11 12 13 14 15 16 17 18

1

M

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18

M (v)

(v)

(e) Case 5 – φ1

(f) Case 5 – φ2

0.75 N=20 N=22 N=24 N=26 N=28

0.7

(i)

0.65

0.6

0.55

0.5

0.45 1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18

M (g) Case 5 – φ(i) Fig. 3. Continued

Algorithm 2: Calculate the equilibrium joining fractions when service is visible. Input: R0 , p, C, η, λ, θ , β , ε , M, N Output: φ (v ) (= (φ1(v ) , φ2(v ) ))

1

2 3

Set the step length 0.001, change the values of φ1(v ) from 0 to min(1, Nλθ ), change the values of φ2(v ) from

max(φ1(v ) , Mλθ ) to 1; Substituting the values of φ (v ) into the right hand of (60), yields a new series of values, denoted by f (φ (v ) );   Find the value arg min(φ (v ) − f (φ (v ) ) ); ) φ (v



1

(v ) − f (φ (v ) ) ) satisfies the steady condition, if so, output the result, otherwise, delete the case. 4 Check if arg min (φ 1 φ (v ) 5

end

324

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

For case 2, consider the system with invisible service, φ (i) is monotonically increasing in θ , it’s intuitive that the faster the servers arrive, the less time the customers have to wait, so the higher the actual arrival rate. And the monotonicity of φ1(v) can be interpreted in the same way. However, the values of φ2(v) with respect to θ are almost the same, because there must be at least one server in the system when service is available, so the customers’ decisions will be basically unaffected by the change of the arrival rate of servers θ . For case 3, φ1(v ) , φ2(v ) and φ (i) are all monotonically decreasing with respect to p. It’s intuitive the higher the price, the fewer actual arriving customers there will be. For case 4, first of all, we discuss the influence of the servers’ impatience on the sojourn time of customers in both service states. There are two aspects of the influence. I: On the one hand, when β increases, which means servers become more impatient, thus the more likely the server is to leave the system early, this results in a reduction in the second part of sojourn time. II: On the other hand, the increasing of β is equivalent to the decreasing of the service rate. Because servers are supposed to serve N customers, but now, the number of customers served in a batch may be less than N owing to his impatience, which leads to the declining of service rate equivalently, finally results in an increasing of the matching time. Therefore, in figure (d), φ2(v ) is monotonically increasing in β because of the influence I. φ1(v ) increases in β because the influence I is dominant when the service is available for a system with visible service. And φ (i) decrease in β because the influence II is dominant for a system with invisible service in our case. We also find that in this case, as β → ∞, which can be interpreted as that once the number of customers the server has matched with reaches M, they leave the system immediately, φ (i) , φ1(v ) and φ2(v ) all converge to their corresponding limit values shown in figure (d).

For case 5, the fractions of customers who will join the system φ (i) and φ1(v ) are increasing in N for a fixed M. Of course, in practice, when service capacity N increases, the threshold M will also increase. Despite this, according to figure (e) (g), φ (i) (M + i, N + j ) > φ (i) (M, N ), i, j > 0 and φ1(v) (M + i, N + j ) > φ1(v) (M, N ), i, j > 0, so increasing service capacity N has a significant effect on increasing the actual arrival rate of customers when service is invisible or visible and unavailable. And in figure (f), we find φ2(v ) is decreasing in M and N, so φ2(v ) (M + i, N + j ) < φ2(v ) (M, N ), i, j > 0, that is to say, increasing N will lead to a decline of the actual arrival rate of customers when service is visible and available, which is understandable, larger N or M is, the more time new arrival customers have to wait when the service is available. Moreover, we find the curve of φa(vve) almost coincides with the curve of φ (i) in figure (a)-(d), that is to say, whether the service state is visible or invisible, the joining fractions are approximately equal. 7. Optimization analysis of service provider’s utility In this section, we focus on the profit of server’s provider. We assume the servers charge each customer who joins the system an admission fee denoted by p in the above. And the service provider charges part of the profit for each transaction. In addition to the earnings, there are several types of expenditures for the service provider, such as the management fees and some other fixed fees. Combined with the above considerations, we investigate the effects of main parameters θ , p, M and N on the expected utility of service provider per unit time, which is given by

F (θ , p, M, N ) = C0 (P {X ≥ 0}λ1 p + P {X < 0}λ2 p) − C1 E (Ls ) − C2 E (L|X ≥ 0 )P {X ≥ 0} − C3 E (L|X < 0 )P {X < 0} − C4 ,

(61)

where the cost elements appear in provider’s expected utility function (61) are defined as

C0 = the proportion of profit the service provider charges for each transaction; C1 = holding cost per unit time of a server wait in the system; C2 = management fees per unit time should pay for one customer who is waiting for his server; C3 = management fees per unit time should pay for one customer who has been matched and wait to leave; C4 = provider’s fixed cost per unit time. Based on whether the service is visible or not, we discuss service provider’s utility in two cases as well. When service is visible, λ1 = λφ1(v ) (N, M, θ , p) < λ2 = λφ2(v ) (N, M, θ , p). When service is invisible, λ1 = λ2 = λφ (i ) (N, M, θ , p). Then some numerical experiments are conducted to investigate the effects of main parameters (i.e. θ , p, M and N) on the service provider’s utility per unit time. We present three numerical examples with some specific values of these main parameters for both cases as follows:

I: R = 30, C = 5, η = 1, β = 5, λ = 30, M = 14, N = 24, p = 15, and

θ varies in [0.2, 1.6],

II: R = 30, C = 5, η = 1, β = 5, λ = 30, M = 14, N = 24, θ = 1, and p varies in [5, 26.5], III: R = 30, C = 5, η = 1, β = 5, λ = 30, p = 15, θ = 1, and various M for N = 20, 22, 24, 26, 28, 30. Furthermore, set C0 = 0.5, C1 = 20, C2 = 5, C3 = 2, C4 = 1. The numerical results are presented in Fig. 4. For convenience of observation, we translate the negative value of utility into 0. Now, according to the figures of service provider’s utility versus various main parameters, there are some facts can be observed. We articulate the exposition into two parts of invisible service and visible service, separately.

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

250

325

350 invisible service visible service

invisible service visible service

300

service provider's utility

service provider's utility

200

150

100

250 200 150 100

50 50 0 0.2

0.4

0.6

0.8

1

1.2

1.4

0

1.6

5

10

(a) Model 1

p

20

25

(b) Model 2

120

250 N=20 N=22 N=24 N=26 N=28 N=30

100

80

60

40

20

service provider's utility (visible service)

service provider's utility (invisible service)

15

200

N=20 N=22 N=24 N=26 N=28 N=30

150

100

50

0

0 5

10

15

20

25

5

10

(c) Model 3 – invisible

15

20

25

M

M

(d) Model 3 – visible

Fig. 4. The service provider’s utility for both cases versus various main parameters.

We firstly focus on the case of invisible service. For example I and II shown in Fig. 4-(a)–(c), the functions of service provider’s utility with respect to θ , p and M are all unimodal within the given range. For low values of θ or p, the utility is negative, the function starts out monotonically increasing till the peak, then monotonically decreases. For a fixed value of N we set, the value of utility starts out increasing slowly from a positive value till the maximum, then decreases sharply and changes from positive into negative. Moreover, we observe that the maximum points of provider’s utility for N = 20, 22, 24, 26, 28, 30, are located on the auxiliary dotted line given in the Fig. 4-(c) and the maximum utility among the range of M is increasing in N within the given range. Then we consider the case of visible service with the same parameters. Referring to Fig. 4-(a)–(d), the regularity of the utility function with respect to single parameter is not easy to characterize relative to the case of invisible service. In spite of this, we can still find the maximum point within a given parameter range. And compared with the results under the case of invisible service, it not difficult to find that the values of the service provider’s utility for visible service case are basically larger than the values for invisible case. Based on the above analysis and under this utility model we provide, the effect of single parameter on the provider’s utility is clear at a glance. And we provide a reference for the service provider who wants to pursue a maximum profit within a range of single parameter. Moreover, when the parameters are fixed, the service provider may choose to expose their service state to capture greater benefits if permitted. We have discussed the impact of single parameter on service provider’s utility, and obtain some partial optimum results within the given range of these single parameters. However, in some actual systems, there are more than one parameter needs to be considered. To solve the multi-parameter optimization problem of our model, we apply PSO (particle swarm

326

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

110

X: 5 Y: 74.85

76

X: 8 Y: 105.6

105

72

X: 3 Y: 71.51

service provider's utility

service provider's utility

74

X: 4 Y: 71.51

70 68 66 X: 2 Y: 64.91

64

X: 7 Y: 103.3

100

X: 4 Y: 96.41

95

X: 5 Y: 96.41 X: 3 Y: 93.84

90

X: 2 Y: 86.84

85 80

62 60

X: 6 Y: 101.7

X: 1 Y: 75.37

X: 1 Y: 61.31

2

75 4

6

8

10

12

14

2

4

6

iterations

8

10

12

14

iterations

(a) Invisible case

(b) Visible case

X: 21 Y: 1.2 Z: 74.85

80

service provider's utility (visible)

service provider's utility (invisible)

Fig. 5. Iterative convergence of PSO algorithm.

60 40 20 0 2

X: 22 Y: 1.2 Z: 105.6

120 100 80 60 40 20 0 2 1.5

30

1.5

25

1

20 0.5

15

p

1 0.5

10

(a) Invisible case

0

5

10

15

20

25

30

p

(b) Visible case Fig. 6. The service provider’s utility versus various (p, θ ).

optimization) algorithm which is originally attributed to Kennedy and Eberhart [27], Eberhart and Shi [28] and used to solve continuous nonlinear optimization problems. PSO does not use the gradient, so it is quite suitable for our model. For example, the system of passenger departure station, the capacity of each coach has been determined, so the passenger transport companies can only make a reasonable price and a reasonable coach arrival frequency to maximize their interests. Thus, in order to deal with such real problems, we should search for the (θ ∗ , p∗ ) which leads to optimal utility of service provider for the case of fixed N, M. So that the service provider’s problem is formulated as maxθ ,p F (θ , p). We consider the case N = 8, M = 4, λ = 10, R0 = 30, C = 5, η = 1, β = 5, C0 = 0.5, C1 = 10, C2 = 5, C3 = 2, C4 = 1. To cut down the amount of calculation, we discretize the value of the variable θ and p, let θ = 0.1, 0.2, . . . , 2 and p = 0, 1, 2, . . . , 30. In order to verify the correctness of the results obtained by PSO, two three-dimensional graphs which depict the service providers utility versus various (p, θ ) are presented in Fig. 6. And the procedure for searching the optimum solution by PSO is summarized in Algorithm 3. Applying the PSO algorithm, we set the number of particles to 20 and set the maximum iterations to 20, then repeated 10 times and each result is identical, which can be illustrated that PSO is convergent for our model. For invisible case, we present one of the iterative processes shown in Fig. 5-(a). We observe that the iterative process eventually converges on the value of 74.85 and corresponding coordinate is (1.2,21), which is consistent with the results shown in Fig. 6-(a). And for the system with invisible service, the service provider will capture the maximum utility per unit time when the arrival rate of servers is set to 1.2 and the price to 22. Similarly, we obtain the optimal value of provider’s utility for the system with visible service is 105.6, and corresponding coordinate is (1.2,22), which is also consistent with the results shown in

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

327

Algorithm 3: Searching for the optimum solution (θ ∗ , p∗ ) by PSO.

1 2 3 4 5

6

Input: R0 , C, η, λ, β , M, N Output: (θ ∗ , p∗ ) Set the number of particles n, range of (θ , p) and maximum iterations. Generate n particles randomly, denoted as {X1 , X2 , . . . , Xn },where Xi = (θi , pi ). while iterations < maximum iterations do Calculate the fitness value (i.e. the utility value) of each particle. Update the coordinates of the local optimum {P best1, P best2, . . . P bestn} andthe coordinates of the global optimum Gbest. Update every particle coordinates by the formulas:

vki = vki −1 + 2 ∗ rand ∗ (Pbestik − Xik−1 ) + 2 ∗ rand ∗ (Gbestik − Xik−1 ), Xik = Xik−1 + vki , 7

where k represents the kth iteration, i represents the ith particle, and rand is a random variable in [0,1]. end

Fig. 6-(b), which implies that the provider should set the arrival rate of servers to 1.2 and the price to 22 to maximize their profits. And referring to Fig. 6, we can observe that PSO applied in our model requires few iterations but leads to an accurate result, which exactly imply the superiority of the PSO algorithm. Furthermore, for some other actual systems, the service provider may consider other combinations of the parameters, and may encounter higher dimensional optimization problems. For these various optimization problems with multidimensional parameters, employing the PSO algorithm is still a well choice which will reduce quite a lot amount of computation. And the specific algorithm flows are the same as Algorithm 3 except the target parameters, so they won’t be covered here. 8. Conclusions In this paper, we have researched on a batch matching system with impatient servers and boundedly rational customers. To formulate the system and solve the steady-state probabilities of system sates, we combine the generating function method and the G-matrix method. Consequently, we obtain the expressions for some significant performance measures in terms of the effective arrival rate of customers. Taking the bounded rationality of customers into account, we numerically capture the equilibrium joining fraction of customers. Further some sensitivity analyses are conducted to explore the impact of target parameters on the equilibrium joining fraction of customers, and reveal the local monotonicity with regard to each target parameter. In our examples, the results show that the joining fraction of customers is less than equilibrium joining probability of fully rational customers. Finally we construct a utility function of the service provider, and employ the PSO algorithm to search the maximum point of provider’s utility, which can bring some management implications to the service provider. It’s found that the server’s impatience doesn’t necessarily lead to a loss but sometimes increases benefits. Therefore, the provider can utilize the mechanism of servers’ impatience to acquire greater benefits. Future work might extend this model, for instance, it would make sense to consider the case of limited capacity of customers area or servers area. And the fully observable case is remained to be discussed in future research. Acknowledgements This work is supported by the National Natural Science Foundation of China (No. 61773014), and Postgraduate Research & Practice Innovation Program of Jiangsu Province, China (No. KYCX18_0373). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]

D.G. Kendall, Some problems in the theory of queues, J. R. Stat. Soc. 13 (2) (1951) 151–185. J.M. Dobbie, Letter to the editor doubled-ended queuing problem of Kendall, Oper. Res. 9 (5) (1961) 755–757. B.R.K. Kashyap, The double-ended queue with bulk service and limited waiting space, Oper. Res. 14 (5) (1966) 822–834. B.R.K. Kashyap, Further results for the double ended queue, Metrika 11 (1) (1967) 168–186. A.D. Crescenzo, V. Giorno, B.K. Kumar, A.G. Nobile, A double-ended queue with catastrophes and repairs, and a jump-diffusion approximation, Methodol. Comput. Appl. Probab. 14 (4) (2012) 937–954. Y. Shi, Z. Lian, Optimization and strategic behavior in a passenger-taxi service system, Eur. J.Oper. Res. 249 (3) (2016) 1024–1032. F. Wang, J. Wang, Z.G. Zhang, Strategic behavior and social optimization in a double-ended queue with gated policy, Comput. Ind. Eng. 114 (2017) 264–273. B. Prabhakar, N. Bambos, T.S. Mountford, The synchronization of poisson processes and queueing networks with service and synchronization nodes, Adv. Appl. Probab. 32 (3) (20 0 0) 824–843. M.F. Neuts, A general class of bulk queues with poisson input, Ann. Math. Stat. 38 (3) (1967) 759–770. M.L. Chaudhry, J.G.C. Templeton, A first course in bulk queues, J. Oper. Res. Soc. 34 (12) (1983) 1222–1223. T. Jiang, L. Liu, Y. Zhu, Analysis of a batch service polling system in a multi-phase random environment, Methodol. Comput. Appl. Probab. 20 (2) (2018) 699–718.

328

X. Chai, L. Liu and B. Chang et al. / Applied Mathematics and Computation 354 (2019) 308–328

[12] T. Jiang, L. Liu, Analysis of a batch service multi-server polling system with dynamic service control, J. Ind. Manag. Optim. 14 (2) (2018) 743–757. [13] O. Bountali, A. Economou, Equilibrium joining strategies in batch service queueing systems, Eur. J. Oper. Res. 260 (2017) 1142–1151. [14] Q. Bu, L. Liu, T. Jiang, A limited clearing queueing model with an orbit and non-persistent customers, Qual. Technol. Quant. Manag. (2017), doi:10. 1080/16843703.2017.1396955. [15] D.Y. Barrer, Queuing with impatient customers and indifferent clerks, Oper. Res. 5 (5) (1957) 644–649. [16] E. Altman, U. Yechiali, Analysis of Customers’ Impatience in Queues with Server Vacations, Queueing Syst. 52 (4) (2006) 261–279. [17] A. Brandt, M. Brandt, On the M(n)/M(n)/s queue with impatient calls, Perform. Eval. 35 (1-2) (1999) 1–18. [18] P. Naor, The regulation of queue size by levying tolls, Econometrica 37 (1) (1969) 15–24. [19] U. Yechiali, On Optimal Balking Rules and Toll Charges in the GI/M/1 Queuing Process, Oper. Res. 19 (2) (1971) 349–370. [20] X. Su, Bounded rationality in newsvendor models, Manuf. Serv. Oper. Manag. 10 (4) (2008) 566–589. [21] T. Huang, G. Allon, A. Bassamboo, Bounded Rationality in Service Systems, Manuf. Serv. Oper. Manag. 15 (2) (2013) 263–279. [22] M.F. Neuts, Matrix-Geometric Solutions to Stochastic Models, Springer Berlin Heidelberg, 1984. [23] G. Latouche, V. Ramaswami, Introduction to Matrix Analytic Methods in Stochastic Modeling, ASA-SIAM Series On Statistics and Applied Probability, SIAM, Philadelphia, PA, 1999. [24] Neuts, Structural stochastic matrices of M/G/1 type and their applications, Acta Appl. Math. 24 (2) (1991) 205–206. [25] T. Burr, Introduction to Matrix Analytic methods in Stochastic Modeling, Technometrics 43 (1999) 379–380. [26] J.D. Hey, C. Orme, Investigating generalizations of expected utility theory using experimental data, Econometrica 62 (6) (1994) 1291–1326. [27] J. Kennedy, R. Eberhart, Particle Swarm Optimization, Springer, US, 2011. [28] Y. Shi, R. Eberhart, Modified particle swarm optimizer, in: Proceedings of the IEEE ICEC Conference, 1998, pp. 69–73. Anchorage.