Determining Trunk Lines in Call Centers with Nonstationary Arrivals and Lognormal Service Times

Two important resources in a call center are the number of staff and the number of trunk lines required. In this paper, we focus on the decision of the number of trunk lines that a call center should have. The current practice is to use the Erlang B or the M/M/s/0 queueing model which assumes Poisson arrivals, exponential service times, s servers and no places in queue, i.e. no customers can wait. In this paper, we improve on the state of practice in determining the required number of trunk lines, by including two realistic features present in call centers. The first realistic feature is to consider nonstationarity of arrivals. The second feature is to consider the lognormal service time distribution instead of the exponential distribution. There is extensive empirical evidence for both features. In order to carry out our computations we use the results of a paper by Massey and Whitt, Operations Research, 44(6), 1996. We have two main findings. Firstly, we find numerically that in our nonstationary Erlang loss model, Mt/G/s/0, an insensitivity result holds. The blocking probability of arrivals at the call center depends only on the mean of the lognormal service time distribution and not on its variance. Our second finding is that current practice is quite robust. In particular, we find the number of trunk lines required using a stationary Poisson approximation. This approximation assumes stationary Poisson arrivals with an appropriately chosen arrival rate and exponential service times. The approximation does quite well in predicting the number of trunk lines required.


Introduction
Two important resources in a call center are the number of staff and the number of trunk lines required to connect calls to the call center. The first resource is certainly very important since it accounts for a major portion of the costs incurred in operating a call center. However the second resource, the number of trunk lines, has also got to be considered. The number of trunks in place determines the routing of calls to the staff and therefore has an important part to play in the realization of service level measures.
There are quite a few service level measures. One is ASA (Average Speed of Answer), i.e. how many seconds does a call have to wait on average before being answered. This includes those calls that do not have to wait at all. A second service level measure is, Delay of Delayed Calls. As mentioned previously, some calls do not have to wait at all. For those that have to wait, what is the average wait? A third measure is Service Level, what x% of calls are answered within a fixed time interval, y seconds. A common Service Level measure in call centers is to have 80% of calls answered within 20 seconds.
Usually, the cost of a trunk line is less than the cost of a single staff member or Customer Service Representative (CSR). If there are too few trunk lines, then of course calls get blocked from entering the call center. But if there are too many trunk lines relative to agents then that also can be a problem. Calls simply get put on hold and customer wait times increase significantly, resulting in bad service levels. Also many times the trunk lines are toll-free lines, so it is the call center instead of the customer, that is paying for the cost of waiting. This can also significantly increase costs for the call center.
Current practice in determining the number of trunk lines in a call center, is to first set a service level target. The service level target here is what percentage of calls can be allowed to be blocked. A typical service level target as mentioned in Reynolds [1] is to allow 2% blocked calls.
Then the Erlang B or the M/M/s/0 queueing model is used. This model implies Poisson arrivals, exponential service times, s servers and no places in queue, i.e., no customers can wait. In this model, servers refer to trunk lines and do not have anything to do with CSRs (or Customer Service Representatives). A formula for the blocking probability in the Erlang B model is used, see for example, Gross et al. [2]. The blocking probability from this formula is equated to the required service level (e.g. 2% blocked calls). Given this, the value of s, or the number of trunk lines is found.
For a more extensive discussion of call center trunking requirements, refer Reynolds [1], which is a good reference on call center workforce management.
In this paper, the state of practice is improved, in determining the required number of trunk lines, by including two realistic features present in call centers. The first realistic feature is to consider nonstationarity of arrivals. It has been repeatedly found in call center data that the number of arrivals vary quite significantly over the period of time the call center operates. One such small data set from Reynolds [1], is considered in our analysis here. The second realistic feature is to consider the lognormal distribution as the distribution of service times in the call center, instead of the exponential distribution. Empirical analysis of service time data in call centers has shown that service times follow the lognormal distribution. This evidence is discussed below. So the usual assumption of exponential service times may not be a good one. We next discuss more on these two realistic features.
Nonstationarity of arrivals is highly prevalent in call centers. Green, Kolesar and Whitt [3] plot hourly arrival rates for a financial services call center. There is significant variation in arrivals by time of day. Reynolds [1], reiterates this by mentioning, `The most accurate approach for call center forecasting involves time series analysis, which takes into account both trend and seasonality. It is the approach used in most call centers and serves as the basis for most of the automated workforce management forecasting models'. Brown et al [4] perform extensive statistical analysis of call center data. Their data comprises a complete operational history of a small Israeli banking call center, call by call, over a full year. Their plot of arrivals in calls/hour by time of day shows clear nonstationarity of arrivals.
Next evidence for the lognormal service time distribution is discussed. Firstly, a service time refers to the time spent by a CSR in talking to the customer, the time spent `on hold' while the CSR is processing the customer's request and lastly the time spent after the customer hangs up but while the CSR is still doing work related to the customer's request. Gans et al [5] provide extensive evidence of the lognormal distribution as the distribution of service times in a call center. Confirmations of the lognormal fit are provided in Bolotin [6], Chlebus [7], and Mandelbaum et al [8]. The authors also mention that the lognormal distribution has been found in unpublished call center data of a Dutch bank.
We now quote from Brown et al [4]. `Looking at the figure we see that the distribution of service times is clearly not exponential, as is assumed by standard queueing theory. In fact, after separating the calls with very short service times, our analysis reveals a remarkable fit to the lognormal distribution'. For all these reasons, the lognormal distribution as the service time distribution has been chosen.
The queueing model we consider in the paper, is the nonstationary Erlang loss model, M t /G/s/0, which has nonstationary Poisson arrivals, s servers (i.e. trunk lines) in parallel, no extra waiting spaces and i.i.d. service times which follow the lognormal distribution. Massey and Whitt [9] have analyzed this model for a general service time distribution. Their results have been used in this paper to do computations for determining the number of trunk lines required.
Massey and Whitt [9] in their paper, approximate a queueing model with a nonstationary arrival process with a queueing model with a stationary arrival process. The authors consider a fixed time interval and divide it into subintervals. They consider approximations for the blocking probabilities over subintervals by replacing the nonstationary arrival process over the subinterval by a stationary arrival process. The authors act as if the nonstationary M t arrival process were a stationary G arrival process and then try to approximate the stochastic variability. The approximation they propose is based on a heavy traffic peakedness formula. We discuss more on the results from their paper, as we use them in our calculations, in Section 3.
As mentioned before, in this paper we find the number of trunk lines required in call centers using nonstationary arrival data and the lognormal service time distribution. For our analysis, we consider two different lognormal distributions. Both have the same mean, but the second distribution has a variance double that of the first.
We have two main findings. Firstly, we find numerically that in our nonstationary Erlang loss model, an insensitivity result holds. The blocking probability of arrivals at the call center depends only on the mean of the lognormal service time distribution and not on its variance. In particular, both the lognormal service time distributions predict the same requirement of trunk lines.
Our second finding is that current practice is quite robust. In particular, we find the number of trunk lines required using a stationary Poisson approximation. This approximation assumes stationary Poisson arrivals with an appropriately chosen arrival rate and exponential service times (the Erlang B formula). The approximation predicts the same requirement of trunk lines as the original model. That this happens is not too surprising, given a statement made in this regard in the paper by Davis, Massey and Whitt [10]. However, it is still worthwhile to go through the numerical analysis and verify that the original model and the Poisson approximation match quite closely.
Kim and Park [11] consider 'two-stage' call centers where some incoming calls are completed by first service while others require an additional second service. The authors develop an effective outsourcing strategy in 'two-stage' call centers. To this end, they model a 'two-stage' service system and propose several call routing structures. The structures are compared through numerical testing.
Kilincli and Zhang [12] consider staff scheduling in call centers with cross-trained workers. Call centers face demand variations over time across multiple service categories and typically employ a cross-trained workforce with flexible schedules to hedge against these fluctuations. In practice, it is often impossible to cross-train agents in each category, thus partial and limited cross-training are the norm. To solve the problem an integer program that addresses cross-training, shift schedule, days off and break assignments across multiple service categories is proposed. The model is hard to solve and a two-phase sequential approach is developed. Experimental results with data from a call center with nine categories clearly demonstrate the significance of crosstraining. The authors find that partial limited cross-training, where 30% of staff is cross-trained with two skills or 10% of staff is cross-trained with three skills, could result in considerable cost saving. However, these savings could diminish quickly with the increase of efficiency loss in secondary skills.
Yu et al [13] consider delay announcements for call centers with hyperexponential patience modelling. Using real call center data, the patience distribution is modeled by the hyperexponential distribution. A state-dependent Markovian approximation is applied for computing abandonment. An empirical study shows good fit to real call center data. Simulation experiments indicate that the model and approximation are reasonable. The authors conclude that the approach could be applied to a voice response system of real call centers.
Li et al [14] state that efficient management of call centers needs accurate modeling of customer waiting behavior. This customer waiting behavior contains important information about customer patience (how long a customer is willing to wait) and service quality (how long a customer needs to wait to get served). The authors develop a two-way functional hazards model to study customer waiting behavior. The authors analyze data from a US Bank call center and provide insights about customer patience and service quality. The findings also provide inputs for call center agent staffing and scheduling.
Bimpikis and Markakis [15] motivated by applications such as call centers and healthcare, consider service systems that process two types of tasks that are unknown beforehand. There are also two kinds of servers with different skillsets. The service provider wants to maximize the systems long term throughput. Given this, what should be the resource allocation policy, i.e. how to assign tasks to severs over time? The authors show that the performance loss of the system due to uncertainty in task types can be significant. The authors also show that one optimal design could be a hierarchical structure. Tasks are initially routed to the least skilled servers and then progressively moved to more skilled servers.
The paper is organized as follows. In Section 2, the performance measure, i.e. call congestion is discussed. In Section 3, we discuss the results of Massey and Whitt [9] as used in the paper for our calculations. Section 4 describes the lognormal distribution while in Section 5 Boole's rule on numerical integration is presented. In Section 6, we present the data and the application of the method discussed in Section 3. Section 7 considers the stationary Poisson approximation and Section 8 presents Conclusions.

Performance Measure
The arrival process is a nonstationary Poisson process and is described by the deterministic arrival rate function, ( ) As in Massey and Whitt [9], we assume that the mean service time of a customer is 1. We also assume that T is not too small or too large. For a justification of this, please refer Massey and Whitt [9]. The authors recommend that T should be between 6 and 20. We take T =10. We assume that the average length of a call at the call center is 3 minutes and the planning interval under consideration is 30 minutes. Thus, if 3 minutes is taken as a service time of 1 time unit, we have T = 10 time units.
Let Q(t) be the number of busy servers at time t. The blocking probability is, which is the probability that all s servers (trunk lines) are busy at time t. The performance measure that is considered is the call congestion, c β . This is the ratio of the expected number of lost customers to the expected number of arrivals. As in Massey and Whitt [9], let B(t) be the number of blocked calls in the interval [0, t] and A(t) be the number of arrivals in the interval [0, t]. The ratio is, We next discuss the results of Massey and Whitt [9].

Results of Massey and Whitt [9] as Used in the Paper
As mentioned before, a queueing model with a nonstationary Poisson arrival process is approximated by a queueing model with a general stationary arrival process. In particular, the distribution of Q(t) in the nonstationary M t /G/s/0 model over the interval (0, T] is approximated by the distribution of Q(t) in the stationary G/G/s/0 model. Here, in our case, the distribution G is the lognormal distribution.
The first approach in Massey and Whitt [9] is to construct a stationary point process from the nonstationary Poisson arrival process. Let ≥ be the stationary point process constructed, please see Massey and Whitt [9], for details. The authors characterize the variability of the stationary point process using the index of dispersion for counts, The authors define 2 (1) c I = and find the value of 2 c for the special case of the linear arrival rate function, i.e. when ( ) We need this expression for 2 c in the approximation for the blocking probability that we use.
The second expression that we need for approximating the blocking probability is the expression for peakedness. According to Massey and Whitt [9], peakedness is defined as `the ratio of the variance to the mean of the steady-state number of busy servers in an associated infinite-server model with the same service time distribution and the same arrival process'.
The authors consider the limiting behaviour of the peakedness as the arrival rate grows. The heavy traffic peakedness is, (see Massey and Whitt [9]) with E[S]=1 and as mentioned before, G(.) being the cdf of the lognormal distribution. Let B(s, a) be the Erlang blocking formula for s servers and offered load a, extended to nonintegral s. Then, as the final result, we have the Hayward approximation. The call congestion is approximated by, Here we need to specify the expression for B(s, a). For nonintegral s, where ( ) is the incomplete gamma function, which is the standard Erlang loss formula for integral s, see Gross et al. [2]. We next discuss the lognormal distribution.

