Mathematical Model of Call Center in the Form of Multi-Server Queueing System

: The paper considers the model of a call center in the form of a multi-server queueing system with Poisson arrivals and an unlimited waiting area. In the model under consideration, incoming calls do not differ in terms of service conditions, requested service, and interarrival periods. It is assumed that an incoming call can use any free server and they are all identical in terms of capabilities and quality. The goal problem is to ﬁnd the stationary distribution of the number of calls in the system for an arbitrary recurrent service. This will allow us to evaluate the performance measures of such systems and solve various optimization problems for them. Considering models with non-exponential service times provides solutions for a wide class of mathematical models, making the results more adequate for real call centers. The solution is based on the approximation of the given distribution function of the service time by the hyperexponential distribution function. Therefore, ﬁrst, the problem of studying a system with hyperexponential service is solved using the matrix-geometric method. Further, on the basis of this result, an approximation of the stationary distribution of the number of calls in a multi-server system with an arbitrary distribution function of the service time is constructed. Various issues in the application of this approximation are considered, and its accuracy is analyzed based on comparison with the known analytical result for a particular case, as well as with the results of the simulation.


Introduction
One of the fastest-growing segments of the telecommunications market is the use of call centers. In this field, most of the costs of maintaining a call center are the salaries of employees who provide various kinds of information services to clients. The number of operators is chosen based on some standard measures, for example, the proportion of customers for which the waiting time for the start of service does not exceed the average time of the request servicing. The solution to the problem depends on the cost relationships between the components of the system. Therefore, the development of mathematical models and methods in this area is an necessary task. Traditional approaches of such studies lay in the field of queueing theory and related topics.
An excellent description of call centers and general issues in applying queueing models to studies were presented in [1]; some other general issues on the topic may be found in [2][3][4][5]. In this paper, we consider a call center model in the form of a multi-server queue with an incoming Poisson flow of calls, recurrent service, and an unlimited number of waiting places (M/G/N). A number of articles by different authors are devoted to the study of models of this type, for example: [6][7][8][9][10]. Furthermore, there are studies of models with non-Poisson incoming flows, for example: [11][12][13][14][15][16][17][18]. The aforementioned works make it possible to obtain characteristics of the operation of the system mainly using numerical, approximate, or asymptotic methods. Moreover, most of the papers only evaluate mean values of waiting time or queue length, or the distribution of waiting time, or the distribution of queue length, but only under asymptotic conditions. Our paper proposes one more approximate approach to solving the problem of studying the M/G/N model, which allows one to obtain an accurate enough approximation of the probability distribution of the number of calls in the system for a wide class of service time distributions. Unlike other publications, our approach allows one to obtain a distribution of the number of calls, and this is performed directly without any asymptotic conditions, as in [8,16]. In addition, the proposed approach is based on substituting an approximation before performing the analysis instead of adapting known results for M/M/N queues for the more general cases as was made in [9,10,12,14]. This allows us to obtain the final result in a simpler way with very high precision.
The main idea of the approach is in the approximation of the service time distribution by a two-phase hyperexponential one. The idea is taken from [19], where the problem of waiting time distribution was considered. We apply the idea to find the stationary distribution of the number of calls in this paper. For the model with hyperexponential service time, we use the matrix-geometric method to find the distribution of the number of calls in the system. To estimate the parameters of the hyperexponential distribution approximating a given distribution of the service times, in most cases, we use an approach based on three raw moments, which is in a certain sense similar to the results of [20]. However, for some cases, we have to use an approach based on two moments, and the results show us very high accuracy, even for such an approximation.
The only known exact analytical result is represented by the Pollaczek-Khinchine formula for a system M/G/1 with a single server [21,22]. In this paper, we compare the obtained approximation with this formula for the particular case N = 1 in the M/G/N model. For the other cases (when the number of servers N > 1), we use the simulation results for the comparison. The performed numerical experiments prove the high accuracy of the approach for models M/G/N in a wide class of service times.
The rest of the paper is organized as follows. In Section 2, we propose the mathematical model for the problem domain and formulate the goal of the study. In Section 3, we solve the problem for the partial case of the model with two-phase hyperexponential service times. To obtain the results for this case, we first derive the system of balance equation; after that, we apply the matrix-geometric method to solve it and derive the goal stationary probabilities at the end of the section. In Section 4, we propose a technique to obtain the results for the model M/G/N; also, we consider methods of the approximation and various issues of their applications. The results of estimating the accuracy of the proposed approach are presented in Section 5.

