Fluid approximation analysis of a call center model with time-varying arrivals and after-call work

Kawai, Yosuke; Takagi, Hideaki Article Fluid approximation analysis of a call center model with time-varying arrivals and after-call work Operations Research Perspectives Provided in Cooperation with: Elsevier Suggested Citation: Kawai, Yosuke; Takagi, Hideaki (2015) : Fluid approximation analysis of a call center model with time-varying arrivals and after-call work, Operations Research Perspectives, ISSN 2214-7160, Elsevier, Amsterdam, Vol. 2, pp. 81-96, https://doi.org/10.1016/j.orp.2015.03.003 This Version is available at: https://hdl.handle.net/10419/178252 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. http://creativecommons.org/licenses/by-nc-nd/4.0/ Operations Research Perspectives 2 (2015) 81–96 Contents lists available at ScienceDirect Operations Research Perspectives journal homepage: www.elsevier.com/locate/orp Fluid approximation analysis of a call center model with time-varying arrivals and after-call work Yosuke Kawaia, Hideaki Takagib,∗ aGraduate School of Systems, Information and Engineering, University of Tsukuba, Tsukuba Science City, Ibaraki 305-8573, Japan bFaculty of Engineering, Information and Systems, University of Tsukuba, Tsukuba Science City, Ibaraki 305-8573, Japan article info Article history: Available online 4 April 2015 Keywords: Fluid approximation Multiserver queue Call center Time-varying arrival process After-call work Abandonment abstract Important features to be included in queueing-theoretic models of the call center operation are multiple servers, impatient customers, time-varying arrival process, and operator’s after-call work (ACW). We propose a fluid approximation technique for the queueing model with these features by extending the analysis of a similar model without ACW recently developed by Liu and Whitt (2012). Our model assumes that the service for each quantum of fluid consists of a sequence of two stages, the first stage for the conversation with a customer and the second stage for the ACW. When the duration of each stage has exponential, hyperexponential or hypo-exponential distribution, we derive the time-dependent behavior of the content of fluid in each stage of service as well as that in the waiting room. Numerical examples are shown to illustrate the system performance for the cases in which the input rate and/or the number of servers vary in sinusoidal fashion as well as in adaptive ways and in stationary cases. ©2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction Queueing models have been widely used to model the performance of call centers with impatient customers [1–4], which means that customers in the waiting line may leave before getting service. The multiserver queue M/M/swith impatient customers is called Erlang-A model, ‘‘A’’ for ‘‘abandonment’’, in contrast with the well-known Erlang-B model (M/M/s/s) and the Erlang-C model (M/M/swith only patient customers). Through the measurements at real call centers, however, we observe that operators usually spend sizable amount of time to complete additional work after finishing conversation with customers. For example, they enter customer profiles and summary of conversation into the customer management database after conversation. Such extra work of operators is called the after-call work (ACW). Cleveland and Harne [5, Section 8] describe: The ACW is the work that is necessitated by and immediately follows an inbound transaction. Often includes entering data, filling out forms and making outbound calls necessary to complete the transaction. The agent is unavailable to receive another inbound call while in this mode. ∗Corresponding author. Tel.: +81 29 853 5414; fax: +81 29 853 7291. E-mail address: [email protected] (H. Takagi). The ACW is also called ‘‘post call activity’’ [6–8], ‘‘wrap-up times’’ [1], ‘‘after-hung-up times’’ [9], and ‘‘postservice activity’’ [10,11]. Harris and Phillips [6] mention: The post call activity is a phase in which the operator may fill out dockets, make supplementary phone calls or perform other clerical activities before pressing a key to indicate that he/she is able to accept another call from the queue (if such a call is present). Takagi and Taguchi [12] study a two-dimensional birth-anddeath process for the M/M/K/Jqueue with ACW, where K, the number of servers, represents the total number of operators working in the call center and J, the maximum number of customers accommodated in the system, stands for the number of incoming telephone lines. Unlike usual queueing models, we do not necessarily assume that J≥K, because servers may be working on ACW while some customers are present in the waiting room. Phung-Duc and Kawanishi [13] present a matrix-geometric analysis for a queueing model with retrial arrivals of blocked and abandoned customers. All models in these pieces of work assume the steady state of the system. Another realistic feature of call center operation is that the call input process is time-varying. However, the exact stochastic analysis of a queueing model with multiple servers with generally distributed service times and/or time-varying arrival process is not easy. The fluid approximation technique has been exploited to http://dx.doi.org/10.1016/j.orp.2015.03.003 2214-7160/©2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4. 0/). 82 Y. Kawai, H. Takagi / Operations Research Perspectives 2 (2015) 81–96 Fig. 1. System model and state variables in the underloaded state. Fig. 2. System model and state variables in the overloaded state. deal with such models traditionally [14]. More recently, the fluid approximation is applied to stationary multiserver queues with impatient customers [15,16] as well as those with the time-varying input rate and number of servers [17–19]. In this paper, we present the fluid approximation for a multiserver queueing model with impatient customers, two stages of service time (representing the conversation and ACW in a call center), and time-varying input rate and number of servers. Our approach is an extension of the method for the Mt/GI/st+GI model originally developed by Liu and Whitt [20,19] (they later extended the analysis to networks of fluid queues [21,22]). In this notation, ‘‘Mt’’ means a Poisson arrival process with time-varying arrival rate, the first ‘‘GI’’ an independent, generally distributed service time, ‘‘st’’ a time-varying number of servers, and ‘‘+GI’’ a general abandonment-time distribution. We show that the system state alternates between the underloaded interval in which there are idle servers and the overloaded interval in which arriving fluid quanta must wait for service because all servers are busy. We study the dynamics of the fluid content in service in both underloaded and overloaded states. We also study the dynamics of the fluid content in the waiting room and the waiting time of a fluid quantum that arrives in the overloaded state. Our analysis is applied to several illustrative cases with time-varying input rate and number of servers. If the number of servers is determined adaptively to cope with only the load of conversation, the system is always overloaded but it remains stable. If the number of servers is determined adaptively in accordance with the load of both conversation and ACW, the system is always underloaded. To the best of the authors’ knowledge, this paper is the first work in which the fluid approximation is applied to the Mt/GI/st+ GI model with two stages of service time as a model of the call center operation with ACW. This paper is partly based on the Master Thesis of the first author [23] submitted to the Graduate School of Systems, Information and Engineering of the University of Tsukuba, Japan. 2. Fluid model of call center operation with after-call work In this section, we introduce a fluid model of the call center operation with ACW by extending the model and analysis by Liu and Whitt [19,16]. 2.1. Definition of the system model We consider a fluid queueing system with multiple servers where incoming calls in a call center are modeled by quanta of fluid. We assume that the service time a server, representing an operator, spends on each fluid quantum consists of a sequence of two stages, called ‘‘service 1’’ for the conversation with a customer and ‘‘service 2’’ for ACW, each having independent duration. We assume that the same server continues to provide service 2 immediately after service 1 for each fluid quantum. Let there be s(t)servers in the system at time t≥0. The staffing function s(t)is given exogenously or adaptively somehow depending on the input rate of fluid. At any time, each server is either in service or being idle such that si(t)servers are engaged in service i(i=1,2), where s(t)≥s1(t)+s2(t). The input of fluid quantum directly enters service 1 if there is a server available; this state is called underloaded. Otherwise, the input flows into the ‘‘waiting room’’ for service 1; this state is called overloaded. No waiting room is needed for service 2 because the same server takes care of the ACW for the fluid quantum that he has just given service 1. The server who has finished service 2 can start service 1 for another fluid quantum if any in the waiting room, or he becomes idle otherwise. The fluid quanta leave the system either by completing service 2 or by abandonment while being in the waiting room. The fluid quanta never leave the system during services 1 and 2. For a system with time-varying staffing function, we assume that the time variation in the total number of servers is solely turned to the time variation in the number of servers assigned to service 1. We also assume that the number of servers assigned to each service never goes below the level of fluid content in that service at any moment so that no fluid quanta are forced out of the system once they have entered service. This model is schematically depicted along with relevant state variables in Figs. 1 and 2for the underloaded and overloaded states, respectively. The state variables are introduced in the following subsection. 2.2. Definition of state variables and their relations We assume that the fluid quanta arrive at service 1 according to a deterministic process with time-varying rate λ(t), t≥0. We denote by F(x)and f(x)the distribution function and the probability density function (pdf), respectively, for the abandonment time of each fluid quantum in the waiting room. Also, we denote by Gi(x) and gi(x)the distribution and density functions, respectively, for the service time of each fluid quantum in service i(i=1,2). Thus we have F(x):= x 0 f(u)du,Gi(x):= u 0 gi(u)du x ≥0,i=1,2. Furthermore, let F(x)and Gi(x)be their complimentary distribution functions (CDF’s) defined by F(x):= 1−F(x), Gi(x):= 1−Gi(x)x≥0,i=1,2. These functions are assumed to be given in the model. At time t(≥0), we denote by Q(t,x)the fluid content that has been waiting for the time units less than or equal to xin the waiting room. Similarly, we denote by Bi(t,x)the fluid content in service that has been in service ifor the time units less than or Y. Kawai, H. Takagi / Operations Research Perspectives 2 (2015) 81–96 83 Fig. 3. Alternation of underloaded and overloaded intervals. equal to x(i=1,2). Furthermore, let q(t,x)and bi(t,x)be their corresponding density functions as Q(t,x)=x 0 q(t,u)du,Bi(t,x)=x 0 bi(t,u)du x≥0,i=1,2.(1) Then the fluid content in the waiting room and that in service iat time tare respectively given by Q(t):= Q(t,∞), Bi(t):= Bi(t,∞)i=1,2.(2) Let us define w(t)as the largest waiting time for the fluid in the waiting room: w(t):= inf {y≥0:q(t,x)=0 for all x>y}t≥0. If we assume that fluid quanta enter service in the order of arrival, i.e., according to the first-come first-served (FCFS) discipline, w(t) is the waiting time of a fluid quantum that enters service at time t. Since q(t,x)=0x> w(t), w(t)is called the boundary waiting time by Liu and Whitt [19]. The above-defined density functions for the fluid content satisfy the following fundamental evolution equations [19, Assumption 6]: q(t+y,x+y)=q(t,x)F(x+y) F(x)0≤x< w(t)−y,(3) bi(t+y,x+y)=bi(t,x)Gi(x+y) Gi(x)x≥0,i=1,2.(4) Eq. (3) simply says that the fluid quanta present in the waiting room at time tthat have not abandoned for xtime units and do not do so for ymore time units remain there at time t+y. Similarly, Eq. (4) says that the fluid quanta in service iat time tthat have not completed service for xtime units and do not do so for ymore time units remain in service iat time t+y(i=1,2). We define the instantaneous ending rates, i.e., the hazard-rate functions, for the abandonment time and the service times by hF(x):= f(x) F(x),hGi(x):= gi(x) Gi(x)x≥0,i=1,2. Then the rate α(t)of fluid that abandons at time tand the output rate σi(t)of fluid that service iis completed at time tare respectively given by α(t)=∞ 0 q(t,x)hF(x)dx =w(t) 0 q(t,x)hF(x)dx 0≤w(t)≤t,(5) σi(t)=∞ 0 bi(t,x)hGi(x)dx t ≥0,i=1,2.(6) 2.3. Alternation of underloaded and overloaded intervals We assume that the system state, started at time t0, alternates between the underloaded and overloaded intervals. Let us denote by {t2n;n=0,1,2, . . .}a sequence of epochs at which underloaded intervals are started, and denote by {t2n+1;n=0,1,2, . . .} a sequence of epochs at which overloaded intervals are started such that [t2n,t2n+1] : underloaded interval; [t2n+1,t2n+2] : overloaded interval n=0,1,2,.... See Fig. 3 for the alternation of underloaded and overloaded intervals. The system is said to be in the underloaded state if a fluid quantum that arrives enters service 1 immediately because there are more servers than the total fluid content in the system. Therefore, if the system is underloaded at time t∈ [t2n,t2n+1], we have Q(t)=0,s(t) > B1(t)+B2(t)t∈ [t2n,t2n+1].(7) The underloaded interval [t2n,t2n+1]ends when the total fluid content in services 1 and 2 becomes equal to the total number of servers for the first time after t2n. Thus the termination epoch t2n+1 of the underloaded interval is determined by the condition t2n+1=inf{t≥t2n:s(t)=B1(t)+B2(t)}n=0,1,2,.... (8) The system is said to be in the overloaded state if a fluid quantum that arrives cannot enter service 1 because there are no servers available. Therefore, if the system is overloaded at time t∈ [t2n+1,t2n+2], we have Q(t) > 0,s(t)=B1(t)+B2(t)t∈ [t2n+1,t2n+2].(9) The overloaded interval [t2n+1,t2n+2]ends when the fluid content in the waiting room vanishes. Thus the termination epoch t2n+2of the overloaded interval is determined by the condition t2n+2=inf{t≥t2n+1:Q(t)=0}n=0,1,2,.... (10) 2.4. Fluid in service in the underloaded state Let us first study the fluid content in service i(i=1,2)at time twhen the system is in the underloaded interval. For the simplicity of notation, we assume in this subsection that the underloaded interval of our concern is started at time 0 without loss of generality. For service 1, from Eq. (4) with i=1, we can derive the following transport partial differential equation for b1(t,x)(Liu and Whitt [19, online version, Appendix B]): ∂b1(t,x) ∂t+∂b1(t,x) ∂x= −hG1(x)b1(t,x)(11) with initial condition b1(0,x)and the boundary condition b1(t,0)=λ(t)t≥0,(12) which is the rate of fluid going into service 1. Then the solution is given by (Liu and Whitt [19, Proposition 2]): b1(t,x)=G1(x)λ(t−x)1{x≤t}+G1(x) G1(x−t)b1(0,x−t)1{x>t},(13) where 1Eis the indicator function for event Edefined by 1E:= 1 if Eis true, 0 otherwise. See Fig. 4 for the two cases x≤tand x>tof contribution to b1(t,x). If x≤t, the new arrival of fluid quanta with rate λ(t−x)at time t−xthat stay longer than xtime units contributes to b1(t,x). If 84 Y. Kawai, H. Takagi / Operations Research Perspectives 2 (2015) 81–96 (a) x≤t. (b) x>t. Fig. 4. Two cases for calculating b1(t,x)at time tin the underloaded state. x>t, the fluid content b1(0,x−t)at time 0 that stays longer than xtime units contributes to b1(t,x). From Eq. (13) we can derive B1(t)=t 0 G1(x)λ(t−x)dx +∞ 0 G1(x+t) G1(x)b1(0,x)dx.(14) From Eqs. (6) and (13), the output rate from service 1 is given by σ1(t)=t 0 g1(x)λ(t−x)dx +∞ 0 g1(x+t) G1(x)b1(0,x)dx. For service 2, by similar manipulation of Eq. (4) with i=2, we get the partial differential equation for b2(t,x): ∂b2(t,x) ∂t+∂b2(t,x) ∂x= −hG2(x)b2(t,x)(15) with initial condition b2(0,x)and the boundary condition b2(t,0)=σ1(t)t≥0.(16) Then we can derive b2(t,x)=G2(x)σ1(t−x)1{x≤t}+G2(x) G2(x−t)b2(0,x−t)1{x>t}, B2(t)=t 0 G2(x)σ1(t−x)dx +∞ 0 G2(x+t) G2(x)b2(0,x)dx, σ2(t)=t 0 g2(x)σ1(t−x)dx +∞ 0 g2(x+t) G2(x)b2(0,x)dx.(17) During the underloaded interval, there are s(t)−B1(t)−B2(t) > 0 idle servers. 2.5. Fluid in service and in the waiting room in the overloaded state We next consider the fluid content present in the system at time twhen it is in the overloaded interval. Again we assume in this subsection that the overloaded interval of our concern is started at time 0 without loss of generality. From Eq. (4), the fluid content in service i(i=1,2)is governed by b1(t,x)=G1(x)b1(t−x,0)1{x≤t} +G1(x) G1(x−t)b1(0,x−t)1{x>t},(18) b2(t,x)=G2(x)b2(t−x,0)1{x≤t} +G2(x) G2(x−t)b2(0,x−t)1{x>t}.(19) In these equations, while b1(0,x)and b2(0,x)are given as initial conditions, we must know the boundary conditions b1(t,0)and b2(t,0). This can be done in principle by extending the method of Liu and Whitt [19, Section 6.2] for the solution to the fixed-point equation for our Mt/GI/st+GI model as follows. Since the input rate to service 2 equals the output rate of service 1, we have b2(t,0)=σ1(t)=t 0 g1(x)b1(t−x,0)dx +∞ 0 g1(x+t) G1(x)b1(0,x)dx.(20) In the overloaded state, the fluid entrance rate into service 1 is the sum of the rate at which the number of servers is increased and the output rate of service 2. Thus we have b1(t,0)=s′(t)+σ2(t)=s′(t)+t 0 g2(x)b2(t−x,0)dx +∞ 0 g2(x+t) G2(x)b2(0,x)dx.(21) Eqs. (20) and (21) provide the set of simultaneous integral equations for the functions b1(t,0)and b2(t,0)given s′(t)and the initial conditions b1(0,x)and b2(0,x). During the overloaded interval, all servers are busy, each being engaged either in service 1 or in service 2, so that the fluid content in service iequals the number of servers for service i: Bi(t)=si(t)i=1,2;s(t)=s1(t)+s2(t). (22) There are no idle servers in the overloaded interval. Let us derive the equations governing the number si(t)of servers in service i(i=1,2). The fluid quantum whose service 1 is finished at time tenters service 2 immediately with rate σ1(t). At this moment, the number of servers in service 1 decreases with rate σ1(t)and the number of servers in service 2 increases with the same rate. The fluid quantum whose service 2 is finished at time t leaves the system with rate σ2(t). The server then turns to service 1 immediately. At this moment, the number of servers in service 2 decreases with rate σ2(t)and the number of servers in service 1 increases with the same rate. Since the time variation in the total number of servers is assumed to occur in service 1, we have the following set of ordinary differential equations for {s1(t), s2(t)}: ds1(t) dt = −σ1(t)+σ2(t)+s′(t), ds2(t) dt =σ1(t)−σ2(t). (23) Recall our assumption that the number of servers for service i never goes below the level of fluid content in service i(i=1,2). These conditions are expressed as s1(t)+s′(t)∆t+σ2(t)∆t−σ1(t)∆t≥B1(t)−σ1(t)∆t and s2(t)+σ1(t)∆t−σ2(t)∆t≥B2(t)−σ2(t)∆t, where the left-hand side of each equation shows the number of servers while the right-hand side shows the fluid content at time t+∆t. The second inequality is always satisfied as σ1(t)≥0. From the first inequality, a sufficient condition for the feasibility is given by [19, online version, Appendix G.2] b1(t,0)=s′(t)+σ2(t)≥0.(24) Therefore, we impose the condition in Eq. (24) for the set of equations in Eq. (23). We next consider the fluid content in the waiting room during the overloaded interval. From Eq. (3), we can derive the partial differential equation for q(t,x): ∂q(t,x) ∂t+∂q(t,x) ∂x= −hF(x)q(t,x)0≤x≤w(t)(25) Y. Kawai, H. Takagi / Operations Research Perspectives 2 (2015) 81–96 85 with initial condition q(0,x)and the boundary condition q(t,0)=λ(t)t≥0.(26) Then the density function of the fluid content in the waiting room is given by q(t,x)=F(x)λ(t−x)1{x≤min{w(t),t}} +F(x) F(x−t)q(0,x−t)1{t<x≤w(t)}.(27) This expression has the following meaning; see Fig. 4 for a similar situation in the underloaded state. For x≤t,q(t,x)dx =F(x)λ(t− x)dx is the quantity of fluid quanta that arrive during [t−x−dx,t− x]and do not abandon for xtime units. For x>t,q(t,x)dx = [F(x)/F(x−t)]q(0,x−t)dx is the fluid content in the waiting room at time twhich comes from the fluid content q(0,x−t)dx at time 0 that does not abandon for xtime units. From Eq. (27) we can derive Q(t)=min{w(t),t} 0 F(x)λ(t−x)dx +1{w(t)>t}w(t)−t 0 F(x+t) F(x)q(0,x)dx. 2.6. Waiting time of a fluid quantum in the overloaded state Let us now derive the differential equation for the boundary waiting time w(t)by following Liu and Whitt [19, online version, Appendix D.1]. To do so, for 0 ≤w(t)≤t, we note Q(t)=w(t) 0 λ(t−x)F(x)dx =t t−w(t) λ(x)F(t−x)dx.(28) Differentiating the rightmost side of this equation with respect to t, we get dQ (t) dt =λ(t)F(0)−λ[t−w(t)]F[w(t)][1−w′(t)] −t t−w(t) λ(x)f(t−x)dx =λ(t)−q[t, w(t)][1−w′(t)] − α(t), (29) where, from Eq. (5) we have used α(t)=w(t) 0 q(t,x)hF(x)dx =w(t) 0 λ(t−x)f(x)dx =t t−w(t) λ(x)f(t−x)dx.(30) On the other hand, the fluid content in the waiting room varies as dQ (t) dt =λ(t)−α(t)−b1(t,0). (31) Comparing Eqs. (29) and (31), we obtain the relation b1(t,0)=q[t, w(t)][1−w′(t)], which leads to the first-order nonlinear differential equation for w(t): dw(t) dt =1−b1(t,0) λ[t−w(t)]F[w(t)]t≥0 (32) with initial condition w(0)=0, where b1(t,0)is given as the solution to the set of Eqs. (20) and (21). In addition, Liu and Whitt [19, Section 7.3] consider the potential waiting time v(t). This is defined as the virtual waiting time of an arriving fluid quantum at time twhich elects never to abandon. Since the waiting time of the fluid quantum that is entering service 1 at time tis w(t), then this quantum must have entered the waiting room w(t)time units ago. This implies that the potential waiting time at t−w(t)is w(t). It follows that v[t−w(t)] = w(t)or v(t)=w[t+v(t)].(33) Then the differential equation for v(t)is given by dv(t) dt =λ(t)F[v(t)] b1[t+v(t), 0]−1t≥0 (34) with initial condition v(0)=0. 2.7. Stationary model We show the results for the stationary fluid model of the current model, which is again an extension of the analysis by Whitt [16, Section 3]. In this case we assume that λ(t)≡λ, si(t)≡sit≥0,i=1,2 and that the distribution functions F(x),G1(x), and G2(x)are given. We first consider the case in which both services 1 and 2 are underloaded. In this case, for service 1 we have the density of fluid content given by b1(x)=λG1(x)x≥0, in particular, b1(0)=λ. Thus we get the total fluid content in service 1: B1=∞ 0 b1(x)dx =λ∞ 0 G1(x)dx =λE[G1]. Then the rate of fluid at which service 1 is completed equals the input rate: σ1=∞ 0 b1(x)hG1(x)dx =λ∞ 0 g1(x)dx =λ. Similarly, for service 2 we have the density of fluid content given by b2(x)=σ1G2(x)x≥0, in particular, b2(0)=σ1=λ. Thus we get the total fluid content in service 2: B2=∞ 0 b2(x)dx =σ1∞ 0 G2(x)dx =σ1E[G2] = λE[G2]. Then the rate of fluid at which service 2 is completed also equals the input rate: σ2=∞ 0 b2(x)hG2(x)dx =σ1∞ 0 g2(x)dx =σ1=λ. This case occurs if Bi<si(i=1,2), or λ < min s1 E[G1],s2 E[G2]. We next consider the overloaded case for service 1. Due to the abandonment of fluid in the waiting room, there is certainly a stationary state in the overloaded case for service 1. There is no stationary state in the overloaded case for service 2. In this case, the density of the fluid content in the waiting room is given by q(x)=λF(x)1{x≤w}x≥0, 86 Y. Kawai, H. Takagi / Operations Research Perspectives 2 (2015) 81–96 in particular, q(0)=λ. We determine wlater by Eq. (37). Then the total fluid content in the waiting room is given by Q=∞ 0 q(x)dx =λw 0 F(x)dx. The rate of fluid that abandons in the waiting room is given by α=∞ 0 q(x)hF(x)dx =λw 0 f(x)dx =λF(w). (35) Considering the stationary case in Eq. (31), we get the balance relation λ=α+b1(0). (36) We can find b1(0)as follows. Since b1(x)=b1(0)G1(x)for x≥0, we get B1=b1(0)∞ 0 G1(x)dx =b1(0)E[G1]. It follows from the overloaded condition B1=s1that b1(0)=s1 E[G1]. Substituting this result into Eqs. (35) and (36), we obtain the equation to solve for w: F(w) =s1 λE[G1].(37) The rate of fluid that service 1 is completed is given by σ1=∞ 0 b1(x)hG1(x)dx =b1(0)∞ 0 g1(x)dx =b1(0)=s1 E[G1]. Then the rate of abandonment is given by α=λF(w) =λ1−F(w)=λ−s1 E[G1]. Note that service 2 must be underloaded because there is no waiting room for service 2. In this case, we again have b2(x)= σ1G2(x)for x≥0, which leads to σ2=σ1=b1(0);B2=σ1E[G2] = b1(0)E[G2] = s1E[G2] E[G1]. This case occurs if α > 0, or λ > s1 E[G1]. 3. System with exponentially distributed conversation and ACW times We study a special system with exponentially distributed conversation and ACW times for more explicit analysis. For such a system, we assume Gi(x)=1−e−µixx≥0,i=1,2,(38) which means that the associated hazard-rate functions are constant: hGi(x)=µi(i=1,2). Then the mean conversation time is 1/µ1and the mean ACW time is 1/µ2. The solution in this case is straightforward as we show below for the underloaded interval [t2n,t2n+1]and the overloaded interval [t2n+1,t2n+2],n= 0,1,2, . . . , defined in Section 2.3. 3.1. Solution for the underloaded state Let us consider the underloaded state at time t∈ [t2n,t2n+1]. From Eq. (13), we have b1(t,x) =   G1(x)λ(t−x)0≤x≤t−t2n, G1(x) G1(x−t+t2n)b1(t2n,x−t+t2n)x>t−t2n =e−µ1xλ(t−x)0≤x≤t−t2n, e−µ1(t−t2n)b1(t2n,x−t+t2n)x>t−t2n. In particular, we can confirm Eq. (12). Therefore, or from Eq. (14), we obtain the fluid content B1(t)in service 1 at time tas B1(t)=t−t2n 0 G1(x)λ(t−x)dx +∞ 0 G1(x+t−t2n) G1(x)b1(t2n,x)dx =t−t2n 0 λ(t−x)e−µ1xdx +e−µ1(t−t2n)B1(t2n). (39) The output rate of service 1 is simply given by σ1(t)=∞ 0 b1(t,x)hG1(x)dx =µ1∞ 0 b1(t,x)dx =µ1B1(t). (40) Similarly, from Eq. (17), we have b2(t,x) =   G2(x)σ1(t−x)0≤x≤t−t2n, G2(x) G2(x−t+t2n)b2(t2n,x−t+t2n)x>t−t2n =µ1e−µ2xB1(t−x)0≤x≤t−t2n, e−µ2(t−t2n)b2(t2n,x−t+t2n)x>t−t2n. Thus we obtain the fluid content B2(t)in service 2 at time tas B2(t)=t−t2n 0 σ1(t−x)e−µ2xdx +e−µ2(t−t2n)∞ t−t2n b2(t2n,x−t+t2n)dx =µ1t−t2n 0 B1(t−x)e−µ2xdx +e−µ2(t−t2n)B2(t2n). (41) The output rate of service 2 is given by σ2(t)=µ2B2(t). (42) We note a great merit of exponentially distributed conversation and ACW times that we can obtain Bi(t)in Eqs. (39) and (41) and then σi(t)in Eqs. (40) and (42) without finding bi(t,x),i=1,2. The initial values Bi(t2n), i=1,2, for the underloaded interval [t2n,t2n+1]are given from the solution for the preceding overloaded interval [t2n−1,t2n]as B1(t2n)=s(t2n)−B2(t2n), B2(t2n)=e−(µ1+µ2)(t2n−t2n−1)B2(t2n−1) +µ1e−(µ1+µ2)t2nt2n t2n−1 e(µ1+µ2)us(u)du from Eq. (47) in the sequel. The termination epoch t2n+1of the underloaded interval [t2n,t2n+1]is found by the condition in Eq. (8). Y. Kawai, H. Takagi / Operations Research Perspectives 2 (2015) 81–96 87 3.2. Solution for the overloaded state Let us consider the overloaded state at time t∈ [t2n+1,t2n+2]. The output rate of fluid at which service iis completed at time tis obtained from Eqs. (6) and (22) as σi(t)=∞ 0 bi(t,x)hGi(x)dx =µi∞ 0 bi(t,x)dx =µiBi(t)=µisi(t)i=1,2.(43) Substituting Eq. (43) into Eq. (23), we have the following set of simultaneous linear differential equations for {B1(t), B2(t)}: dB1(t) dt = −µ1B1(t)+µ2B2(t)+s′(t), (44) dB2(t) dt =µ1B1(t)−µ2B2(t). (45) Using s(t)=B1(t)+B2(t), we get the first-order differential equation for B2(t): dB2(t) dt =µ1s(t)−(µ1+µ2)B2(t)t2n+1≤t≤t2n+2,(46) whose solution is given by B2(t)=e−(µ1+µ2)(t−t2n+1)B2(t2n+1) +µ1e−(µ1+µ2)tt t2n+1 e(µ1+µ2)us(u)du,(47) with the initial condition B2(t2n+1)from the fluid content B2(t)at the end of the preceding underloaded interval [t2n,t2n+1]. We then get B1(t)=s(t)−B2(t). In the differential equation (32) for w(t), we have noted that b1(t,0)should be given as the solution to the fixed-point equations (20) and (21) in general. However, it can be obtained easily in the present case in which both conversation and ACW times are exponentially distributed. Since service 1 is started by the servers who have completed service 2 as well as by the servers who are added to the system, we have the relation b1(t,0)=µ2B2(t)+s′(t), (48) which is assumed to be nonnegative from the feasibility condition in Eq. (24). Substituting Eq. (48) into Eq. (32), we get the differential equation for w(t): dw(t) dt =1−µ2B2(t)+s′(t) λ[t−w(t)]F[w(t)]t2n+1≤t≤t2n+2(49) with initial condition w(t2n+1)=0. Substituting Eq. (48) into Eq. (34), we get the differential equation for v(t): dv(t) dt =λ(t)F[v(t)] µ2B2[t+v(t)] + s′[t+v(t)]−1 t2n+1≤t≤t2n+2(50) with initial condition v(t2n+1)=0, which is consistent with w(t2n+1)=0. Once w(t)is obtained, we get the fluid content Q(t)in the waiting room at time tby Eq. (28). The abandonment rate α(t) is calculated by Eq. (30). The termination epoch t2n+2of the overloaded interval [t2n+1,t2n+2]is determined by the condition in Eq. (10), or equivalently t2n+2=inf{t≥t2n+1:w(t)=0}.(51) 4. Numerical examples for systems with exponentially distributed abandonment time The nonlinear differential equations (49) for w(t)and (50) for v(t)are to be solved numerically when the function F(t)is given. Therefore, a simple form for F(t)does not help much for analytical solution. Nevertheless, we show some numerical examples for systems with the exponentially distributed conversation and ACW times and the exponentially distributed abandonment time. Thus we assume Eq. (38) for the service times along with F(x)=1−e−θx,f(x)=θe−θxx≥0, where θis the constant hazard rate of the abandonment time, and 1/θ is the mean abandonment time. Green et al. [18] mention that it often suffices to work with an exponentially distributed time-toabandon approximation as it was empirically justified. In this case, the rate α(t)of fluid that abandons at time tis proportional to the amount Q(t)of the fluid in the waiting room during overloaded intervals. Indeed, from Eqs. (28) and (30), we have α(t)=w(t) 0 λ(t−x)f(x)dx =θw(t) 0 λ(t−x)e−θxdx =θw(t) 0 λ(t−x)F(x)dx =θQ(t), which is reasonable as θis the rate at which each quantum of fluid in the waiting room abandons. In all the numerical examples in this section, we specifically assume that G1(x)=1−e−5x/4,G2(x)=1−e−5x, F(x)=1−e−2xx≥0 which means that µ1=5 4, µ2=5 (meaning that the total mean service time is 1/µ1+1/µ2=1), and θ=2 in the above formulation. If the unit of time is 1 h, the mean conversation time is 48 min, the mean ACW time is 12 min, and the mean abandonment time is 30 min. 4.1. Sinusoidal input rate and constant number of servers In the first example, we consider a system in which the input rate changes as a sinusoidal function of time twhile the total number of servers is kept constant: λ(t)=100[1+0.6 sin(t)],s(t)=100 t≥0.(52) Green et al. [18] mention that the dynamic behavior of the demand function (for call centers) is reasonably characterized by a sinusoidal function. The input rate function λ(t)in Eq. (52) is used as a base example by Liu and Whitt [19, Section 2]. According to them, by making the mean input rate coincide with the fixed number of servers, we ensure that the system will alternate between underloaded and overloaded states. The result of numerical analysis is shown in Table 1 and Fig. 5. Table 1 shows a sequence of epochs at which the state changes from underloaded to overloaded and vice versa. Fig. 5(a) shows the sinusoidal arrival rate λ(t)given above. Fig. 5(b) plots the fluid content Bi(t)in service i(i=1,2)(dashed line for i=1 and dashed-and-dotted line for i=2) as well as B(t)=B1(t)+B2(t) (solid line) at time t. We observe that B1(t)=80 and that B2(t)= 20 if the system is overloaded at time t. The ratio B1(t)/B2(t)=4 coincides with the ratio of mean service times (1/µ1)/(1/µ2)=4. We also plot the fluid content Q(t)in the waiting room and the total fluid content in the system B(t)+Q(t)with superposed input rate λ(t)in Figs. 5(c) and (d), respectively. We observe that B(t)+Q(t)lags in time behind λ(t). The boundary and potential waiting times are shown in Fig. 5(e). The fluid entrance rate b1(t,0) into service 1 is shown in Fig. 5(f). 88 Y. Kawai, H. Takagi / Operations Research Perspectives 2 (2015) 81–96 Table 1 Epochs of state change in the system with sinusoidal input rate λ(t)=100[1+0.6 sin(t)]and constant number of servers s(t)=100. t0t1t2t3t4t5t6t7t8 0 1.12283 3.60165 6.99994 9.88785 13.28309 16.17103 19.56628 22.45422 (a) Input rate λ(t). (b) Fluid in service Bi(t)and B(t). (c) Waiting fluid Q(t). (d) Total fluid B(t)+Q(t). (e) Boundary waiting time w(t)(thick) and potential waiting time v(t)(thin). (f) Fluid entrance rate into service b1(t,0). Fig. 5. Performance of the system with sinusoidal input rate λ(t)=100[1+0.6 sin(t)]and constant number of servers s(t)=100. 4.2. Constant input rate and sinusoidal number of servers In the second example, we consider a system in which the input rate is constant while the total number of servers changes in time: λ(t)=100,s(t)=100[1+0.6 sin(t)]t≥0, where the non-integer number of servers is thought of as approximation. The result of numerical analysis is shown in Table 2 and Fig. 6. Table 2 shows a sequence of epochs at which the state changes from underloaded to overloaded and vice versa. Fig. 6(a) shows the sinusoidal number of servers s(t)given above. Fig. 6(b) plots the fluid content Bi(t)in service i(i=1,2)as well as B(t)= B1(t)+B2(t)at time twith superposed s(t). We observe that B(t)=s(t)if the system is overloaded at time t. We also plot the fluid content Q(t)in the waiting room and the total fluid content B(t)+Q(t)in Figs. 6(c) and (d), respectively. The boundary and potential waiting times are shown in Fig. 6(e). The fluid entrance rate b1(t,0)into service 1 is shown in Fig. 6(f), which is always positive so that the sufficient condition for feasibility in Eq. (24) is satisfied. 4.3. Identical sinusoidal input rate and number of servers In the third example, we consider a system in which the input rate and the total number of servers change as identical sinusoidal Y. Kawai, H. Takagi / Operations Research Perspectives 2 (2015) 81–96 95 with the feasibility condition b1(t,0)=µ22B22(t)+s′(t)≥0. The initial values B11(t2n+1), B12(t2n+1), B21(t2n+1), and B22(t2n+1) are given from their corresponding values at the end of the preceding underloaded interval [t2n,t2n+1]. The solution can be easily obtained by the conventional method of Laplace transform for solving a set of simultaneous linear differential equations with constant coefficients. The differential equation for the boundary waiting time w(t)is given by dw(t) dt =1−µ22B22(t)+s′(t) λ[t−w(t)]F[w(t)]t2n+1≤t≤t2n+2 with initial condition w(t2n+1)=0. The differential equation for the potential waiting time v(t)is given by dv(t) dt =λ(t)F[v(t)] µ22B22[t+v(t)] + s′[t+v(t)]−1t2n+1≤t≤t2n+2 with initial condition v(t2n+1)=0. Once w(t)is obtained, we get the fluid content Q(t)in the waiting room and the abandonment rate α(t)by Eqs. (28) and (30), respectively. The termination epoch t2n+2of the overloaded interval [t2n+1,t2n+2]is determined by the condition in Eq. (51). 5.3. Numerical example Let us show the numerical results for an Mt/(H2+E2)/s+M system with the sinusoidal input rate function λ(t)and constant number of servers s(t)given in Eq. (52). However we assume that the service 1 has a hyperexponential distribution and that service 2 has a hypo-exponential distribution with parameters µ11 =5 3, µ12 =5 6, µ21 =20, µ22 =20 3,p=2 3,q=1 3 so that g1(x)=10 9e−5 3x+5 18e−5 6x,E[G1] = 4 5,C2 G1=5 4>1, g2(x)=10(e−20 3x−e−20x), E[G2] = 1 5,C2 G2=5 8<1. We also assume the exponentially distributed abandonment time with mean 1/θ =0.5: F(x)=1−e−2xx≥0. The result of numerical analysis is shown in Table 4 and Fig. 12. Table 4 shows a sequence of epochs at which the state changes from underloaded to overloaded and vice versa. These values are only slightly different from those in Table 1 for the system with exponentially distributed conversation and ACW times with the same means. Figs. 12(a) and (b) show the pdf’s of the hyperexponential distribution for the conversation time and the hypo-exponential distribution for the ACW time, respectively. Fig. 12(c) plots the fluid content Bi(t)=Bi1(t)+Bi2(t)in service i(i=1,2)as well as B(t)=B1(t)+B2(t)at time t. These curves are not much different from the corresponding curves in Fig. 5(b). Other performance measures are also shown in Figs. 12(c)–(g). 6. Concluding remarks In this paper, we have studied the fluid approximation to the Mt/GI/st+GI model with two stages of service, each being independent and exponentially distributed, for the call center operation. It is true that the same model can also be handled exactly by the Mt/GI/st+GI model in [19] with a single stage of service consisting of the sum of two exponentially distributed service times. Then, however, one must solve a fixed-point equation for the function b(t,0), which is not straightforward as mentioned in Section 6.2 of [19]. By our method of separating the service to exponentially distributed stages, we can do without bi(t,0)to obtain the fluid content Bi(t). This merit is already pointed out in Section 2 of [19]. It is well-known that most distributions can be approximated precisely enough by a Coxian distribution for which the Laplace transform of pdf is written as a rational function in the transform parameter. Therefore, we should be able to enjoy the above-mentioned merit by decomposing a non-exponential distribution into the serial–parallel combination of exponential distributions. In the present case, the two stages of service correspond to the conversation and ACW times in the real operation of call centers. Therefore, we can know the number of operators engaged in the conversation and the number of those in the ACW individually at each moment. This is another merit that cannot be obtained if the service times of the two stages are treated together. We have shown an example in which the number of servers is determined to cope with the load of only conversation time (service 1). Then the system is always overloaded but it seems to remain stable. On the other hand, if the number of servers is determined to cope with the load of both conversation and ACW times (services 1 and 2), the system is always underloaded. In this case, we need much more servers than in the former example. Therefore, the staffing rule in the former example may be useful as efficient utilization of the resource (servers). Acknowledgments We are grateful to Professor Ward Whitt in the Department of Industrial Engineering and Operations Research of Columbia University, USA, for his valuable comments on the draft of this paper. We also appreciate constructive suggestions by the reviewers of the manuscript who encouraged us to elaborate on the case of non-exponentially distributed conversation and ACW times. References [1] Fischer MJ, Garbin DA, Gharakhanian A. Performance modeling of distributed automatic call distribution systems. Telecommun Syst 1998;9(2):133–52. http://dx.doi.org/10.1023/A:1019139721840. [2] Gans N, Koole G, Mandelbaum A. Telephone call centers: tutorial, review, and research prospects. Manuf Serv Oper Manag 2003;5(2):79–141. http://dx.doi.org/10.1287/msom.5.2.79.16071. [3] Koole G, Mandelbaum A. Queueing models of call centers: an introduction. Annals of Operations Research 2002;113(1–4):41–59. http://dx.doi.org/10.1023/A:1020949626017. [4] Mandelbaum A, Zeltyn S. Service engineering in action: The Palm/Erlang-A queue, with applications to call centers. In: Spath D, Fähnrich K-P, editors. Advances in services innovations. Berlin: Springer; 2007. p. 17–45. [5] Cleveland B, Harne D, editors. Call center operations management. In: Handbook and study guide, Version 2.1. Colorado Springs: ICMI Press; 2004. [6] Harris R.J., Phillips M.J.. Extensions to a model for post call activity in ACD systems, Preprint, December 1989. http://www-ist.massey.ac.nz/rharris/publicationfiles/pre2000/acdsystemsfoatrs.pdf. [7] Jolley W.M., Harris R.J.. Analysis of post-call activity in queueing systems, In: Proceedings of the 9th international teletraffic congress. (Torremolinos); 1979. p. 1–9. [8] Jolley WM, Harris RJ. Analysis and optimal design of queueing systems with post-call activity. Austral Telecommun Res 1981;15(1):12–26. [9] Harris CM, Hoffman KL, Saunders PB. Modeling the IRS telephone taxpayer information system. Oper Res 1987;35(4):504–23. http://dx.doi.org/10.1287/opre.35.4.504. [10] Kawanishi K. On the counting process for a class of Markovian arrival processes with an application to a queueing system. Queueing Sys 2005;49(2):93–122. http://dx.doi.org/10.1007/s11134-005-6478-7. 96 Y. Kawai, H. Takagi / Operations Research Perspectives 2 (2015) 81–96 [11] Kawanishi K. Waiting time distribution of a queueing system with postservice activity. Comput Math Appl 2006;51(2):209–18. http://dx.doi.org/10.1016/j.camwa.2005.11.021. [12] Takagi H, Taguchi Y. Analysis of a queueing model for a call center with impatient customers and after-call work. Int J Pure Appl Math 2014;90(2): 205–37. http://dx.doi.org/10.12732/ijpam.v90i2.10. [13] Phung-Duc T, Kawanishi K. Multiserver retrial queues with after-call work. Numer Algebra Control Optim 2011;1(4):639–56. http://dx.doi.org/10.3934/naco.2011.1.639. [14] Hall RW. Queueing systems. In: For services and manufacturing. Prentice-Hall; 1991. [15] Bassamboo A, Randhawa RS. On the accuracy of fluid models for capacity sizing in queueing systems with impatient customers. Oper Res 2010;58(5): 1398–413. http://dx.doi.org/10.1287/opre.1100.0815. [16] Whitt W. Fluid models for multiserver queues with abandonments. Oper Res 2006;54(1):37–54. http://dx.doi.org/10.1287/opre.1050.0227. [17] Feldman Z, Mandelbaum A, Massey WA, Whitt W. Staffing of time-varying queues to achieve time-stable performance. Manag Sci 2008;54(2):324–38. http://dx.doi.org/10.1287/mnsc.1070.0821. [18] Green LV, Kolesar PJ, Whitt W. Coping with time-varying demand when setting staffing requirements for a service system. Prod Oper Manag 2007;16(1): 13–39. http://dx.doi.org/10.1111/j.1937-5956.2007.tb00164.x. [19] Liu Y, Whitt W. The Gt/GI/st+GI many-server fluid queue. Queueing Syst 2012;71(4):405–44. http://dx.doi.org/10.1007/s11134-012-9291-0. Longer online version with appendix in http://www.columbia.edu/ww2040/allpapers.html. [20] Liu Y, Whitt W. Fluid models for multiserver queues with abandonments. Oper Res 2006;54(1):37–54. http://dx.doi.org/10.1287/opre.1050.0227. [21] Liu Y, Whitt W. A network of time-varying many-server fluid queues with customer abandonment. Oper Res 2011;59(4):835–46. http://dx.doi.org/10.1287/opre.1110.0942. [22] Liu Y, Whitt W. Algorithms for time-varying networks of many-server fluid queues. INFORMS J Comput 2014;26(1):59–73. http://dx.doi.org/10.1287/ijoc.1120.0547. [23] Kawai Y. Fluid approximation analysis for the queueing model of a call center [Master thesis], Tsukuba, Japan: Master of Engineering, Graduate School of Systems, Information and Engineering, University of Tsukuba; 2013 [in Japanese]. [24] Eick S, Massey WA, Whitt W. The physics of the Mt/G/∞queue. Oper Res 1993;41(4):731–42. http://dx.doi.org/10.1287/opre.41.4.731. [25] Eick S, Massey WA, Whitt W. Mt/G/∞queues with sinusoidal arrival rates. Manag Sci 1993;39(2):241–52. http://dx.doi.org/10.1287/mnsc.39.2.241. [26] Bertsimas D, Mourtzinou G. Transient laws of non-stationary queueing systems and their applications. Queueing Syst 1997;25(1):115–55. http://dx.doi.org/10.1023/A:1019100301115.