Network Stochastic Call Center Model Tatiana V. Rusilko Yanka Kupala State University of Grodno, 22 Ozheshko St, Grodno, 230023, Belarus Abstract A closed exponential queueing network with impatient customers is considered as a model of a call center. The process of handling calls in the call center is described. A detailed description of the model is given. The stay of the customer at a specific network node and its routing between the nodes correspond to the customer’s call status in the call center and the process of its routing between telephone agents of different levels. The process of changing the number of customers at the nodes is studied under the asymptotic assumption of a large number of customers. It is proved that its probability density function satisfies the Fokker–Planck– Kolmogorov equation. The system of differential equations for the expected value of this process is derived. The solution of this system allows predicting the dynamics of the expected number of customers in the network nodes and regulating the parameters of the call center operation. The asymptotic technique used is applicable both in transient and steady-state. The problem of finding the optimal staffing list of agents to maximize call center performance is investigated. The areas of implementation of research results are the design of call centers and improving their performance. Keywords 1 Call center, queueing network, impatient customers, mathematical modeling 1. Introduction A call center is a specialized organization or a centralized department located within a company in which telephone calls from customers are handled by a team of advisors, otherwise known as agents. Call centers can handle inbound and outbound calls via voice communication channels, although some may specialize in one or the other. The article [1] concentrates on review the state of research on telephone call centers, this article begins with a tutorial on how call centers function and proceed to survey academic research devoted to the management of their operations. The quality and operational efficiency of these telephone services must be exceptional to meet the needs of the customers [2]. Agents of a large best-practice call center have to cater to thousands of phone callers per hour. The waiting time of delayed calls has not to exceed a few seconds. To achieve high levels of service quality and efficiency it is a necessary understanding of the scientific principles underlying call-center operations. Undoubtedly sophisticated approaches are needed to accurately describe the reality of contact-center operations, and mathematical modeling of this reality can improve call center performance significantly [1, 3, 4]. Call center capabilities also include the accumulation of a call queue in standby mode, waiting time information, call distribution between agents, call classification, and interaction strategy determination depending on the client's status. The operational service level of call centers is typically quantified in terms of some efficiency indices or performance measures. Among them, the following can be noted: workload, placing calls by agents, expected call processing time, percentage of successfully processed calls, duration of the service cycle of one call, the fraction of customers that abandon the queue before being served, queue lengths, waiting time. Proceedings of VI International Scientific and Practical Conference Distance Learning Technologies (DLT–2021), September 20-22, 2021, Yalta, Crimea EMAIL: tatiana.rusilko@gmail.com ORCID: 0000-0002-4880-0619 ©️ 2021 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) 91 Having analyzed the above organization of call centers, we conclude that queueing theory can be used to analyze their efficiency, predict the load and optimize the structure. Performance measures listed above also underscore the natural fit between queueing models and call centers [5]. The emergence of the queueing theory is associated with the ordering and optimization of the operation of telephone exchanges. Queueing systems, in particular, Markovian ones, are widely used for call center modeling and performance analysis [6, 7, 8, 9, 10]. Queueing systems with impatient customers [7, 10], retrial queueing systems [6, 9], and many others have originated from the need for mathematical modeling of the behavior of customers of telephone services. Call centers often have a branched or multilevel structure. The functions of the first-level call center agents include call registration and provision of general information. When the topic of a client request goes beyond the specialization of current level agents, the client is redirected to a higher level, whose specialists are competent in resolving the issue. Thus, calls are routed according to their specificity between different levels of agents. Hence the call center organization structure can be graphically depicted by a network diagram showing the hierarchy of call center agents, routing, and call forwarding according to their specifics between different levels of advisors. Thus, queueing networks are suitable for modeling multilevel or multi-stage call centers. A queueing network is a collection of interconnected queueing systems, between which customers (jobs, requests) circulate. Each queueing system, which is also called a node, has at least one server and a queue. The main classification of queueing networks is based on the location of the customer source, the flow of arriving customers, the service time probability distribution, the scheduling policy at queueing nodes, and the customer routing. The purpose of this paper is call center mathematical modeling, optimization, and analysis of the effectiveness of the call processing process using closed exponential queueing networks with impatient customers. The network is studied under the asymptotic assumption of a large number of customers [11, 12]. Thus, the greater the number of customers in the network, the greater the precision in the calculation [13]. 2. Model Description. Problem Formulation Our focus in this article is on the inbound call center of a large company with a limited number of clients. We denote this number of customers by 𝐾, 𝐾 ≫ 1. Inbound call centers handle incoming calls that are initiated by outside callers calling into a center. Typically, these types of centers provide customer support, help-desk services, reservation and sales support for airlines and hotels, and order- taking functions for catalog and Web-based merchants. A closed queueing network is proposed as a model. The network consists of 𝑛 + 1 Markov 𝑚𝑖 -linear queueing systems 𝑆𝑖 , 𝑖 = 1, 𝑛 + 1, and 𝐾-linear hypothetical system 𝑆0 , which plays the role of an external environment and is considered as a dependable source of arriving customers. Each call to the call center corresponds to the customer's arrival in the queueing network. The stay of the customer at a specific network node and its routing between the nodes correspond to the customer’s call status in the call center and the process of its routing between telephone agents of different levels. Let's describe the process of handling calls in the call center. Consider a company’s two-level call center that handles calls on 𝑛 topics. Each customer can be in one of the following states:  𝑆0 – no need to turn to the call center,  𝑆𝑛+1 – call is required service by one of 𝑚𝑛+1 first-level agents who register and distribute calls depending on their topic,  𝑆𝑖 – call is required service by one of the 𝑚𝑖 second-level agents consulting on the 𝑖th topic, 𝑖 = 1, 𝑛. Accordingly, each customer of the queueing network can be located in one of the queueing systems 𝑆𝑖 , 𝑖 = 0, 𝑛 + 1. The customer transition from state 𝑆0 to state 𝑆𝑛+1 corresponds to a customer's call to the call center, which occurs at a random time, regardless of the state of other customers, but depends on the time of day. The probability of this transition in the time interval [𝑡, 𝑡 + ∆𝑡] is equal to 𝜆(𝑡)Δ𝑡 + 𝜊(Δ𝑡), where 92 𝜆(𝑡) is the arrival rate (the expected time between each call arriving). The inbound call flow is non- stationary. From state 𝑆𝑛+1 the call is routed to the state 𝑆𝑖 with probability 𝑝(𝑛+1)𝑖 , 𝑖 = 1, 𝑛, or with probability 𝑝(𝑛+1)0 to stage 𝑆0 in case of an erroneous customer request. After the end of call service in the state 𝑆𝑖 , 𝑖 = 1, 𝑛, the customer goes to 𝑆0 with probability 𝑝𝑖0 , if his request is completely fulfilled or redirected to specialists on another topic to stage 𝑆𝑗 , 𝑗 = 1, 𝑛, with probability 𝑝𝑖𝑗 , 𝑖 ≠ 𝑗, if the agents of the state 𝑆𝑖 have a lack of competence to answer the customer's questions. The served customer transition matrix is stochastic (𝑝𝑖𝑗 )(𝑛+2)×(𝑛+2) with the following nonzero entries: 𝑛 𝑝0(𝑛+1) = 1; 𝑝(𝑛+1)𝑖 , 𝑖 = 0, 𝑛, ∑ 𝑝(𝑛+1)𝑖 = 1; 𝑖=0 𝑛 𝑝𝑖𝑗 , 𝑖 = 1, 𝑛, 𝑗 = 0, 𝑛, 𝑖 ≠ 𝑗, ∑ 𝑝𝑖𝑗 = 1. 𝑗=0 The customer routing in the network model of the call center is depicted with a solid line in Figure 1. Figure 1: The state-transition diagram of served and impatient customers If at least one of the agents of the system 𝑆𝑖 is idle, then the incoming call immediately begins to be served at the moment of transition to state 𝑆𝑖 , 𝑖 = 1, 𝑛 + 1. Calls are served on a first-come, first-served basis. For the queueing network, this means that service requests are selected accordingly to the discipline FIFO. All 𝑚𝑖 agents (servers) at the state (node) 𝑆𝑖 are identical and they have exponentially distributed service time with rate 𝜇𝑖 , i.e., 𝜇𝑖−1 is the mean service time, 𝑖 = 1, 𝑛 + 1. If all 𝑚𝑖 agents in 𝑆𝑖 are busy, then the customer waits in an unlimited queue. However, the client waiting time in the queue of the system 𝑆𝑖 is limited by an exponential random variable 𝜏 with an expected value of 𝜂𝑖−1 , 𝑖 = 1, 𝑛 + 1. Under this condition, customers in the queueing network are called impatient. If during time 𝜏 the call is not answered, then it is lost, i.e., with a probability equal to one, the customer after this time passes into the external environment 𝑆0 . The impatient (unfulfilled) customer transition matrix is (𝑞𝑖𝑗 )(𝑛+1)×(𝑛+2) with the following nonzero entries: 𝑞𝑖0 = 1, 𝑖 = 1, 𝑛 + 1. Let 𝑘𝑖 (𝑡) denote the number of customers in the state 𝑆𝑖 at the time 𝑡, 𝑖 = 1, 𝑛 + 1. It's obvious that ∑𝑛+1 𝑖=0 𝑘𝑖 (𝑡) = 𝐾. The allocation of customers according to possible states at time 𝑡 fully describes the situation in the call center at that time. Accordingly, the allocation of customers by queueing nodes completely determines the state of the queueing network. Thus, the state of the network model understudy at time 𝑡 is represented by a vector 𝑘(𝑡) = (𝑘1 (𝑡), 𝑘2 (𝑡), . . . , 𝑘𝑛+1 (𝑡)). Taking into account the above-described, the process 𝑘(𝑡) of changing the number of customers at queueing nodes is a continuous-time Markov chain with a finite state space. It is clear, the number of customers in the external node 𝑆0 is 𝑘0 (𝑡) = 𝐾 − ∑𝑛+1 𝑖=1 𝑘𝑖 (𝑡). 93 The following problem is of interest: taking into account the benefits of customer calls and considering agents wage fund, it is required to find such a staffing list (a set of servers in the network nodes) that would ensure the maximum call center performance, provided that 𝐾, 𝜆(𝑡), 𝜇𝑖 , 𝜂𝑖 , 𝑝𝑖𝑗 are fixed 𝑖 = 1, 𝑛 + 1, 𝑗 = 0, 𝑛 + 1. Let us introduce the following notation:  𝐹𝑖 is the monetary effect (or some kind of efficiency) generated by servicing one call in state 𝑆𝑖 per unit of time, 𝑖 = 1, 𝑛 + 1;  𝐿𝑖 is the wage of one agent in the system 𝑆𝑖 per unit of time, 𝑖 = 1, 𝑛 + 1;  𝐹0 is the monetary effect created by one customer when he does not need to turn to the call center (it can be equal to zero). Then the total monetary effect of the call center at time 𝑡 is given by the following expression 𝑛+1 𝑛+1 𝑛+1 𝐹(𝑡) = 𝐹0 (𝐾 − ∑ min(𝑚𝑖 , 𝑘𝑖 (𝑡))) + ∑ 𝐹𝑖 min(𝑚𝑖 , 𝑘𝑖 (𝑡)) − ∑ 𝐿𝑖 𝑚𝑖 . (1) 𝑖=1 𝑖=1 𝑖=1 Since 𝑘𝑖 (𝑡) are components of a random process, the productivity 𝐹(𝑡) is also a random process. Applying the mathematical expectation operator E to both sides of the expression (1), we obtain the following expression for the expected monetary effect of a call center: 𝑛+1 𝑛+1 E𝐹 (𝑡) = 𝐹0 𝐾 + ∑(𝐹𝑖 − 𝐹0 )min(𝑚𝑖 , E𝑘𝑖 (𝑡)) − ∑ 𝐿𝑖 𝑚𝑖 . (2) 𝑖=1 𝑖=1 It is assumed that call center operation is worthwhile if E𝐹 (𝑡) > 0, and its efficiency increases with E𝐹 (𝑡). It is interesting to consider the problem of maximizing the functional (2) by the number of servers 𝑚𝑖 in queueing systems 𝑆𝑖 , 𝑖 = 1, 𝑛 + 1, provided that the statistical characteristics of customer behavior are fixed. In the simplest case of the assumption about the absence of queues, the mathematical problem of maximizing the average efficiency at the time 𝑡 is given by 𝑛+1 𝑛+1 E𝐹 (𝑡) = 𝐹0 𝐾 + ∑(𝐹𝑖 − 𝐹0 )E𝑘𝑖 (𝑡) − ∑ 𝐿𝑖 𝑚𝑖 → max , 𝑖=1 𝑖=1 𝑚𝑖 ,𝑖=1,𝑛+1 (3) {E𝑘𝑖 (𝑡) ≤ 𝑚𝑖 , 𝑖 = 1, 𝑛 + 1. It is clear that to solve this problem, first of all, it is necessary to determine the probability distribution of the process 𝑘(𝑡). Only then can the expected values E𝑘𝑖 (𝑡) be determined, 𝑖 = 1, 𝑛 + 1. The distribution of 𝑘(𝑡) parametrically depends on the number of servers in queueing systems. 3. Asymptotic Analysis of the Network Model Asymptotic analysis implies an approximation method of queueing network study under the assumption of a large number of customers. The mathematical approach used in this paper is based on a discrete model of a continuous Markov process described in many books on the theory of diffusion Markov processes [14, 15]. 𝑖 Let us introduce a (𝑛 + 1)-vector of the form 𝐼𝑖 = (0, ⏞ 0, . . . , 0, 1 , 0, . . . , 0) and Heaviside function 1, 𝑥 > 0, 𝐻(𝑥) = { As mentioned above, the process 𝑘(𝑡) is a continuous-time Markov chain with a 0, 𝑥 ≤ 0. finite state space. The assumptions made in the formulation of the problem determine that in the short time Δ𝑡 the Markov process 𝑘(𝑡) = (𝑘, 𝑡) can make one of the following transitions to the state (𝑘, 𝑡 + Δ𝑡):  from the state (𝑘 − 𝐼𝑛+1 , 𝑡) the process transfer to (𝑘, 𝑡 + Δ𝑡) with probability 𝑛+1 𝜆(𝑡) (𝐾 − ∑ 𝑘𝑖 (𝑡) + 1) Δ𝑡 + 𝜊(Δ𝑡), 𝑖=1 that corresponds to the customer transition from state S0 to state Sn+1 ;  from the state (𝑘 + 𝐼𝑖 , 𝑡) the process transfer to (𝑘, 𝑡 + Δ𝑡) with probability 𝜂𝑖 (𝑘𝑖 (𝑡) + 1 − 𝑚𝑖 )𝐻(𝑘𝑖 (𝑡) + 1 − 𝑚𝑖 )Δ𝑡 + 94 +𝜇𝑖 𝑝𝑖0 min( 𝑚𝑖 , 𝑘𝑖 (𝑡) + 1)Δ𝑡 + 𝜊(Δ𝑡), that corresponds to the impatient customer transition from state 𝑆𝑖 to state S0 , 𝑖 = 1, 𝑛 + 1;  from the state (𝑘 + 𝐼𝑖 − 𝐼𝑗 , 𝑡) the process transfer to (𝑘, 𝑡 + Δ𝑡) with probability 𝜇𝑖 𝑝𝑖𝑗 min( 𝑚𝑖 , 𝑘𝑖 (𝑡) + 1)Δ𝑡 + 𝜊(Δ𝑡), that corresponds to the served customer transition from state 𝑆𝑖 to state S𝑗 , 𝑖 = 1, 𝑛 + 1, 𝑗 = 1, 𝑛;  from the state (𝑘, 𝑡) to (𝑘, 𝑡 + Δ𝑡) with probability 𝑛+1 𝑛+1 1 − (𝜆(𝑡) (𝐾 − ∑ 𝑘𝑖 (𝑡)) + ∑ 𝜇𝑖 𝑝𝑖𝑗 𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡)) + 𝑖=1 𝑖,𝑗=1 𝑛+1 + ∑(𝜇𝑖 𝑝𝑖0 𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡)) + 𝜂𝑖 (𝑘𝑖 (𝑡) − 𝑚𝑖 )𝐻(𝑘𝑖 (𝑡) − 𝑚𝑖 ))) Δ𝑡 + 𝜊(Δ𝑡), 𝑖=1 that corresponds to no customers transfer;  from other states the process transfer to (𝑘, 𝑡 + Δ𝑡) with probability 𝜊(Δ𝑡). The probabilities 𝑝𝑖𝑗 are entries of the state transition matrix (𝑝𝑖𝑗 )(𝑛+2)×(𝑛+2) that was specified in the model description section. For example, 𝑝𝑖(𝑛+1) = 0, 𝑖 = 1, 𝑛 + 1. Having regard to the transitions listed above in the short time Δ𝑡, using the law of total probability, the following system of equations is valid for the probability 𝑃(𝑘, 𝑡) = 𝑃(𝑘(𝑡) = 𝑘): 𝑛+1 𝑃(𝑘, 𝑡 + 𝛥𝑡) = 𝑃(𝑘 − 𝐼𝑛+1 , 𝑡)𝜆(𝑡) (𝐾 − ∑ 𝑘𝑖 (𝑡) + 1) Δ𝑡 + 𝑖=1 𝑛+1 + ∑ 𝑃(𝑘 + 𝐼𝑖 , 𝑡)(𝜇𝑖 𝑝𝑖0 𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡) + 1) + 𝜂𝑖 (𝑘𝑖 (𝑡) + 1 − 𝑚𝑖 )𝐻(𝑘𝑖 (𝑡) + 1 − 𝑚𝑖 ))Δ𝑡 + 𝑖=1 𝑛+1 + ∑ 𝑃(𝑘 + 𝐼𝑖 − 𝐼𝑗 , 𝑡)𝜇𝑖 𝑝𝑖𝑗 𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡) + 1)Δt + 𝑖,𝑗=1 𝑛+1 𝑛+1 +𝑃(𝑘, 𝑡) (1 − (𝜆(𝑡) (𝐾 − ∑ 𝑘𝑖 (𝑡)) + ∑ 𝜇𝑖 𝑝𝑖𝑗 𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡)) + 𝑖=1 𝑖,𝑗=1 𝑛+1 + ∑(𝜇𝑖 𝑝𝑖0 𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡)) + 𝜂𝑖 (𝑘𝑖 (𝑡) − 𝑚𝑖 )𝐻(𝑘𝑖 (𝑡) − 𝑚𝑖 ))) Δ𝑡) + 𝜊(Δ𝑡). 𝑖=1 We set 𝛥𝑖0 𝑃(𝑘, 𝑡) = 𝑃(𝑘 + 𝐼𝑖 , 𝑡) − 𝑃(𝑘, 𝑡), 𝛥0(𝑛+1) 𝑃(𝑘, 𝑡) = 𝑃(𝑘 − 𝐼𝑛+1 , 𝑡) − 𝑃(𝑘, 𝑡), 𝛥𝑖𝑗 𝑃(𝑘, 𝑡) = 𝑃(𝑘 + 𝐼𝑖 − 𝐼𝑗 , 𝑡) − 𝑃(𝑘, 𝑡), 𝑖 = 1, 𝑛 + 1, 𝑗 = 1, 𝑛. On the assumption that Δ𝑡 → 0, we pass to the Kolmogorov difference-differential equation for these probabilities: 𝑛+1 𝜕𝑃(𝑘, 𝑡) = 𝜆(𝑡) (𝐾 − ∑ 𝑘𝑖 (𝑡)) 𝛥0(𝑛+1) 𝑃(𝑘, 𝑡) + 𝜆(𝑡)𝑃(𝑘 − 𝐼𝑛+1 , 𝑡) + 𝜕𝑡 𝑖=1 𝑛+1 + ∑(𝜇𝑖 𝑝𝑖0 𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡)) + 𝜂𝑖 (𝑘𝑖 (𝑡) − 𝑚𝑖 )𝐻(𝑘𝑖 (𝑡) − 𝑚𝑖 ))𝛥𝑖0 𝑃(𝑘, 𝑡) + 𝑖=1 𝑛+1 + ∑ 𝜇𝑖 𝑝𝑖0 (𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡) + 1) − 𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡)))𝑃(𝑘 + 𝐼𝑖 , 𝑡) + (4) 𝑖=1 95 𝑛+1 + ∑ 𝜂𝑖 𝑝𝑖0 ((𝑘𝑖 (𝑡) + 1 − 𝑚𝑖 )𝐻(𝑘𝑖 (𝑡) + 1 − 𝑚𝑖 ) − (𝑘𝑖 (𝑡) − 𝑚𝑖 )𝐻(𝑘𝑖 (𝑡) − 𝑚𝑖 ))𝑃(𝑘 + 𝐼𝑖 , 𝑡) + 𝑖=1 𝑛+1 + ∑ 𝜇𝑖 𝑝𝑖𝑗 𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡))𝛥𝑖𝑗 𝑃(𝑘, 𝑡) + 𝑖,𝑗=1 𝑛+1 + ∑ 𝜇𝑖 𝑝𝑖𝑗 (𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡) + 1) − 𝑚𝑖𝑛( 𝑚𝑖 , 𝑘𝑖 (𝑡)))𝑃(𝑘 + 𝐼𝑖 − 𝐼𝑗 , 𝑡). 𝑖,𝑗=1 The last equation defies analytical solution for large 𝑛. Call centers undoubtedly handle a large number of customer calls. In connection with this, consider the important asymptotic case of a large number of customers 𝐾 ≫ 1 and find the probability distribution of the state vector of the call center model with an accuracy of 𝛰(𝜀 2 ), where 𝜀 = 1⁄𝐾. The technique of asymptotic analysis of queueing networks described in the papers [11, 12, 16] is used. Suppose that we are interested in the properties of the process 𝜉(𝑡) = 𝑘(𝑡)/𝐾 as 𝐾 becomes very large. Vector 𝜉(𝑡) indicates which part of the company's customers contacted the call center and how the calls were distributed by model states at time 𝑡. During the time Δ𝑡 → 0, the possible changes in the vector 𝜉(𝑡) are 𝑒𝑖 , where 𝑒𝑖 = 𝐼𝑖 𝜀 = 𝐼𝑖 /𝐾. In the case when 𝐾 → ∞ we have 𝑒𝑖 → 0 and the vector 𝑘(𝑡) 𝑘1 (𝑡) 𝑘2 (𝑡) 𝑘𝑛+1 (𝑡) 𝜉(𝑡) = =( , ,…, ) 𝐾 𝐾 𝐾 𝐾 will be continuous-time continuous-state Markov processes with a probability density function 𝑃(𝑥1 ≤ 𝜉1 (𝑡) < 𝑥1 + 𝜀, 𝑥2 ≤ 𝜉2 (𝑡) < 𝑥2 + 𝜀, . . . , 𝑥𝑛+1 ≤ 𝜉𝑟 (𝑡) < 𝑥𝑛+1 + 𝜀) 𝑝(𝑥, 𝑡) = lim = 𝜀→0 𝜀 𝑛+1 𝑃(𝐾𝑥1 ≤ 𝑘1 (𝑡) < 𝐾𝑥1 + 1, 𝐾𝑥2 ≤ 𝑘2 (𝑡) < 𝐾𝑥2 + 1, . . . , 𝐾𝑥𝑛+1 ≤ 𝑘𝑛+1 (𝑡) < 𝐾𝑥𝑛+1 + 1) = lim , 𝜀→0 𝜀 𝑛+1 i.e., 𝐾 𝑛+1 𝑃(𝑘, 𝑡) → 𝑝(𝑥, 𝑡), (5) 𝐾→∞ where 𝑥 ∈ 𝑋, 𝑛+1 𝑋 = {𝑥 = (𝑥1 , 𝑥2 , . . . , 𝑥𝑛+1 ): 𝑥𝑖 ≥ 0, 𝑖 = 1, 𝑛 + 1, ∑ 𝑥𝑖 ≤ 1}. 𝑖=1 Realizing the asymptotic transition (5) for the equation (4), we obtain the following partial differential equation: 𝑛+1 𝜕𝑝(𝑥, 𝑡) = 𝜆(𝑡)𝐾 (1 − ∑ 𝑥𝑖 (𝑡)) 𝛥0(𝑛+1) 𝑝(𝑥, 𝑡) + 𝜆(𝑡)𝑝(𝑥 − 𝑒𝑛+1 , 𝑡) + 𝜕𝑡 𝑖=1 𝑛+1 +𝐾 ∑(𝜇𝑖 𝑝𝑖0 𝑚𝑖𝑛( 𝑙𝑖 , 𝑥𝑖 (𝑡)) + 𝜂𝑖 (𝑥𝑖 (𝑡) − 𝑙𝑖 )𝐻(𝑥𝑖 (𝑡) − 𝑙𝑖 ))𝛥𝑖0 𝑝(𝑥, 𝑡) + (6) 𝑖=1 𝑛+1 𝑛+1 𝜕 𝑚𝑖𝑛( 𝑙𝑖 , 𝑥𝑖 ) 𝜕(𝑥𝑖 (𝑡) − 𝑙𝑖 )𝐻(𝑥𝑖 (𝑡) − 𝑙𝑖 )) + ∑ 𝜇𝑖 𝑝𝑖0 𝑝(𝑥 + 𝑒𝑖 , 𝑡) + ∑ 𝜂𝑖 𝑝𝑖0 𝑝(𝑥 + 𝑒𝑖 , 𝑡) + 𝜕𝑥𝑖 𝜕𝑥𝑖 𝑖=1 𝑖=1 𝑛+1 𝑛+1 𝜕 𝑚𝑖𝑛( 𝑙𝑖 , 𝑥𝑖 ) +𝐾 ∑ 𝜇𝑖 𝑝𝑖𝑗 𝑚𝑖𝑛( 𝑙𝑖 , 𝑥𝑖 (𝑡))𝛥𝑖𝑗 𝑝(𝑥, 𝑡) + ∑ 𝜇𝑖 𝑝𝑖𝑗 𝑝(𝑥 + 𝑒𝑖 − 𝑒𝑗 , 𝑡). 𝜕𝑥𝑖 𝑖,𝑗=1 𝑖,𝑗=1 If 𝑝(𝑥, 𝑡) is twice continuously differentiable function for 𝑥, then 96 𝜕𝑝(𝑥, 𝑡) 𝜀 2 𝜕 2 𝑝(𝑥, 𝑡) 𝑝(𝑥 ± 𝑒𝑖 , 𝑡) = 𝑝(𝑥, 𝑡) ± 𝜀 + + 𝑜(𝜀 2 ), 𝜕𝑥𝑖 2 𝜕𝑥𝑖2 𝜕𝑝(𝑥, 𝑡) 𝜕𝑝(𝑥, 𝑡) 𝑝(𝑥 + 𝑒𝑖 − 𝑒𝑗 , 𝑡) = 𝑝(𝑥, 𝑡) + 𝜀 ( − )+ 𝜕𝑥𝑖 𝜕𝑥𝑗 𝜀 2 𝜕 2 𝑝(𝑥, 𝑡) 𝜕 2 𝑝(𝑥, 𝑡) 𝜕 2 𝑝(𝑥, 𝑡) + ( − 2 + ) + 𝑜(𝜀 2 ), 𝑖, 𝑗 = 1, 𝑛 + 1. 2 𝜕𝑥𝑖2 𝜕𝑥𝑖 𝜕𝑥𝑗 𝜕𝑥𝑗2 Substituting this series into equation (6) and grouping terms, we obtain that 𝑝(𝑥, 𝑡) satisfies the multidimensional Fokker–Planck–Kolmogorov equation (7): 𝑛+1 𝑛+1 𝜕𝑝(𝑥, 𝑡) 𝜕 𝜀 𝜕2 = −∑ (𝐴𝑖 (𝑥, 𝑡)𝑝(𝑥, 𝑡)) + ∑ (𝐵 (𝑥, 𝑡)𝑝(𝑥, 𝑡)) , (7) 𝜕𝑡 𝜕𝑥𝑖 2 𝜕𝑥𝑖 𝜕𝑥𝑗 𝑖𝑗 𝑖=1 𝑖,𝑗=1 with drifts 𝑛+1 𝐴𝑖 (𝑥, 𝑡) = ∑ 𝜇𝑗 (𝑝𝑗𝑖 − 𝛿𝑗𝑖 )𝑚𝑖𝑛( 𝑙𝑗 , 𝑥𝑗 (𝑡)) − 𝜂𝑖 (𝑥𝑖 (𝑡) − 𝑙𝑖 )𝐻(𝑥𝑖 (𝑡) − 𝑙𝑖 ) + (8) 𝑗=1 𝑛+1 +𝜆(𝑡)𝑝0𝑖 (1 − ∑ 𝑥𝑗 (𝑡)), 𝑗=1 and diffusion coefficients 𝑛+1 𝐵𝑖𝑖 (𝑥, 𝑡) = ∑ 𝜇𝑗 (𝑝𝑗𝑖 + 𝛿𝑗𝑖 )𝑚𝑖𝑛( 𝑙𝑗 , 𝑥𝑗 (𝑡)) + 𝜂𝑖 (𝑥𝑖 (𝑡) − 𝑙𝑖 )𝐻(𝑥𝑖 (𝑡) − 𝑙𝑖 ) + 𝑗=1 𝑛+1 +𝜆(𝑡)𝑝0𝑖 (1 − ∑ 𝑥𝑘 (𝑡)), 𝑗=1 𝐵𝑖𝑗 (𝑥, 𝑡) = −𝜇𝑖 𝑝𝑖𝑗 𝑚𝑖𝑛( 𝑙𝑖 , 𝑥𝑖 (𝑡)), 𝑖 ≠ 𝑗; 𝛿𝑖𝑗 is the Kronecker delta, 𝑖, 𝑗 = 1, 𝑛 + 1. The error in this approximation is no more than 𝜀 2 = (1/𝐾)2 . It is known [14, 17] that the expected value E𝜉𝑖 (𝑡) are determined with an accuracy of 𝛰(1/𝐾 2 ) from a system of ordinary differential equations: 𝑑E𝜉𝑖 (𝑡) = 𝐴𝑖 (E𝜉𝑖 (𝑡), 𝑡), 𝑖 = 1, 𝑛 + 1, 𝑑𝑡 where 𝐴𝑖 (𝑥, 𝑡) is drift given by formula (8). The proof of this statement, based on the application of characteristic functions, is presented in the article [18]. The system of differential equations for the components of vector 𝑘(𝑡) = 𝐾𝜉(𝑡) is given by 𝑑E𝑘𝑖 (𝑡) 1 = 𝐾𝐴𝑖 ( E𝑘𝑖 (𝑡), 𝑡) , 𝑖 = 1, 𝑛 + 1. (9) 𝑑𝑡 𝐾 Taking into account (9), (8), and the transition matrix (𝑝𝑖𝑗 )(𝑛+2)×(𝑛+2) of the network model, the expected number of customers E𝑘𝑖 (𝑡) in state 𝑆𝑖 , 𝑖 = 1, 𝑛 + 1, can be found by solving the following system: 𝑛+1 𝑑E𝑘𝑖 (𝑡) = ∑ 𝜇𝑗 𝑝𝑗𝑖 𝑚𝑖𝑛( 𝑚𝑗 , E𝑘𝑗 (𝑡)) − 𝜇𝑖 𝑚𝑖𝑛( 𝑚𝑖 , E𝑘𝑖 (𝑡)) − 𝑑𝑡 𝑗=1 −𝜂𝑖 (E𝑘𝑖 (𝑡) − 𝑚𝑖 )𝜃(E𝑘𝑖 (𝑡) − 𝑚𝑖 ), 𝑖 = 1, 𝑛, (10) 𝑛+1 𝑑E𝑘𝑛+1 (𝑡) = 𝜆(𝑡) (𝐾 − ∑ E𝑘𝑖 (𝑡)) − 𝜇𝑛+1 𝑚𝑖𝑛( 𝑚𝑛+1 , E𝑘𝑛+1 (𝑡)) − 𝑑𝑡 𝑖=1 −𝜂𝑛+1 (E𝑘𝑛+1 (𝑡) − 𝑚𝑛+1 )𝜃(E𝑘𝑛+1 (𝑡) − 𝑚𝑛+1 ). The expected number of customers in the external environment is 97 𝑛+1 E𝑘0 (𝑡) = 𝐾 − ∑ E𝑘𝑖 (𝑡). 𝑖=1 Solution of the system (10) with a certain initial condition, firstly, makes it possible to predict the average number of calls in each model state, and, secondly, it is necessary to solve the problem of maximizing the call center performance (3). To solve (10), it is possible to use numerical methods or to solve the system analytically in the linearity regions of the right-hand side. Having solved (10) for different sets of the number of agents 𝑚 = (𝑚1 , 𝑚2 , . . . , 𝑚𝑛+1 ), substituting this solution into (2) and comparing the values of the expected productivity E𝐹 (𝑡) for each 𝑚 = (𝑚1 , 𝑚2 , . . . , 𝑚𝑛+1 ), we find the optimal set 𝑚 for which E𝐹 (𝑡) takes the largest value at a fixed time 𝑡. It is well established that under certain conditions, the network reaches "steady-state", the stationary value of the expected productivity is E𝐹 = lim E𝐹 (𝑡). The steady-state means E𝑘𝑖 = lim E𝑘𝑖 (𝑡), where 𝑡→∞ 𝑡→∞ 𝑖 = 1, 𝑛 + 1, can be found from a system of algebraic balance equations of the form (10), where the left-hand sides are equal to zero. In this case, the arrival rate should be constant: 𝜆(𝑡) = 𝜆. The relative number E𝑘𝑖 /𝐾 is numerically equal to the probability of the customer staying at the state (node) 𝑆𝑖 , 𝑖 = 1, 𝑛 + 1, in the steady-state. In the steady-state, problem (3) turns into a linear programming problem. 4. Numerical Example Consider the company call center that handles calls on three topics, 𝑛 = 3. Let the number of customers of the company be equal to 𝐾 = 20 000. Let the work of the call center be specified by the following parameters: 0 0 0 0 1 0.8 0 0.1 0.1 0  the transition matrix is 𝑃 = 0.9 0.03 0 0.07 0 , 0.75 0.1 0.15 0 0 ( 0.1 0.45 0.3 0.15 0) 2𝜋𝑡  the arrival rate is 𝜆(𝑡) = 0.002sin ( + 2) + 0.0083 or 𝜆 =0.01; 24  the number of servers in network nodes 𝑚1 = 3, 𝑚2 = 4, 𝑚3 = 3, 𝑚4 = 2;  the service rates 𝜇1 = 40, 𝜇2 = 20, 𝜇3 = 15, 𝜇4 = 135;  the mean waiting time 𝜂1−1 = 2, 𝜂2−1 = 2, 𝜂3−1 = 2.5, 𝜂4−1 = 10;  the initial placement of customers E𝑘0 (𝑡) = 20 000, E𝑘𝑖 (𝑡) = 0, 𝑖 = 1,4. Let us solve the system (10) by numerical methods and present its solution graphically in Figure 2. Figure 2: Dynamics of E𝑘𝑖 (𝑡), 𝑖 = 1,4, when the arrival rate depends on the time 98 Figure 2 shows the dynamics of the expected number of customers in the network nodes 𝑆𝑖 , 𝑖 = 2𝜋𝑡 1,4, when the arrival rate is the function of time: 𝜆(𝑡) = 0.002sin ( + 2) + 0.0083. In this situation, 24 it is advisable to vary the number of agents depending on the arrival rate. In peak situations, node 𝑆3 has a queue. Figure 3 shows the dynamics of the expected number of customers in the queueing systems 𝑆𝑖 , 𝑖 = 1,4, when the arrival rate is constant: 𝜆 = 0.01. Figure 3 demonstrates that the process quickly reaches a steady state. The largest number of customers accumulates in node 𝑆2 . The servers in the nodes of the network are idle in the mean, and in this situation, the number of agents handling calls is sufficient for the absence of queues. Figure 3: Dynamics of E𝑘𝑖 (𝑡), 𝑖 = 1,4, when the arrival rate is constant Let us examine the expected monetary effect of a call center E𝐹 (𝑡), which is calculated using the formula (2). Let us calculate the value (2) under the assumption that 𝐹0 = 0, 𝐹1 = 𝐹2 = 20, 𝐹3 = 40, 𝐹4 = 10, 𝐿1 = 𝐿2 = 0,05, 𝐿3 = 0.07, 𝐿4 = 0.01 in monetary units. Figure 4 shows the dynamics of E𝐹 (𝑡) in the cases of the non-stationary arrival flow. Figure 4: Dynamics of E𝐹 (𝑡) when the arrival rate depends on the time The results of calculating the expected monetary effect of a call center for different sets of servers 𝑚 = (𝑚1 , 𝑚2 , 𝑚3 , 𝑚4 ) at the same time 𝑡 = 25 and the stationary arrival flow with 𝜆 = 0.01 are presented in Table 1. As is seen from the calculation results, the most optimal set is m = (3, 4, 3, 2), in which there are no queues and the monetary efficiency is maximal. 99 Table 1 The average monetary effect of a call center at time t=25 No. 𝑚1 𝑚2 𝑚3 𝑚4 E𝐹 (25) 1 3 4 3 2 251.59 2 4 4 3 2 251.54 3 3 3 3 2 229.97 4 2 4 3 2 234.67 5 3 4 2 2 207.09 Thus, the calculation results can be useful in analyzing and making decisions regarding the operation of the call center with different parameters. 5. Conclusion In this paper, the network stochastic model of a call center was presented as a queueing network with impatient customers. The model was investigated in the asymptotic case of a large number of customers. This allowed us to calculate the expected number of customers at each queueing node and determine the optimal staffing list of agents to maximize call center performance. The method is applicable both in transient and steady-state. Call center performance indicators can be found using mathematical methods for calculating the nodal characteristics of queueing networks. The areas of implementation of these research results are the design of call centers, the improving their performance, and the using queueing networks as models. In further studies, it is planned to use HM-networks, this type of queueing network first proposed by M. Matalytski [19, 20]. HM-networks assume revenues from changes in network states. 6. References [1] N. Gans, G. Koole, A. Mandelbaum, Telephone call centers: tutorial, review, and research prospects, Manufacturing & service operations management 5(2) (2003) 79-141. DOI: 10.1287/msom.5.2.79.16071. [2] L. Bennington, J. Cummane, P. Conn, Customer satisfaction and call centers: an Australian study, International journal of service industry management 11(2) (2000) 162-173. DOI: 10.1108/09564230010323723. [3] R. Stolletz, Performance Analysis and Optimization of Inbound Call Centers, Springer, Berlin, Heidelberg, 2003. DOI:10.1007/978-3-642-55506-0. [4] R. Srinivasan, J. Talim, J. Wang, Performance analysis of a call center with interactive voice response units, Top 12 (2004) 91–110. DOI:10.1007/BF02578926. [5] G. Koole, A. Mandelbaum, Queueing models of call centers: an introduction, Annals of operations research 113 (2002) 41–59. DOI:10.1023/A:1020949626017. [6] S. Aguir, F. Karaesmen, O. Akşin, F. Chauvet, The impact of retrials on call center performance, OR Spectrum 26 (2004) 353–376. DOI:10.1007/s00291-004-0165-7. [7] S. Zeltyn, A. Mandelbaum, Call centers with impatient customers: many-server asymptotics of the M/M/n + G queue, Queueing Systems 51 (2005) 361–402. DOI:10.1007/s11134-005-3699-8. [8] C. Kim, O. Dudina, A. Dudin, S. Dudin, Queueing system MAP/M/N as a model of the call center with the call-back option, in: K. Al-Begain, D. Fiems, JM. Vincent (Eds.), Proceedings of 19th international conference on Analytical and Stochastic Modeling Techniques and Applications, ASMTA 2012, Lecture Notes in Computer Science, LNCS, Springer, Berlin, Heidelberg, 2012, pp. 1–15. DOI:10.1007/978-3-642-30782-9_1. [9] S. V. Pustova, Investigation of call centers as retrial queueing systems, Cybernetics and System Analysis 46 (2010) 494–499. DOI:10.1007/s10559-010-9224-z. [10] F. Iravani, B. Balcıog̃lu, On priority queues with impatient customers, Queueing Systems 58 (2008) 239. DOI:10.1007/s11134-008-9069-6. 100 [11] G. A. Medvedev, Closed queueing systems and their optimization, Proceedings of the USSR Academy of Sciences. Engineering Cybernetics 6 (1978) 199–203. [12] M. Matalytski, T. Rusilko, A. Pankov, Asymptotic analysis of the closed queueing structure with time-dependent service parameters and single-type messages, Journal of applied mathematics and computational mechanics 12(2) (2013) 73–80. DOI:10.17512/jamcm.2013.2.09. [13] T. V. Rusilko, D. Ya. Kopats, Differential equations for the moments of the state vector of a closed queueing network with impatient customers, Herald of Dagestan State University. Series 1. Natural sciences 36(2) (2021) 20–30. DOI:10.21779/2542-0321-2021-36-2-20-30. [14] V. I. Tikhonov, M. A. Mironov, Markov Processes, Soviet Radio, Moscow, 1977. [15] K. V. Gardiner, Stochastic Methods in Natural Sciences, Mir, Moscow, 1986. [16] M. A. Matalytski, T. V. Rusilko, Approximate Methods for Analysis of Networks with a Central Queueing System and Their Applications, GrSU, Grodno, 2003. [17] Y. I. Paraev, Introduction to Statistical Dynamics of Management and Filtering, Soviet Radio, Moscow, 1976. [18] T. V. Rusilko, The first two orders moments determination method for the state vector of the queueing network in the asymptotic case, Vesnik of Yanka Kupala State University of Grodno. Series 2 11(2) (2021) 152–161. [19] M. A. Matalytski, Forecasting anticipated incomes in the Markov networks with positive and negative customers, Automation and remote control 78(5) (2017) 815–825. DOI:10.1134/S0005117917050046. [20] M. A. Matalytski, Analysis and forecasting of expected incomes in Markov networks with unreliable servicing systems, Automation and remote control 76 (2015) 2179–2189. DOI:10.1134/S0005117915120073. 101