Mathematical Model
Consider a call center where we have several service devices (operators, servers) ( Figure 1). An incoming call is redirected instantly to any free server or should be placed in a waiting queue if all the servers are busy. Calls in the queue are ordered according to the FCFS rule (First Come First Serve). After the completion of the service at any server, the server takes the first call from the queue for its service. Consider a model of the system in the following form. Let calls arrive at the system according to the stationary Poisson process with intensity λ, let there be an unlimited number of places for waiting calls, and also, let there be N operators (servers) working independently from each other with service times distributed with distribution function B(x). For the goals of the study, we suppose that at least two finite raw moments exist for this distribution (the reasons may be found in Section 4).
Denote the stationary probabilities that we have exactly i calls in the system (both under the service and in the queue) by π(i), i ∈ {0, 1, 2, . . . }. The goal of the study is to find this probability distribution. The condition of existence of the stationary regime is obvious and standard for such a class of models: where µ is a rate of service in one server (the value that is inverse to the mean of a random variable with distribution function B(x)).
The approach that we use is based on the approximation of the distribution B(x) by the hyperexponential distribution with two phases. So, we first try to obtain the solution for the model with hyperexponential service times. After that, we propose methods for the approximation of an arbitrary probability distribution by the hyperexponential one. At the end of the paper, we estimate errors of the solution of the original problem on the basis of comparison with known analytical results or simulation results.
The choice of hyperexponential distribution for the goals of the study is made based on the fact that this distribution allows studies to be carried out on queueing systems where it is used, and it can also approximate many other types of random variables distributions with a quite good accuracy.

Solution of the Problem for the Case of Hyperexponential Service Time
Let the service times be hyperexponentially distributed with the CDF: where q 1 = q, q 2 = 1 − q, 0 < q < 1. Actually, this hyperexponential distribution is a mixture of two exponential distributions with parameters µ 1 and µ 2 chosen with probabilities q 1 and q 2 , respectively. So, actually, a server performs the service of a single call according to the exponential distribution with a chosen parameter. We name the servers that are currently working with parameter µ 1 as working at the first phase, and the servers that are currently working with parameter µ 2 as working at the second phase of the hyperexponential distribution. Denote the number of servers that are working at the first phase in time moment t by n 1 (t) and servers working at the second phase by n 2 (t). At any moment, the number of free servers equals N − n 1 (t) − n 2 (t). Furthermore, denote the number of calls in the queue in moment t by i(t). Notice that it may be greater than zero only if n 1 (t) + n 2 (t) = N.
Let us introduce the following notations for the stationary probabilities: when not all servers are busy), • P n (i) = Pr{n 1 (t) = n, n 2 (t) = N − n, i(t) = i} when all servers are busy. Here, 0 ≤ n ≤ N and i = 0, 1, 2, . . . So, the goal of the study for this model is to find these probabilities.

Balance Equations
We can construct the following system of the balance equations for the stationary regime. When the queue is empty, we have −P(n 1 , n 2 )(λ + n 1 µ 1 + n 2 µ 2 ) + P(n 1 − 1, n 2 )λq 1 + P(n 1 , n 2 − 1)λq 2 + P(n 1 , n 2 + 1)(n 2 + 1)µ 2 + P(n 1 + 1, n 2 )(n 1 + 1)µ 1 = 0 (3) for 1 ≤ n 1 + n 2 ≤ N − 1; for the non-empty queue, we can write We propose an algorithm for solving this system of equations with two stages. The first stage is based on the matrix-geometric method [23,24]. In this case, solution P n (i) of part (5) of the system for i > 0 will be expressed using term P n (0). This fact allows us to rewrite the rest of the system (2)-(4) in a closed form of the system of (N + 1)(N + 2)/2 equations with the same number of variables. The second stage of the proposed algorithm will solve this system.