The Lognormal Distribution
For a discussion of the lognormal distribution, refer Wikipedia. The lognormal distribution is a continuous probability distribution. The lognormal random variable is such that its logarithm follows the normal distribution.
Define, Y Z µ σ = + , where Z has the standard normal distribution. Then Y has a normal distribution with mean µ and standard deviation σ . We have, Y X e = has the lognormal distribution. Its support is the positive real line. The pdf is, The mean is, and the variance is, The cdf is 1 1 ln ( ) erf 2 2 2 [ ] In the above erf(.) or the error function is a special function of sigmoid shape.
We have, For non-negative values of x, the error function has the following interpretation, as given in Wikipedia. For a random variable X that is normally distributed with mean 0 and variance ½, erf(x) is the probability of X lying in the range [-x, x].
Next numerical integration is discussed.

Numerical Integration: Boole's Rule
In the expression for the peakedness z, the following integral needs to be evaluated numerically, is the cdf of the lognormal distribution. This is done using Boole's rule, which is a quadrature formula. Before stating the rule, we state some definitions. For a discussion of Numerical Integration, refer Mathews and Fink [16].
Next the data and the application of the method above is discussed.

Data and Application of the Method
Reynolds [1] has the following data (Chapter 3, page 35), which are samples from a call center that takes calls from 8:00 AM to 6:00 PM daily. The data represents half-hourly call volumes for the previous three Mondays.  From Figure 1, it is seen that there is significant variation in call volumes by time of day. The maximum number of calls in a half-hour are 405 from 10:30-11:00, while the minimum calls during a half-hour are 100 from 5:30-6:00. This variation indicates a significant amount of nonstationarity in the arrival process.
It is assumed that the 10 hour interval is divided into 20 half-hour periods and that ( )  As regards the service time distribution, it is lognormal with mean 1. For our numerical work, took two lognormal distributions were taken, one with variance 1 and one with variance 2. The first lognormal distribution has a SCV (squared coefficient of variation, variance divided by the square of the mean) of 1, the same as an exponential, and so that is an important benchmark. The second lognormal distribution has a SCV of 2, much higher than the exponential. Since performance tends to degrade in a system with higher variance, we decided to use this distribution. In particular, in our application, degraded performance would imply that a higher number of trunk lines are required to achieve the same low blocking probability. We wanted to investigate if that is the case.
For a lognormal distribution with mean 1 and variance 1, it follows from equations (12) and (13) Although in the equation for peakedness, z, the upper limit of the integral is infinity, for practical purposes, it suffices to evaluate the integral upto an upper limit of 6 to get sufficient numerical accuracy. The upper limit of 6 in this case corresponds to Mean + 5*Standard Deviation.
Similarly, for the case of the lognormal distribution with variance 2, we have, It is again found that an upper limit of 8 gives sufficient numerical accuracy and corresponds to roughly Mean + 5.6*Standard Deviation, similar to the previous case. Now that 2 c and the two integrals have been evaluated, one each for each lognormal distribution, we can find the value of the peakedness z using equation (5) for each of the 20 time intervals and for each distribution. We are now ready to do the final calculation of the number of trunk lines required for each time interval and for each of the two distributions. Based on Reynold's [1], we set a blocking probability of 0.02. That is, number of trunk lines to make the blocking probability just below 0.02 is chosen. This is done as follows. We start with a low integral number of trunk lines, s, and keep increasing them till the blocking probability, ( ) / , / B s z z λ as given by equation (7), just falls below 0.02 for the first time. The value of s for which this happens is chosen as the number of trunk lines. This is done for each time interval and for each distribution.
Results are shown below. For each distribution, the number of trunk lines required is chosen as the maximum over all the 20 time intervals. Thus for both the lognormal distributions, 50 trunk lines are required to achieve a low 2% blocking probability.
If we look at Table 6 above, we find that the same number of trunk lines are required for each of the two service time distributions for each time interval. This is so barring one time interval. The two service time distributions have the same mean and differ significantly in their variance. For the standard M/M/s/0 Erlang loss model, it is known that an insensitivity result holds. That is, the blocking probability depends only on the mean of the service time distribution and not on the second moment and higher moments.
This insensitivity result need not hold for the nonstationary Erlang loss model. Davis, Massey and Whitt [10] find that `the service time distribution beyond the mean can have a significant impact on the time-dependent blocking probability in the nonstationary Erlang loss model'.
The authors consider the nonstationary M t /PH/s/0 model with a phase type (PH) service time distribution consisting of two phases. In the paper, the authors consider five different service time distributions, all of them being special two phase PH distributions. These are an Erlang (E 2 ), an exponential (M) and three hyperexponential (H 2 ) distributions. All the five distributions have mean 1. The SCVs (squared coefficient of variation) are 0.5 for E 2 , 1 for M and 4 for the three H 2 distributions. The authors find that the peak time-dependent blocking probabilities for the different service time distributions can differ by as much as a factor of 3.5. So the service time distribution beyond the mean can affect the blocking probability in the nonstationary Erlang loss model.
However the above situation does not always have to hold for the nonstationary Erlang loss model. In our case, both the lognormal service time distributions have mean 1. The first distribution has variance 1, the second has variance 2. We numerically find that in our case, in the M t /G/s/0 model, an insensitivity result holds. The blocking probability depends only on the mean of the lognormal service time distribution and not on its variance.
Next the stationary Poisson approximation is discussed.

Stationary Poisson Approximation
In this approximation, the nonstationary Poisson process is replaced with a stationary Poisson process having arrival rate λ . Again the 10 hour time period is divided into 20 half hour intervals. The interval which has the highest average arrival rate is found. This is the arrival rate λ that we use in the Erlang B formula of the Erlang loss model, stated below. The service rate is 1. The Erlang B blocking formula is given below for the blocking probability ( ) The calculation of the Erlang B formula can cause numerical problems, because of the large values of s!, as s is reasonably large. It is discussed in Gross et al [2] that ( ) , B s λ can be computed using the following iterative method.
The above iterative method is used to compute ( ) , B s λ .
We substitute 40.25 λ = , which is the highest average arrival rate for time interval 7 (see Table 2 ). The value of s such that ( ) Therefore using the stationary Poisson approximation, the number of trunk lines required to achieve a 2% blocking probability is 50.
We can compare this result with the result in Table 7. We thus find in our numerical work that the stationary Poisson approximation is robust and works quite well.
In this paper, we have considered finding the number of trunk lines required using a nonstationary Poisson process and the lognormal service time distribution. The nonstationary Poisson arrival process has been observed in many call centers and the lognormal service time distribution is also justified by empirical work. Existing practice for determining the number of trunk lines assumes a stationary Poisson arrival process and exponential service times and uses the Erlang loss model. Our numerical work finds that existing practice is quite robust and works well in determining the number of trunk lines required, despite the limiting assumptions in the analysis.
This finding of the robustness of existing practice, also finds support in a statement made in Davis, Massey and Whitt [10]. According to the authors, a common engineering approach is to use the stationary Erlang loss model with a constant arrival rate during the time interval over which the system is most heavily loaded. With this approach, the assumed arrival rate in the model is usually greater than the real arrival rate the majority of the time. According to the authors, this approach has been very successful in designing systems with a fixed number of servers that must be able to satisfy demand at any time.
Given the above, it is not too surprising that the stationary Poisson approximation works so well, in this case.
The computations discussed in this section and the previous section, were carried out in Excel and MATLAB.
We next present Conclusions.

Conclusions
Two important resources in a call center are the number of staff and the number of trunk lines required. In this paper, we focus on the decision of the number of trunk lines to have. The current practice is to use the Erlang B or the M/M/s/0 queueing model which assumes Poisson arrivals, exponential service times, s servers and no places in queue, i.e. no customers can wait. In this paper, we improve on the state of practice in determining the required number of trunk lines, by including two realistic features present in call centers. There is extensive empirical evidence for both features as found in the papers by Gans et al [5] and Brown et al [4].
In order to carry out our computations we use the results of a paper by Massey and Whitt [9]. The authors approximate a queueing model with a nonstationary arrival process with a queueing model with a stationary arrival process. In particular, the distribution of the number of busy servers in the nonstationary M t /G/s/0 model is approximated by the distribution of the number of busy servers in the stationary G/G/s/0 model. Here, in our case, the distribution G is the lognormal distribution.
We have two main findings. Firstly, we find numerically that in our nonstationary Erlang loss model, M t /G/s/0, an insensitivity result holds. The blocking probability of arrivals at the call center depends only on the mean of the lognormal service time distribution and not on its variance. In particular, both the lognormal service time distributions predict the same requirement of trunk lines. In the stationary M/M/s/0 Erlang loss model, the insensitivity result theoretically holds, as is well known. Davis, Massey and Whitt [10] have shown that in the nonstationary Erlang loss model, the insensitivity result need not hold. We find numerically find that in our model, it holds.
Our second finding is that current practice is quite robust. In particular, we find the number of trunk lines required using a stationary Poisson approximation. This approximation assumes stationary Poisson arrivals with an appropriately chosen arrival rate and exponential service times. The approximation does quite well in predicting the number of trunk lines required. Based on a statement in Davis, Massey and Whitt [10], this finding is not too surprising. However, it is worthwhile to go through the numerical analysis and come to this conclusion.
Future work could consider determining call center staffing levels under these more realistic assumptions. We are presently working along those lines.