=Paper= {{Paper |id=Vol-3057/paper9.pdf |storemode=property |title=Network Stochastic Call Center Model |pdfUrl=https://ceur-ws.org/Vol-3057/paper9.pdf |volume=Vol-3057 |authors=Tatiana V. Rusilko }} ==Network Stochastic Call Center Model== https://ceur-ws.org/Vol-3057/paper9.pdf
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