The First Stage: Matrix-Geometric Method
Consider part (5) of system (2)-(5) for i > 0. Denote the row vector and vector P = P(0) to separate it from other components of the solution. Then, we can rewrite (5) in the matrix form: where matrices D and A have the form We write solution P(i) of system (6) in the matrix-geometric form: Substituting it into (6), we derive a quadratic equation for the matrix R: This equation can be solved by the iterative method starting from zero matrix O: until absolute values of the entries of matrix R 2 n A − R n D + λI become small enough. To find vector P, which is important for (7), we consider vector P(1) = P · R, whose entries can be written in the form Substituting these expressions into system (2)-(5), we obtain (9) we rewrite the system in the following form: This is a system of (N + 1)(N + 2)/2 linear equations with the same number of variables P(0, 0), P(n 1 , n 2 ) for 1 ≤ n 1 + n 2 ≤ N − 1, and P n (0). Due to the fact that this system is a homogeneous, its non-trivial solution can be determined in a form with a constant multiplier which can be obtained from the normalization condition. In the next section, we propose an algorithm that reduces system (10)-(12) consisting of (N + 1) (N + 2)/2 equations to the homogeneous system with N + 1 equations with variables P n (0), n = 0, 1, . . . , N that are entries of vector P, which we wish to derive.
Using these notations, we rewrite Equation (10) in the matrix form: where We can write two equations from (11), where n 1 + n 2 = 1 in the matrix form: where Further, for the case n 1 + n 2 = n, we can write corresponding equations from (11) in the form where matrices B n−1 with size n × (n + 1), D n with size (n + 1) × (n + 1) and A n+1 with size (n + 2) × (n + 1) have the following structures: Notation 0 in (14) and (15) means vector of zeros with corresponding sizes. Finally, we write the last N + 1 equations of system (10)-(12) in the form where S is a matrix with entries S kn from (9). So, gathering all matrix Equations (13)-(16), we can rewrite system (10)- (12) in the matrix form: From the first equation of the system, we can express: From the second equation, we express: From the next equations of system (17) for n ≤ N − 1, we obtain P n D n − P n−1 B n−1 = P n+1 A n+1 , P n D n − P n G n B n−1 = P n+1 A n+1 , Consider equality (18) for the case n = N − 1. We can write: because P N = P. Let us substitute expression (19) into the last equation of system (17).
Then we obtain the following equation for vector P: which is the homogeneous linear system for entries P n (0) of vector P. This system can be solved without any problems. So, we obtain vector P.

Stationary Probability Distribution of the Number of Calls in the System
As it is mentioned in Section 3.2, we obtained all vectors P n in the previous section, just in a form with a constant multiplier. To find true values of the vectors, we should apply the normalization condition where e n are column vectors consisting of ones. Denote stationary probabilities: Here, r(n) means the probability that n servers are busy, and p(i) means that all N servers are busy and there are exactly i calls waiting in the queue.

Approximation for the Case of Arbitrary Distributed Service Times
Consider the system shown in Figure 1 with i.i.d. service times with arbitrary distribution function B(x). To approximate probabilities π(i) from (25), we use expressions (23)-(24), supposing that distribution B(x) can be approximated by the hyperexponential distribution in form (1). So, the problem is reduced to establishing parameters µ 1 , µ 2 , q 1 , and q 2 = 1 − q 1 that provide enough accuracy for the approximation of distribution B(x) by the hyperexponential distribution (1).
It is not hard to obtain the expressions for the parameters of the two-phase hyperexponential approximation if we have two or three finite raw moments of the original distribution B(x). These expressions are presented in the following subsections.

Estimating Parameters of Hyperexponential Approximation Basing on Two Moments
Suppose distribution B(x) has two finite raw moments that we denote as b 1 and b 2 , respectively. We may estimate parameters µ 1 , µ 2 , q 1 , and q 2 from the following equalities for the moments: The following expressions may be derived by an easy method: (27)

Estimating Parameters of Hyperexponential Approximation Basing on Three Moments
If we suppose that distribution B(x) has three finite raw moments that we denote as b 1 , b 2 , and b 3 , respectively, then the parameters of the hyperexponential distribution can be derived as the solution of the following system of equations: Solving system (28) is not hard. If we denote we obtain It is clear that the approximation based on three moments is more accurate than one based on two moments. Therefore, it should be preferable if it is possible.

How to Apply the Approximations
Actually, estimation of the parameters in such a way may lead us to the following several types of results, including some unusual cases: Pairs µ 1 , µ 2 and q 1 , q 2 are complex numbers with their conjugates where all real parts are positive. 4.
One of parameters µ 1 or µ 2 has a negative value of its real part.
Obviously, in the cases 2 and 3, such values cannot be values of the parameters of hyperexponential distribution (1), but in the next section, we show that these values may be applied for approximation (25) and the result will have enough good accuracy. Moreover, when we deal with the third case with complex values of the parameters, corresponding function (1) is not a distribution function, but approximation (25) still works.
If you meet case 4, it is better to use estimation scheme (27) using two moments, b 1 and b 2 . This excludes the possibility of obtaining negative values for parameters µ 1 and µ 2 .

Accuracy of the Obtained Approximation
The analytical result for the single-server queue M/G/1 is well known and it is expressed by the Pollaczek-Khinchine formula. Unfortunately, we cannot compare obtained approximation (25) with the Pollaczek-Khinchine formula analytically, but we can compare the results numerically by performing evaluations on the basis of analytical expressions. For the case of a multi-server queue, there are no exact analytical results in the literature, so we compare obtained approximation (25) with the results of the simulation.
Consider single-server queueing system M/G/1. According to the Pollaczek-Khinchine formula, the characteristic function of the number of calls in the system has the form Here, is the Laplace-Stieltjes transform of distribution function B(x) of service time. Applying the inverse Fourier transform we obtain probability distribution Π(i), i = 0, 1, 2, . . . of the number of calls in the system. We determine the accuracy of approximation (25) by using the Kolmogorov distance for discrete distributions For the numerical experiment, we take the following three types of distribution functions B(x) with mean b 1 = 1: 1.
Gamma distribution with shape parameter α and inverse scale parameter β. The raw moments of the distribution are evaluated as follows:

2.
Weibull distribution with parameters γ and β = 1 , and with the raw moments

3.
Lognormal distribution with parameters σ 2 and a = − σ 2 2 . Its raw moments are the following: Consider single-server systems with the Poisson arrival process with intensity λ = 0.8 and service times with distribution functions B(x) from the list above. The results of estimation parameters of hyperexponential approximation (1) using expressions (29) and (30) for these distribution functions and corresponding values of Kolmogorov distance (33) (line ∆ 1 ) for the resulting approximation (25) are presented in Tables 1-3. Here, we compare distribution (25) with analytical results obtained using the Pollaczek-Khinchine formula (31) and (32). For the lognormal distribution in interval 0.4 < σ 2 < 0.7, applying formulas (29) and (30) leads to the case µ 2 < 0; therefore, we use estimation expressions (28) for this interval. The corresponding results are presented in Table 4. Furthermore, we cannot obtain the result in the case of gamma distribution with shape parameter α = 2. We think that this case can be resolved directly, but it requires additional study.    In line ∆ 5 of the tables, you may find the Kolmogorov distance for distribution (25) evaluated on the basis of hyperexponential approximations of service time distribution for the system with five servers (the intensity of the Poisson arrivals is taken equal to 4). The results are compared with the simulation results (note that the simulation results in the experiments are obtained with an error of 0.002 in terms of the Kolmogorov distance). All numerical evaluations for the examples are made using MathCAD software.
In Figure 2, you may find probability mass functions (PMFs) for some of the considered examples for the multi-server model with five servers. We chose the most interesting and the most "strange" cases of the approximation to use as examples. As we see from the presented results, distribution (25), using hyperexponential approximation of the service time, has very high accuracy in almost all cases (the worst results are only obtained for gamma distribution with a small value of the shape parameter). You may find various values of the parameters of the hyperexponential approximation in Tables 1-4, including strange ones (complex values, q > 1), but the final result (25) works properly for these cases too.

Discussion
The study presented in the paper shows that a two-phase hyperexponential approximation can be applied for service time distribution in the M/G/N queue. This approximation provides highly accurate results for the goal probability distribution of the number of calls in the system even in cases where the approximation cannot be used as a distribution function. A comparison of the results with analytical results for the case of a single-server queue and simulation results for the case of a multi-server queue show the high accuracy of the approach. Furthermore, in the paper, the analytical result for the system with hyperexponential service time is derived on the basis of the matrix-geometric approach. The result can be used both directly for the analysis of the M/H 2 /N queues or together with the hyperexponential approximation of the service time distribution for the analysis of the M/G/N queues.
We believe that hyperexponential approximation can be applied in studies that investigate other types of queueing models where it can allow one to obtain a solution of a problem or, at least, obtain it more easily than by other methods. For example, we see good possibilities in analyzing systems of G/G/1 and G/G/N types using the approximation and approach proposed in the paper.