Keywords: call centers, stochastic models, queuing systems, Markov chains, simulation, scheduling policies, quality of service, transient analysis, skill-based routing, multi-channel call centers

Chapter 1

Introduction

In this chapter, we give a general introduction of the manuscript. First, we describe the context and the motivation of the thesis. Second, we highlight the approach, the structure and the main contributions of this work.

Context and Motivation

A call center is a service system. It is a facility designed to support the delivery of some interactive service via telephone communications, email, chat, etc. The definition of a call center is continuously changing with technological development, but the core fundamentals of a customer making a call (via a phone, email, web site, fax or Interactive Voice Response) to a center (collection of resources) will remain constant. This thesis focuses on operations management issues for multi-skill and multichannel call centers. In what follows, we first present the context of the call center industry and the related operations management issues. Second, we focus on the motivation of this work, and present our collaboration with the French consulting company Interact-iv, which was at the origin of most of the addressed problems in this thesis.

Context. Call centers make up a large and growing part of the global economy. They are very labor-intensive operations, employing millions of persons across the globe. Call centers serve as the public face in various areas and industries: insurance companies, emergency centers, banks, information centers, help-desks, tele-marketing, etc. The success of call centers is due to the technological advances in information and communications systems, see [START_REF] Pinedo | Call Centers in Financial Services: Strategies, Technologies, and Operations[END_REF]. The most important call centers equipments are the Interactive Voice Response (IVR), the Automated Call Distributor (ACD), and the Computer Telephone Integration (CTI). These technologies have grown cheaper, more reliable, and more sophisticated. Moreover, these advances enabled various call center tasks which requires multiple skills and channels.

In multi-skill call centers, the call assignment strategy of Skills-Based Routing (SBR), is used to assign incoming calls to the most suitable agent. The report of [START_REF] Holman | The Global Call Center Report: International Perspectives on Management and Employment[END_REF] made on 2500 call centers in 17 countries with 475,000 employees points out that 56% of call centers use SBR strategies. These strategies are an enhancement to ACD systems.

Next, the development of alternative channels goes together with an adaptation to impatient customers with higher expectations. The recent report of ICMI (2013), based on the analysis of 361 large contact centers, presents the increasing use of new channels and the related research issues.

In particular, they point out that outbound tasks require intensive integration with inbound ones in most call centers. Although the inbound calls remain present in most call centers (98%), emails are also widely used (89%). Moreover, outbound calls (76%), Web (70%) and chats (40%) are important and developing channels. We refer the reader to the description of the general context of multi-channel call centers management in Chapter 7 of [START_REF] Koole | Call Center Optimization[END_REF].

Due to the operational difficulties to find better solutions than intuitive ones, managers have a continuous interest in the related research disciplines. The literature on operations management in call centers has focused on the following issues: demand forecasting, quality of service, capacity planning, queueing, call routing, staffing and agents scheduling. The related main academical dis-ciplines are Mathematics and Statistics, Operations Research, Industrial Engineering, Information Technology, Human Resource Management, as well as Psychology and Sociology. The overall purpose of this literature is helping the manager in improving the management of their call centers.

We refer the reader to the complete surveys of the academic literature on call center operations management by Gans et al. (2003) and Akşin et al. (2007).

The goal of the present thesis is to contribute to the operations management research in multiskill and multi-channel call centers. The purpose is to enhance our understanding of such complex systems, so as we obtain useful guidelines for the practitioners.

Motivation.

In what follows, we want to motivate the problems under consideration. These are related to flexibility in multi-skill call centers and also to the routing issues in multi-channel call centers.

The concept of flexibility is related to the ability of a company to efficiently match its capacity to an uncertain demand with multiple types. A wide literature has focused on the distribution of skills par agent. Increasing the number of skills per agent goes together with a better use of the resources but also with more costly resources. A well known studied and efficient configuration is chaining, first pointed out by [START_REF] Jordan | Principles on the Benefits of Manufacturing Process Flexibility[END_REF] in the context of manufacturing systems.

In this configuration each agent has only two skills and the distribution of the skills corresponds to a chain. A precise definition of the chaining architecture will be given in Chapter 2. Developing intelligent configurations such as chaining is very interesting for practitioners. The value of these configurations is that they capture the benefits of pooling by only having a limited flexibility.

However, the robustness of chaining fails in the case of asymmetric parameters [START_REF] Sheikhzadeh | Machine Sharing in Manufacturing Systems: Total Flexibility versus Chaining[END_REF]). The situations with asymmetric parameters arise in practice. The typical example is that of an European multilingual call center where customers call from several countries. For instance in Bluelink (a service provider of Air France KLM), each agent speaks two languages: her own native language and English. The workload is unbalanced ranging from only some few calls from a given country to several thousand of calls from another country. For such cases in practice, it is important to develop new architectures that allow on the one hand to account for demand asymmetry, and on the other hand to capture the benefits of pooling with only a limited flexibility.

Call centers require a very accurate match of demand and supply. Since the volatility of call arrival patterns is high, there is often a mismatch between demand and the scheduled number of inbound agents. Moreover, even if the demand is accurately forecasted, a considerable overcapacity should be scheduled to be able to deal with the random Poisson fluctuations of the demand. To prevent idle overcapacity and to limit the necessity to have extremely accurate forecast, inbound calls are sometimes mixed with other types of channels which have a less strict allowable delay, such as emails or outbound calls. This is referred to as (call) blending. It arises in the context of multi-channel call centers. Next, we describe some motivational examples for blended operations issues.

In practice, we may find situations where a conversation between an agent and a customer contains a natural break. For example, an agent of an internet hotline asks the customer to reboot her modem or her computer which may take some time where no interactions can take place. It is also often the case that a call center agent of an electricity supplier company asks the customer for the serial number of her electricity meter box. Another example is that of commercial call centers with a financial transaction during the call conversation. Inside an underway conversation, the agent is then free to do another task if needed. For an efficient use of the agent time, there might be an opportunity to route the less urgent jobs (emails) to agents, not only when the system is empty of calls, but also during the call conversations. An interesting research question here is how should be the routing rules as a function of the system parameters.

Even in the classical case of a single stage call conversation, the ACD programming is still a complex task for a blending situation. [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF] and [START_REF] Gans | A call-routing problem with service-level constraints[END_REF] show that efficient assignment policies are those with agent reservation for inbound calls. The main complexity comes from the fluctuation of the system parameters, in particular those of the jobs arrival processes. For instance from the statistical analysis of a call center data provided by Interact-iv, we observe that during a given period the arrival volume can triple from one day to another. The reasons of the fluctuations is hard to determine and an observation can hardy be duplicated on a future period. Given this fact, one could focus on developing routing policies with a continuous adaptation of the agent reservation threshold, while using at a minimal level the forecasted system parameters.

In the context of highly congested call centers, the use of a callback option can be proposed to customers so as to balance workload and avoid excessive abandonments. Since a callback option transforms an inbound call into an outbound one, the issue in the management of this option is somewhat similar to that of a blended situation. Some practical specific problem can be pointed out: What should be the routing rules of the jobs in the ACD, in order to optimize the system performance in terms of the waiting times of inbound and outbound jobs? This work is done, in its major part, under a collaboration with the French consulting company Interact-iv. Interact-iv sells software, advice and methods to call centers. The customers of Interactiv are for a large part small multi-channel call centers. Through the collaboration, the purpose of Interact-iv is to provide to its customers (call center managers) solutions that are thoroughly supported quantitatively. We had the opportunity to work on various issues of multi-channel call centers and had access to real call center data. This collaboration offered a wealth of learning opportunities.

Structure and Main Contributions

In this section, we describe the structure and the main contributions of the manuscript. We briefly describe the different chapters separately and give their corresponding submitted or working papers.

The current thesis can be divided into two parts. The topic of the first part is the design of multi-skill call center architectures. It corresponds to Chapter 2. The topic of the second part is the optimal routing in multi-channel call centers. It corresponds to Chapters 3, 4 and 5.

In Chapter 2, we focus on architectures with limited flexibility for multi-skill call centers. The context is that of call centers with asymmetric parameters: unbalanced workload, different service requirements, a predominant customer type, unbalanced abandonments and high costs of crosstraining. The well known architectures with limited flexibility such as chaining fail against such asymmetry. We propose a new architecture referred to as single pooling with only two skills per agent. We provide a comparison framework between chaining and single pooling and demonstrate the efficiency of single pooling under various situations of asymmetry. We also develop analytical results for particular single pooling models, in order to get some sense on the effect of arrival asymmetry on performance. This Chapter is based on [START_REF] Legros | A Flexible Architecture for Call Centers with Skill-Based Routing[END_REF] (under second round revision in International Journal of Production Economics).

In the second part, we focus on routing problems in multi-channel call centers. In Chapter 3, we consider a blended call center with calls arriving over time and an infinitely backlogged queue of emails. The call service is characterized by three successive stages where the second one is a break.

We define parameters of control for the routing of emails between or inside calls treatment. Next, we develop a method based on the analysis of Markov chains in order to derive the performance measures of interest for calls and for emails. We focus on optimizing the email routing parameters.

In addition, we develop an approximation method for the system performance evaluation under the light-traffic regime. We also propose an approximation method to extend the results to the multi-server case. We derive various structural results and conclude that all the time at least one of the two email routing parameters has an extreme value. This chapter is based on [START_REF] Legros | Optimal Email routing in a Multi-Channel Call Center[END_REF] (submitted to Stochatic Systems).

In Chapter 4, we examine a threshold policy that reserves agents for inbound calls. We study a general non-stationary model where the call arrival follows a non-homogeneous Poisson process. The optimization problem consists of maximizing the throughput of outbound tasks under a constraint on the waiting time of inbound calls. We propose an efficient adaptive threshold policy easy to implement in the Automatic Call Distributor (ACD). This scheduling policy is evaluated through a comparison with the optimal performance measures found in the case of a constant stationary arrival rate, and also a comparison with other intuitive adaptive threshold policies in the general non-stationary case. This chapter is based on Legros et al. (2013a) (submitted to IIE Transactions).

In Chapter 5, we consider a call center model with a callback option, which allows to transform an inbound call into an outbound one. The optimization problem consists of minimizing the expected waiting time of the outbound calls while respecting a service level constraint on the inbound ones. We propose a routing policy with two thresholds, one on the reservation of the agents for inbound calls, and another on the number of waiting outbound calls. A curve relating the two thresholds is determined. This chapter is based on the ongoing paper [START_REF] Legros | Call centers with a call back option[END_REF].

In Chapter 6, we close the thesis by giving general concluding remarks and highlighting directions for future research.

Chapter 2 A Flexible Architecture for Call

Centers with Skill-Based Routing

We focus on architectures with limited flexibility for multi-skill call centers. The context is that of call centers with asymmetric parameters: unbalanced workload, different service requirements, a predominant customer type, unbalanced abandonments and high costs of cross-training. The most knowing architectures with limited flexibility such as chaining fail against such asymmetry. In this paper, we propose a new architecture referred to as single pooling with only two skills per agent and we demonstrate its efficiency. We conduct a comprehensive comparison between this novel architecture and chaining. As a function of the various system parameters, we delimit the regions where either chaining or single pooling is the best. Single pooling leads to a better performance than chaining while being less costly under various situations of asymmetry: asymmetry in the number of arrivals, in the service durations, in the variability of service times, or in the service level requirements. It is also shown that these observations are more apparent for situations with a large number of skills, or for those with a large call center size.

Introduction

Context and Motivation. The concept of flexibility is related to the ability of a company to efficiently match its capacity to an uncertain demand with multiple types. The need for flexibility arises in a wide range of manufacturing systems. It also extends to service systems, such as call centers, where different types of customers ask for a quasi-instantaneous processing. Resource flexibility in call centers reduces to cross-training agents, which allows to improve both the utilization and the performance. Since cross-training agents is achieved with higher operating costs, resource flexibility could result in a trade-off between performance and cost. The performance is measured through operational indicators such as the expected waiting time, the probability of waiting, and the waiting time distribution, or also through human resource aspects that result in a higher efficiency of the agents. Cross-training may improve the agent motivation and provides a career path.

In this chapter, we only focus on the operational indicators.

The process flexibility problem have been studied in different directions, such as machine sharing, multi-stage supply chains, queueing systems and flexible workforce scheduling. Here we consider flexibility questions in the context of queueing models for call centers. A wide literature has focused on contrasting two extreme situations. The full flexible architecture (FF) versus the the full dedicated (FD) one. In the FF model, each agent is fully cross-trained for all call types. In most situations in which call types have similar service duration requirements, FF would require less agents than any other architecture, in order to reach a given predefined service level. The reason is that it benefits from the economies of scale, which absorb stochastic variability [START_REF] Borst | Dimensioning Large Call Centers[END_REF]). However the agents in FF are too costly and even sometimes impossible to find.

As commented by [START_REF] Marengo | Skill based routing in multi-skill call center[END_REF], the multilingual Compaq call center certainly could not find or train agents to speak eleven languages! In the other extreme situation of the FD model, an agent is only trained to handle a single call type. Agents are then less costly, but FD would require a larger staffing level to reach the same service level as in FF or any other architecture.

Full flexibility and full dedication, however, are only two extreme situations. A well known and studied intermediate configuration is chaining, first pointed out by [START_REF] Jordan | Principles on the Benefits of Manufacturing Process Flexibility[END_REF]. In the chaining model, each call type can be assigned to one of two adjacent agent teams, and each agent can handle calls from two adjacent types. [START_REF] Sheikhzadeh | Machine Sharing in Manufacturing Systems: Total Flexibility versus Chaining[END_REF], [START_REF] Gurumurthi | Modeling and Analysis of Flexible Queueing Systems[END_REF], and [START_REF] Jordan | Chained Cross-Training of Workers for Robust Performance[END_REF] prove that chaining, with an appropriate linkage between demand and resource types, behaves just as well as full flexibility. In the context of Constant Work in Process (CONWIP) serial production lines, Hopp et al. (2004) showed that the impact of forming a complete chain of skill sets can be substantial in increasing throughput. [START_REF] Wallace | A Staffing Algorithm for Call Centers with Skill-Based Routing[END_REF] consider the problem of routing and staffing in multi-skill call centers. They again confirm the principal that a little flexibility has the potential to achieve the performance of total flexibility. Using simulation they demonstrate that the performance, with an appropriate and limited crosstraining of agents (two skills per agent) such as in chaining, is almost as good as when each agent has all skills.

Developing intelligent configurations such as chaining is very interesting for practitioners. The value of these configurations is that they capture the benefits of pooling by only having a limited flexibility. However, the robustness of chaining fails in the case of asymmetric demand [START_REF] Sheikhzadeh | Machine Sharing in Manufacturing Systems: Total Flexibility versus Chaining[END_REF]. By asymmetric demand, we mean different workload intensities and service time requirements, and also different variabilities in inter-arrival and service times. For such cases in practice, it is important to develop new architectures that allows from the one hand to account for demand asymmetry, and from the other hand to capture the benefits of pooling with only a limited flexibility.

In this chapter, we consider skill-based routing (SBR) call centers with two particular features: demand asymmetry and costly/difficult agent training. The typical example is that of an European multilingual call center where customers call from several countries. It is difficult for managers to find agents speaking more than two languages. For instance in Bluelink (the service provider of Air France KLM), each agent speaks two languages: her own native language and English. Note that Bluelink is more interested in agents speaking two languages than those speaking three or more languages. The reason is that the latter often feel themselves over-qualified. They are therefore likely to leave the company faster than the others, which increases the turnover. The workload is also unbalanced ranging from only some few calls from a given country to several thousand of calls from another country. Another example is that of post-sales service call centers of retailers such as Darty and Fnac which are French distributors of white goods, telecommunications products, information technology, but also internet services, photo services or travel services. We also give the example of retail banking call centers where questions are with regard to savings or stock exchange for examples. The main characteristics in the previous examples are (i) the demand is unbalanced, (ii) the nature of the required agent skills can be very different which make difficult or too costly the agent training, and (iii) one may find a predominant and "easy" type of questions that could be handled by most of the agents without any particular training, for example the English task in a multilingual call center, account information and simple bank tasks in banking, order tracking and payment for retailers, etc.

Main findings. Motivated by this prevalence in practice, we propose in this chapter a new call center architecture that can be used instead of chaining. For such cases, applying chaining is too costly and difficult to implement (many combinations of two tasks per agent are even hard to obtain). Moreover, existing literature have shown that chaining is not appropriate for such demand situations: unbalanced workload of the "difficult" tasks and a predominant "easy task". As proven in [START_REF] Bassamboo | Optimal Flexibility Configurations in Newsvendor Networks: Going Beyond Chaining and Pairing[END_REF], the tailored pairing architecture is efficient for small systems including those with asymmetries. However, this architecture requires an important number of cross-trained teams which might be again difficult to implement in practice. We propose a new organizational model, referred to as single pooling, where we dedicate a team of agents to each difficult type of calls, and the easy type of calls have access to all agents from all teams. Balancing the workload among the agents in this way captures the benefits of pooling without requiring every agent to process every call type.

A concise definition of our model will be given later. We do not claim that our model is better than chaining in all cases, but only in the particular range of parameters as shown in the call center examples above. The value of our architecture is that it has a low degree of flexibility (each agent handles one difficult type and the easy task) while behaving in terms of performance as a fully flexible call center. This is important in practice since additional flexibility often comes at the cost of high operating overhead. Hence, the results of our analysis have significant managerial implications.

Using simulation, we conduct a comprehensive comparison between this novel architecture and chaining. As a function of the various system parameters, we delimit the regions where either chaining or single pooling is the best. Our key findings are highlighted next. Single pooling leads to better performance while being less costly than chaining under various situations of asymmetry between the customer types: asymmetry in the number of arrivals, in the service duration, in the variability of service times, or in the service level requirements. Moreover, we conclude that these observations are more apparent for situations with a large number of skills, or for those with a large call center size. In practice, the issue of limiting the flexibility appears more in large call centers, rather than in small ones with a few number of agents. In small call centers, the number of customer types is often very limited or they are very similar in terms of the required agent skills, so that the agents are usually full-flexible. Hence, there is often no need for managers to deal with cross-training questions. These insights show that there might be opportunities for managers of call centers to improve performance using the single pooling architecture.

The rest of the chapter is organized as follows. In Section 2.2 we review some of the literature related to this chapter. In Section 2.3 we describe chaining and single pooling models, and provide the comparison framework. In Section 2.4, we develop analytical results for particular single pooling models, in order to get some sense on the effect of arrival asymmetry on performance. In Section 2.5, we use simulation to compare between the two call center models under various situations of asymmetry on the system parameters. Section 2.6 concludes the chapter and highlights some future research.

Literature Review

There is an extensive and growing literature on call centers. We refer the reader to Gans et al. (2003) and Akşin et al. (2007) for an overview. We review in what follows some of the literature related to this work.

Impact of Pooling. The value of pooling comes from the creation of flexibility. The general known intuition is that pooled systems are more effective than independent ones. The impact of pooling has been first studied in [START_REF] Smith | Resource Sharing for Efficiency in Traffic Systems[END_REF]. They show that pooling always leads to a better performance in terms of the expected delay in queue. [START_REF] Akşin | Characterizing the performance of process flexibility structures[END_REF] investigate the impact of the call center size on the opportunity to add flexibility. They demonstrate that a small call center will benefit more from adding flexibility than a large one. [START_REF] Benjaafar | Performance Bounds for the Effectiveness of Pooling in Multi-Processing Systems[END_REF] studies the impact of pooling for a variety of manufacturing, telecommunication and computer systems. He considers a multi-processing system consisting of several facilities and shows that in some situations of heterogeneity in the workloads, increasing flexibility can deteriorate performance. [START_REF] Mandelbaum | On Pooling in Queueing Networks[END_REF] consider stochastic service systems modeled as queueing networks. The service of a customer amounts to a collection of tasks. They show that adding flexibility does not automatically improve performance. They point out that adding a partial flexibility could be devastating for a queueing network. Recently, van Dijk and van Der Sluis (2008) show in the context of SBR call centers that without any clever routing rules and under a high variability in the call types and the resources, pooling could deteriorate the performance in terms of the average waiting time. In their work, [START_REF] Tekin | Pooling Strategies for Call Center Agent Cross-Training[END_REF] investigate the efficiency benefits achievable via cross-training in SBR call centers. They conclude that under first come, first served (FCFS), the pooling of the dedicated teams is appropriate for heavy workloaded teams. However, it is not necessarily the case for teams with light workloads. They then use the difference between the arrival rates of the customers as a choice parameter based on which the decision of pooling would be taken or not. Inspired by the results of [START_REF] Smith | Resource Sharing for Efficiency in Traffic Systems[END_REF], they also conclude that pooling teams could be counterproductive if services time means are very different from one customer type to another (for example when one is six times higher than the other ones).

Flexible Architectures. The most fundamental work on flexibility is that by [START_REF] Jordan | Principles on the Benefits of Manufacturing Process Flexibility[END_REF] for the automobile assembly plants, but it can be also applied to broader manufacturing system settings. They conduct an extensive simulation study and conclude that "a little flexibility can achieve almost all the benefits of total flexibility" under a configuration referred to as chaining, with two product types per plant. They demonstrate that the expected shortfall and capacity utilization of chaining resources are close to those under a full flexible configuration. [START_REF] Garavelli | Performance analysis of a batch production system with limited flexibility[END_REF] considers the setting of job shop cellular manufacturing systems used to perform batch production. He finds similar results to those by [START_REF] Jordan | Principles on the Benefits of Manufacturing Process Flexibility[END_REF] through a comparison between the performance of a full dedicated system, a full flexible one and chaining. Similar results are found by [START_REF] Garavelli | Flexibility configurations for the supply chain management[END_REF] in a complex supply chain environment, requiring the coordination of many plants producing good to customers located in different places. Again in the context of cellular manufacturing systems, [START_REF] Albino | Limited flexibility in cellular manufacturing systems: A simulation study[END_REF] analyze the benefits of a limited flexibility.

Starting from two industrial case studies concerning in-house metalworking shops, Nomden and van der Zee (2008) find by simulation that a chained distribution of routes behave very well.

For queueing systems, [START_REF] Gurumurthi | Modeling and Analysis of Flexible Queueing Systems[END_REF] compare different scenarios of adding flexibility under different routing policies. They prove that the value of chaining decreases for an asymmetric demand. Hopp and van Oyen (2004) consider the question of how to cross-train a worker to two skills in the context of serial production lines. They conclude that a novel strategy called skill chaining strategy is more robust against variability than a cherry-picking strategy (a team is full flexible) when demand is symmetric. The cherry-picking strategy in a serial production line can be seen as similar to single pooling, where the customers are the machines, and the bottleneck machine represents the easy type of calls. [START_REF] Tomlin | On the Value of Mix Flexibility and Dual Sourcing in Unreliable Newsvendor Networks[END_REF] consider the context of unreliable supply chains that produce multiple products. They study four canonical supply chain design strategies, where one of them, referred to as dual-source flexible, has been already proposed by [START_REF] Chevalier | Routing and Staffing in Large Call Centers with Specialized and Fully Flexible Servers[END_REF] in the context of call centers. They refine the prevailing intuition that a flexible network is preferable to a dedicated network by proving that this intuition is valid if either the resource investments are perfectly reliable or the firm is risk neutral. In a similar setting to ours, [START_REF] Robbins | Cross Training in Call Centers with Uncertain Arrivals and Global Service Level Agreements[END_REF] introduce an SBR call center queueing model with two customer types, referred to as partial pooling. They consider two dedicated agent teams for each customer type, and one cross-trained team for both types. They show that cross-training a small number of agents can deliver a substantial benefit. They also find the level of cross-training that minimizes staffing costs, while satisfying a service level constraint. [START_REF] Bassamboo | Optimal Flexibility Configurations in Newsvendor Networks: Going Beyond Chaining and Pairing[END_REF] study the flexibility problem with a newsvendor network model of resource portfolio investment. They conduct a comparison between chaining and tailored pairing. They show that a system that combines dedicated and cross-trained agents is asymptotically optimal. They also show using simulation experiments that the "tailored pairing" design is superior for small systems, including systems with asymmetries.

The tailored pairing architecture might be one of the best propositions in the literature to deal with the asymmetric parameters. However as already mentioned above, we can not retain this architecture as a reference in our project. The reason is that for call centers with many skills, working under tailored pairing may lead to non-realizable situations. In such a case, the number of two-skills combinations could be very high. [START_REF] Garnett | An Introduction to Skills-Based Routing and its Operational Complexities[END_REF] argue on the importance of adapting the system architecture to the asymmetry in the customer arrival rates. In summary, chaining is robust according to its ability to support variability. It however fails when demand is asymmetric. It can be also too expensive to train agents on various combinations of two skills. For these situations, we propose and analyze in this project a new efficient configuration of a queueing call center model.

Agent Skills, Staffing and Routing. In an SBR call center, agents can often only be trained for a subset of skills. One key management issue is to determine the subset of skills that will be considered, and the number of agents for each subset of skills. [START_REF] Pinker | The Efficiency-Quality Trade-Off of Cross-Trained Workers[END_REF] build a learning model where the quality of service is related to the employee experience. According to their model, the benefits of flexibility are not guaranteed. It is true that a flexible agent can treat more customers, but the quality of service would not be as good as it would be with a dedicated agent. They also compare between different system sizes and show that specialization is preferred in large systems and complete pooling is preferred in small systems. For medium size systems, a mix of flexibility and specialization would be appropriate. In a call center context, [START_REF] Wallace | A Staffing Algorithm for Call Centers with Skill-Based Routing[END_REF] conclude using an extensive simulation study that, when you add skills to an agent, most of the benefits is taken going from one skill per agent to two skills per agent. These results tend to support the idea of limiting the number of skills per agent.

As for the problems of staffing and routing, we refer the reader to the survey by Gans et al. (2003), where the authors present the square-root staffing rule. [START_REF] Borst | Dimensioning Large Call Centers[END_REF] revisited the square-root rule by including principles of routing based on agent costs. To optimize the staffing level in an SBR call center, [START_REF] Henderson | Rostering by Iterating Integer Programming and Simulation[END_REF] combine simulation and integer programming with cutting plane methods. [START_REF] Atlason | Optimizing Call Center Staffing Using Simulation and Analytic Center Cutting-Plane Methods[END_REF] provide interesting properties of a cutting plane method for staffing and prove that it outperforms traditional staffing heuristics which are based on analytical queueing methods. Some other works have investigated the impact of the type of the agent contract on the required staffing level in order to reach a given service level: [START_REF] Ren | Call Center Outsourcing: Coordinating Staffing Level and Service Quality[END_REF] consider piece-meal and pay-per-call-resolved contracts and propose other contracts that coordinate both staffing and effort. [START_REF] Borst | Robust Algorithms for Sharing Agents with Multiple S kills[END_REF] present a routing heuristic that assigns customers to the available agent with the most specialized set of skills which is a generalization of the "specialist-first" principle. [START_REF] Chevalier | Routing and Staffing in Large Call Centers with Specialized and Fully Flexible Servers[END_REF] show that in terms of performance a 20/80 model (20% of generalists and 80% of specialists) performs almost as good as a full flexible model, moreover, it has lower operating costs. In a manufacturing setting, [START_REF] Sheikhzadeh | Machine Sharing in Manufacturing Systems: Total Flexibility versus Chaining[END_REF] evaluate the difference in performance between strict priority, longest queue first and random priority in a comparison between chaining and full flexibility. They remark that the longest queue first policy is the best one. In order to perform a coherent comparison between chaining and our model, we first choose the appropriate routing rules and also staffing levels.

Problem Setting

We consider call center models with n + 1 call types (types 0, 1, ..., n). Customer types 1, 2, ..., n, referred to as also regular types are those requiring specific agent skills 1, 2, ..., n, respectively, while customers 0 can be handled by any agent without a particular "sophisticated" training as required for the regular types. In other words, skill 0 is an easy skill. The mean arrival, service and abandonment rates of customers type i are λ i , µ i and γ i , respectively (i = 0, 1, ..., n). The agents are organized in homogeneous teams, i.e., all agents from a given team have the same set of skills. In this chapter we only consider agent teams with at most two skills per agent. We define an economic framework as follows. We assume that skill 0 costs 1, and that skill i costs 1+t i (for

i = 1, • • • , n).
For two skills i and j, the cost is 1+t i,j (for i, j ∈ {0, • • • , n}). Since skill 0 is the easy skill, we assume that t i,0 ≤ t i,j (for i, j ∈ {0, • • • , n}).

We are interested in the performance in terms of the expected waiting time in the queue of each customer type i taken in service, denoted by W i , for i = 0, 1, ..., n. We denote the objective service level for a type i by W * i , for i = 0, 1, ..., n. In what follows, we describe the two models that we compare in this chapter: chaining and single pooling. We do not consider the architectures developed by [START_REF] Borst | Robust Algorithms for Sharing Agents with Multiple S kills[END_REF], and [START_REF] Bassamboo | Optimal Flexibility Configurations in Newsvendor Networks: Going Beyond Chaining and Pairing[END_REF], because in our context they are too costly or even very hard to implement by a call center manager due to the specificity of the agent skills.

Chaining and single pooling models are shown in Figures 2.1(a) and 2.1(b), respectively. In the chaining model, the skills of the teams are such that they form a chain. The value of the well known chaining model comes from its capability to smooth the workload over the agent teams.

In the single pooling model, customers type 0 benefit from a complete pooling, whereas the other types have only access to one dedicated team. The value of this architecture is that it allows to appropriately handle situations of asymmetry in demand as we will show later. Single pooling can be seen as a dual of the architecture proposed by [START_REF] Chevalier | Routing and Staffing in Large Call Centers with Specialized and Fully Flexible Servers[END_REF], with dedicated single skill agent teams and one team of agents with all skills. The functioning we consider for chaining and single pooling is intuitive and easy to implement in practice. Under chaining, a customer upon her arrival has access to agents from two teams. If at least an agent is available in one of them, then the customer is routed to the team with the higher proportion of idle agents (number of idle agents in a team over the total number of agents in that team). If this proportion is the same for the two teams, then she is equiprobably routed to one of the two teams. Otherwise if all agents from the two teams are busy upon her arrival, the customer waits in her queue (each customer type has its own queue). An agent can handle customers from two queues. Within each queue, the discipline of service is FCFS. When an agent becomes idle, she selects to service one of the customers that are waiting in the two queues, if any. The priority is given to the customer with the longest waiting time.

For single pooling, the routing rules are as shown in Figure 2.1(b). The discipline of service in each one of the n + 1 queues is FCFS. A customer type i (for i = 1, ...n) can be served by only an agent from its associated team. A customer type 0 however can be served by any agent of any one of the n + 1 teams. Upon arrival, a customer type 0 is in priority handled by an idle agent from team 0, if any. If not, she is handled by an idle agent from one of the teams of the regular customer types, if any. If more than one of those have at least one idle agent, then customer 0 is routed to the team with the higher proportion of idle agents. If many teams have the same highest proportion, then customer 0 is equiprobably routed to one of these teams. If all agents of all teams are busy, then customer 0 is placed in her queue. When an agent from one of the teams of the regular customers becomes free, it can serve either a regular customer or a customer 0. However a regular customer has a non-preemptive priority over a customer 0. This means that the idle agent deals with a call from her regular queue first (the first in line). If the queue of the associated regular type is empty, this agent provides service to a customer 0 (the first in line). We assume in our models that the queues are infinite.

In this project, we compare between the two models chaining and single pooling through simulations. In order to have a coherent comparison we optimize their total staffing cost under the constraints W i ≤ W * i , for i = 0, 1, ..., n. We use greedy heuristics for the simulation based optimization step. We refer the reader to the details in Section A.1 of the appendix. For the staffing optimization of SP, we use an increasing greedy algorithm. Starting from an under-staffed situation (a full dedicated model with customers 0), we increase step by step the arrival rate of customers 0. In each iteration, we increment the number of agents in the various teams such that we strictly reach the service level constraints. For chaining, we develop a decreasing greedy algorithm. The algorithm starts with an over-staffed situation using a full dedicated model, which is the worst for chaining since it ignores the links between the teams. We then use an efficient method suggested by [START_REF] Wallace | A Staffing Algorithm for Call Centers with Skill-Based Routing[END_REF] in order to correct the staffing levels to the chaining setting.

Particular Single Pooling Cases

The analytical analysis of the general case of single pooling is too complex. We consider in this section two particular Markovian cases of single pooling, for which, we develop exact and approximate results. The objective of this analysis is to obtain some sense on the effect of the parameters asymmetry on performance. A more comprehensive analysis of the effect of asymmetry is then conducted in Section 2.5 using simulation.

Three Customer Types

Consider a single pooling model with three customer types 0, 1 and 2. The arrival process of types 0, 1 and 2 is Poisson with rates λ 0 , λ 1 and λ 2 , respectively. There two agent teams 1 and 2 with sizes s 1 and s 2 , respectively. The service rate, denoted by µ, is identical for all customer types.

Using a Markov chain approach, we compute in what follows the expected waiting times for all customer types.

Let us define the stochastic process {(x(t), y(t), z(t))t ≥ 0}, where x(t) and y(t) denote the number of busy agents in team i plus the number of waiting customers in queue i (i = 1, 2), and z(t) denotes the number of waiting customers in queue 0, for an instant t ≥ 0. Since inter-arrival and service times are Markovian, {(x(t), y(t), z(t))t ≥ 0} is a Markov chain. Let us denote the system steady-state probabilities by π x,y,z , for x, y, z ∈ N. Note that we can have z ≥ 1 only when x ≥ s 1 and y ≥ s 2 . From the Markov chain, one may write the following set of equations.

We have (λ 1 + λ 2 + λ 0 )π 0,0,0 = µ(π 1,0,0 + π 0,1,0 ). For x = 0 and y > 0, (λ 1 + λ 2 + λ 0 )π 0,y,0 = min(y + 1, s 2 )µπ 0,y+1,0 + µπ 1,y,0 + λ 2 π 0,y-1,0 . For x > 0 and y = 0, (λ 1 + λ 2 + λ 0 )π x,0,0 = min(x + 1, s 1 )µπ x+1,0,0 + µπ x,1,0 + λ 1 π x-1,y,0 . For 0 < x < s 1 or 0 < y < s 2 , there are 5 cases: If y -1 > s 2 s 1 x,

(λ 1 + λ 2 + λ 0 + min(x, s 1 )µ + min(y, s 2 )µ)π x,y,0 = min(x + 1, s 1 )µπ x+1,y,0 + min(y + 1, s 2 )µπ x,y+1,0 + (λ 1 + λ 0 )π x-1,y,0 + λ 2 π x,y-1,0 . If y -1 = s 2 s 1 x, (λ 1 + λ 2 + λ 0 + min(x, s 1 )µ + min(y, s 2 )µ)π x,y,0 = min(x+1, s 1 )µπ x+1,y,0 +min(y+1, s 2 )µπ x,y+1,0 +(λ 1 +λ 0 )π x-1,y,0 +(λ 2 +λ 0 /2)π x,y-1,0 . If y-1 < s 2 s 1 x and y > s 2 s 1 (x -1), (λ 1 + λ 2 + λ 0 + min(x, s 1 )µ + min(y, s 2 )µ)π x,y,0 = min(x + 1, s 1 )µπ x+1,y,0 + min(y + 1, s 2 )µπ x,y+1,0 + (λ 1 + λ 0 )π x-1,y,0 + (λ 2 + λ 0 )π x,y-1,0 . If y = s 2 s 1 (x -1), (λ 1 + λ 2 + λ 0 + min(x, s 1 )µ + min(y, s 2 )µ)π x,y,0 = min(x + 1, s 1 )µπ x+1,y,0 + min(y + 1, s 2 )µπ x,y+1,0 + (λ 1 + λ 0 /2)π x-1,y,0 + (λ 2 + λ 0 )π x,y-1,0 . If y < s 2 s 1 (x -1), (λ 1 + λ 2 + λ 0 + min(x, s 1 )µ + min(y, s 2 )µ)π x,y,0 = min(x + 1, s 1 )µπ x+1,y,0 + min(y + 1, s 2 )µπ x,y+1,0 + λ 1 π x-1,y,0 + (λ 2 + λ 0 )π x,y-1,0 . For x = s 1 , y = s 2 and z = 0, (λ 1 + λ 2 + λ 0 + (s 1 + s 2 )µ)π s 1 ,s 2 ,0 = s 1 µπ s 1 +1,s 2 ,0 + s 2 µπ s 1 ,s 2 +1,0 + (λ 1 + λ 0 )π s 1 -1,s 2 ,0 + (λ 2 + λ 0 )π s 1 ,s 2 -1,0 + (s 1 + s 2 )µπ s 1 ,s 2 ,1 . For x = s 1 and y > s 2 , (λ 1 + λ 2 + λ 0 + (s 1 + s 2 )µ)π s 1 ,y,0 = s 1 µπ s 1 +1,y,0 + s 2 µπ s 1 ,y+1,0 + (λ 1 + λ 0 )π s 1 -1,y,0 + λ 2 π s 1 ,y-1,0 + s 1 µπ s 1 ,y,1 . For x > s 1 and y = s 2 , (λ 1 + λ 2 + λ 0 + (s 1 + s 2 )µ)π x,s 2 ,0 = s 1 µπ x+1,s 2 ,0 + s 2 µπ x,s 2 +1,0 + λ 1 π x-1,s 2 ,0 + (λ 2 + λ 0 )π x,s 2 -1,0 + s 2 µπ x,s 2 ,1 . For x > s 1 , y > s 2 and z = 0, (λ 1 + λ 2 + λ 0 + (s 1 + s 2 )µ)π x,y,0 = s 1 µπ x+1,y,0 + s 2 µπ x,y+1,0 + λ 1 π x-1,y,0 + λ 2 π x,y-1,0 . For z > 0, (λ 1 + λ 2 + λ 0 + (s 1 + s 2 )µ)π s 1 ,s 2 ,z = s 1 µπ s 1 +1,s 2 ,z + s 2 µπ s 1 ,s 2 +1,z + λ 0 π s 1 ,s 2 ,z-1 + (s 1 + s 2 )µπ s 1 ,s 2 ,z+1 . For x = s 1 , y > s 2 and z > 0, (λ 1 + λ 2 +λ 0 +(s 1 +s 2 )µ)π s 1 ,y,z = s 1 µπ s 1 +1,y,z +s 2 µπ s 1 ,y+1,z +λ 0 π s 1 ,y,z-1 +s 1 µπ s 1 ,y,z+1 . For x > s 1 , y = s 2 and z > 0, (λ 1 + λ 2 + λ 0 + (s 1 + s 2 )µ)π x,s 2 ,z = s 1 µπ x+1,s 2 ,z + s 2 µπ x,s 2 +1,z + λ 0 π x,s 2 ,z-1 + s 2 µπ x,s 2 ,z+1 . For x > s 1 , y > s 2 and z > 0, (λ 1 +λ 2 +λ 0 +(s 1 +s 2 )µ)π x,y,z = s 1 µπ x+1,y,z +s 2 µπ x,y+1,z +λ 0 π x,y,z-1 .
One may intuitively see from the Markov chain how the asymmetry in arrivals increases the performance of single pooling. The counterproductive states are those with waiting customers and idle agents at the same time, i.e., x > s 1 and 0 ≤ y < s 2 , or y > s 2 and 0 ≤ x < s 1 . When being in one of these two "bad" cases, the probabilities to take the direction of leaving them are λ 0 +λ 2 +s 1 µ λ 1 +λ 2 +λ 0 +(s 1 +y)µ and λ 0 +λ 1 +s 2 µ λ 1 +λ 2 +λ 0 +(s 1 +y)µ , respectively. This shows for example that increasing the proportion of customers 0 (one form of asymmetry) increases the performance in single pooling.

A further illustration is given next. The expected waiting times as a function of the steady-state probabilities are given by

W 1 = 1 λ 1   +∞ ∑ x=s 1 +∞ ∑ y=0 (x -s 1 )π x,y,0 + +∞ ∑ x=s 1 +∞ ∑ y=s 2 +∞ ∑ z=1 (x -s 1 )π x,y,z   , W 2 = 1 λ 2 ( +∞ ∑ y=s 2 +∞ ∑ x=0 (y -s 2 )π x,y,0 + +∞ ∑ y=s 2 +∞ ∑ x=s 1 +∞ ∑ z=1 (y -s 2 )π x,y,z ) , W 0 = 1 λ 0 +∞ ∑ z=0 +∞ ∑ x=s 1 +∞ ∑ y=s 2 zπ x,y,z .
The performance measures above are computed numerically. We solve the steady-state equations relating the state probabilities using a state space truncation, with a sufficiently high precision (six digits beyond the decimal point

= µ 2 = µ 0 = 0.2, ∑ 2 i=0 λ i = 4, λ 1 = λ 2 , s 1 = s 2 = 12)

A Fixed Point Approximation

We consider here a Markovian single pooling case with an arbitrary number of skills. There are n + 1 customer types (type 0, and types 1, 2, ..., n), n teams (no team 0), n ≥ 1. The arrival rates are λ 0 and λ i = λ for i = 1, ..., n, and the service rates are µ i = µ for i = 0, 1, ..., n. Since the configuration is symmetric, we consider the same staffing level s in each team. In what follows, we develop an approximation to compute the expected waiting time of regular customers type i, for i = 1, ..., n. The approximation is based on a Markov chain approach and a fixed point algorithm.

One can see that our model can be divided into n identical sub-systems. It suffices then to focus on the performance analysis of one of these sub-systems. A sub-system is a simple queueing system with s servers and an infinite queue. Two types of customers arrive to this sub-system: customers type i with a Poisson process with rate λ and customers type 0 with a general arrival process with mean arrival rate λ 0 n . (The arrival process of customers 0 to the whole system is Poisson. However, it becomes a general process at each sub-system because of the routing rules.)

Recall that customers 0 wait in their own queue before being routed to one of the sub-systems for an immediate processing. Because of the routing rule, customers 0 can be routed to a sub-system only if the number of customers in the sub-system is less or equal to s -1. Also since we route customers 0 to the one of the less busiest sub-systems (with an equiprobable choice), the arrival rate of customers 0 is decreasing in the number of busy servers in a sub-system and it becomes 0 when all the s servers become busy.

Let us now define, for a sub-system, the stochastic process {E(t), t ≥ 0}, where E(t) denotes the number of customers in the system (queue + service). Note that the customers in the queue are only the regular customers, and those in service can be both regular or type 0 customers. We approximate customers 0 inter-arrival times by an exponential distribution with state-dependent rates. Since inter-arrival and service times are Markovian, {E(t), t ≥ 0} is a Markov chain as shown in Figure 2.3. The arrival rate δ k denotes the state-dependent arrival rate of customers 0 when the number of customers in the sub-system is k, for k = 0, ..., s -1 (no customers 0 arrive at the sub-system for k ≥ s).

Assume that exactly s customers are in the sub-system and that a service completion occurs first before the next arrival epoch of a regular customer at this sub-system (Figure 2.3). Therefore, two possibilities may happen. The first possibility corresponds to the case of an empty queue 0.

0 1 2 1 s - s 1 s + 2 s + µ 2µ 3µ ( 1) s µ - (1 ) sµ β - sµ sµ sµ 0 λ δ + 1 λ δ + 2 λ δ + 2 s λ δ - + 1 s λ δ - + λ λ λ
Figure 2.3: Markov chain associated to a sub-system of single pooling

We then move to state s -1. The second one corresponds to the case of a non-empty queue 0. We then stay in state s, because the server who just became idle immediately takes the customer 0 in the head of queue 0 into service. Let us denote by β the probability that queue 0 is not empty.

Then the rate to move from state s to sate s -1 in the Markov chain is sµ(1 -β).

Let us now assume that the stability condition of a sub-system holds, i.e., λ + λ 0 n < sµ, and denote the stationary probabilities of the system states by π k , for k ≥ 0. We may then write

π k = ∏ k-1 i=0 (λ + δ i ) k!µ k π 0 , (2.1)
for 1 ≤ k ≤ s -1, and

π s+k = ( λ sµ ) k ∏ s-1 i=0 (λ + δ i ) s!µ s (1 -β) π 0 , (2.2)
for k ≥ 0. Since all probabilities sum up to one, we obtain

π 0 = 1 1 + ∑ s-1 k=1 ∏ k-1 i=0 (λ+δ i ) k!µ k + ∏ s-1 i=0 (λ+δ i ) s!µ s (1-β) 1 1-λ sµ .
(2.

3)

The difficulty to compute the stationary probabilities is that we do not have the values of δ k (k = 0, ..., s -1) and β. We use a fixed point algorithm to jointly compute them with the stationary probabilities. Let us now write δ 0 , the arrival rate of customers 0 at a given sub-system when this sub-system is empty, as a function of the stationary probabilities of this sub-system. We use here a second approximation. We indeed assume that the states of the sub-systems are independent, which is not true. Assume that our sub-system is the only one that is empty, i.e., each one of the other n -1 sub-systems have at least one customer in the system (queue + service). Using the approximation this occurs with probability (1 -π 0 ) n-1 , then δ 0 is simply λ 0 in that case. Assume now that our sub-system and only another one are empty. Then δ 0 is λ 0 2 (equiprobable routing of customers 0 to one of the less busiest sub-systems). This occurs with probability π 0 (1 -π 0 ) n-2

and there are ( n-1

1

) combinations (where

( n k ) = n! k!(n-k)! for 0 ≤ k ≤ n).
Continuing with the same reasoning and averaging over all possibilities, we obtain

δ 0 = λ 0 n-1 ∑ j=0 1 j + 1 ( n -1 j ) π j 0 (1 -π 0 ) n-1-j . (2.4) Since 1 j+1 ( n-1 j ) = 1 n ( n j+1
) , Equation (2.4) becomes

δ 0 = λ 0 n n-1 ∑ j=0 ( n j + 1 ) π j 0 (1 -π 0 ) n-1-j = λ 0 nπ 0 n ∑ j=1 ( n j ) π j 0 (1 -π 0 ) n-j ,
which leads to

δ 0 = λ 0 1 -(1 -π 0 ) n nπ 0 .
(2.5)

In the same way, we obtain

δ k = λ 0 ( 1 - ∑ k-1 j=1 π j ) n - ( 1 - ∑ k j=1 π j ) n nπ k , (2.6) for 1 ≤ k ≤ s -1.
Let us now give the expression of β as a function of the stationary probabilities π k , k ≥ 0. Since the mean arrival rate of customers 0 at our sub-system is λ 0 n , we have

s-1 ∑ k=0 π k δ k + βsµ = λ 0 n , (2.7) which implies β = λ 0 n - ∑ s-1 k=0 π k δ k sµ .
(2.8)

In summary, from the one hand, Equations (2.1)-(2.3) give the stationary probabilities π k (k ≥ 0) as a function of δ k (0 ≤ k ≤ s -1) and β. From the other hand, Equations (2.5), (2.6) and (2.8) give δ k (0 ≤ k ≤ s -1) and β as a function of π k (k ≥ 0). As a consequence, we have a fixed point. We propose the following fixed point algorithm to compute it. In the first iteration, we choose δ 0 = λ 0 n , δ k = 0 for 1 ≤ k ≤ s -1, and β = 0. Then we compute π k (k ≥ 0) using Equations (2.1)-(2.3). From these π k , we next compute the new values of δ k (0 ≤ k ≤ s -1) and β using Equations (2.5), (2.6) and (2.8). In the second iteration, we use the latter values of δ k and β to compute π k . From these new π k , we compute the new values of δ k and β. We do the same in the third iteration, and so on. We stop the algorithm when the values of π k (k ≥ 0), δ k (0 ≤ k ≤ s -1) and β converge to their limits with a given predefined precision (we have chosen a precision of 10 -6 in the numerical experiments below). Proposition 1 proves the convergence of the fixed point algorithm.

Proposition 1 The fixed point algorithm always converges.

Proof. We use the Brouwer's theorem to prove the convergence. The Brouwer's theorem states that any continuous function from a convex compact subset K of an Euclidean space to itself has at least one fixed point. In what follows, we prove that the conditions of the Brouwer's theorem hold in our context.

After k iterations, the fixed point algorithm gives the vector (π 0 , π 1 , π 2 , • • • , π c ) k belonging to a convex compact, [0; 1] s+1 , that is included in an Euclidean space, R s+1 . From Equations (2.1)-

(2.8), it is obvious to see that the function that allows to calculate (π

0 , π 1 , π 2 , • • • , π c ) k+1 (iteration k + 1) as a function of (π 0 , π 1 , π 2 , • • • , π c ) k is continuous (combination of continuous functions),
for π k ̸ = 0 (k = 0, ..., s -1). In what follows, we prove that this function is continuous in π k = 0 (k = 0, ..., s -1) by prolongation. From Equations (2.5) and (2.6), we have

δ k = λ 0 ( 1 - ∑ k-1 j=1 π j ) n - ( 1 - ∑ k j=1 π j ) n nπ k = λ 0 ( 1 - ∑ k-1 j=1 π j ) n nπ k   1 - ( 1 - ∑ k-1 j=1 π j ) n ( 1 - ∑ k-1 j=1 π j ) n   ,
for k = 0, ..., s -1, where by convention an empty sum is equal to 0. Calculating further, we obtain

δ k = λ 0 ( 1 - ∑ k-1 j=1 π j ) n nπ k ( 1 - ( 1 - π k 1 - ∑ k-1 j=1 π j ) n ) ,
for k = 0, ..., s -1. The Taylor expansion of δ k as a function of π k in the neighborhood of 0 is

δ k = λ 0 ( 1 - ∑ k-1 j=1 π j ) n nπ k (1 -(1 -n o(π k ))) = λ 0   1 - k-1 ∑ j=1 π j   n + o(1),
where o(1) is a function that converges to a finite limit as π k goes to 0, for k = 0, ..., s -1. Since

λ 0 ( 1 - ∑ k-1 j=1 π j
) n is finite, δ k is continuous by prolongation in π k = 0, for k = 0, ..., s -1.

It remains now to focus on the issue for β = 1 in Equation (2.3). This case of β = 1 can not happen. The proof is as follows. Assume that β = 1. Equation (2.7) thus leads to sµ = λ 0 n -

∑ s-1 k=0 π k δ k .
Since δ k and π k (0 ≤ k ≤ s -1) are positive, sµ ≤ λ 0 n . As a consequence the sub-system is unstable, which is absurd. This completes the proof of the convergence of the fixed point algorithm. ✷

Having in hand the stationary probabilities, we next compute for the regular customers the expected waiting time in the queue and the probability of delay. Recall that all sub-systems are identical because of the symmetry in the parameters. Using Little's law, the expected waiting time of a regular customer type i (i = 1, ..., n) is given by

W i = 1 λ ∞ ∑ k=1
kπ s+k , (2.9) for i = 1, ..., n, and its probability of delay denoted by P D,i is

P D,i = ∞ ∑ k=1 π s+k , (2.10)
for i = 1, ..., n. The approximation for both W i and P D,i works very well for the regular customer types, however it does not for customers 0 because of their complex routing. The comparison between the approximate results using the fixed point algorithm and the exact ones using simulation are given in Table 2.1. Note that in the extreme situations of n = 1 or λ 0 = 0, our method gives the exact results.

Table 2.1 reveals that our approximation yields very accurate estimates, while slightly overestimating the performance. It gives lower values for W i and P D,i than those from simulation.

An explanation would be as follows. In our approximation, we assume that all sub-systems are independent one of another. In reality, the routing rule leads to a fair sharing of customers between the sub-systems. Therefore, when a given sub-system is almost busy, the other ones are likely to be almost busy. Thus, the arrival rates δ k (for high values of k close to s -1) should be in reality higher than those in the approximation, which implies that the latter would give better performance (lower waiting and lower probability of delay) than simulation does.

Using the above approximate analysis, we illustrate in Figure 2.4 how the asymmetry in arrivals (by increasing p) improves performance (expected waiting time W i for i = 1, ..., n).

µ i = µ 0 = 0.2, ∑ 2 i=0 λ i = 8, λ i = λ j , s i = 12 for i, j = 1, ..., 4)

Effect of the Parameters Asymmetry

In this section, we present the results of the comparison between chaining and single pooling. We use simulation experiments to optimize the call center staffing. In using simulation for call center operations management, we are following longstanding practice, see for example [START_REF] Wallace | A Staffing Algorithm for Call Centers with Skill-Based Routing[END_REF].

We simplify the cost model such that the SP cost is upper bounded and that of chaining is lower bounded. All the numerical comparisons are based on the lower and upper bounds values. This makes the results pessimistic for SP and optimistic for chaining, i.e., the performance of SP is in reality better than what we present. The cost of single pooling is

∑ n i=0 (1 + t i,0 )s i . This is upper bounded by ( ∑ n i=0 s i ) max i (1+t i,0 ). The cost of chaining is (1+t 0,1 )s 0 +(1+t 1,2 )s 1 +• • • (1+t n,0 )s n
and is lower bounded by (1 + t 0,1 )s 0 + (1 + min i,j (1

+ t i,j ))( ∑ n-1 i=1 s i ) + (1 + t n,0 )s n .
Let us now simplify the problem as follows. An agent with skills 0 and i (i = 1, ..., n) costs 1. An agent with skills i and j (i, j = 1, ..., n and i ̸ = j) costs 1 + t, t ≥ 0. In this simplification, we have max i (1 + t i,0 ) = 1 and min i,j (1 + t i,j ) = 1 + t (i, j = 1, ..., n and i ̸ = j). The parameter t is then the incremental cost of an agent with two regular skills compared to that with a regular skill and skill 0.

Design of Experiments. As we are interested in the effect of asymmetry of the parameters on performance, we propose various forms of asymmetry. For customers 0, we define the parameters p and p ′ to measure the relative importance in arrivals and service durations, respectively. They are given by p = λ 0

∑ n i=0 λ i and p ′ = 1 µ 0 ∑ n i=0 1 µ i
. We measure the asymmetry between the arrival rates of regular customers by

V = λ 1 λ 2 = λ 2 λ 3 = • • • = λ n-1
λn , and that between service durations by

U = 1/µ 1 1/µ 2 = 1/µ 2 1/µ 3 = • • • = 1/µ n-1
1/µn . We also consider for customers 0 the asymmetry in the variability of service times, measured by the coefficient of variation of its distribution and denoted by cv s .

We consider other forms of asymmetry in terms of the required service level and also the time to abandon for customers 0 relatively to those for the regular customers. These effects are studied in the settings of small and large call centers, and also in the settings of small and large number of skills. Although the considered forms of asymmetries do not cover all possibilities, they allow to obtain the main useful conclusions.

The approach to conduct the simulation experiments is as follows. Due to the high number of parameters, we first run experiments by separately treating one parameter at a time. In a systematic way, we vary one parameter while holding all the others constant. Second to see the possible interaction effects, we simultaneously vary the values of more than one of them at a time.

For the values of the parameters, we choose wide ranges that allow to cover most of call center situations in practice. For the rest of the chapter, inter-arrival are assumed to be Markovian.

Service times are also assumed to be Markovian, except in Section 2.5.2. The abandonment rates are assumed to be equal to zero, except for Section 2.5.4.

Asymmetry in Arrival Rates

We want to understand the effect of the asymmetry in the demand. We separate the study into two steps. First, we construct the asymmetry only on the arrival rate of customers 0. Second, we construct it by differentiating between all the arrival rates of all customer types.

Asymmetry on Customers 0

To isolate the impact of p = λ 0 ∑ n i=0 λ i , we assume that all customer types have the same expected service time, and all the arrival rates of the regular customers are the same, λ i = λ for i = 1, ..., n Since any agent in SP has skills 0 and i (i.e., costs 1), the staffing cost of SP does not depend on t. In Table 2.2, the column Crossing value gives the value of t for which the two models chaining and SP are equivalent. Below this threshold chaining is better than SP and viceversa (see Figure 2.5(a)). Consider small values of t.

(V =
= µ 0 = 0.2, W * 0 = W * i = 0.2, ∑ 4 i=0 λ i = 8, i = 1, ..., 4, U = V = 1, p ′ = 20%, n = 4)
= µ 0 = 0.2, W * 0 = W * i = 0.2, ∑ 4 i=0 λ i = 8, i = 1, ..., 4, U = V = 1, p ′ = 20%, n = 4)
from pooling, which highly improves the performance of SP. With t = 0, SP and chaining become equivalent for values of p ≥ 75%.

For higher values of t, SP goes ahead of chaining. The reason is related to the increase of the costs of the agents with two skills i and j (i, j = 1, ..., 4). It suffices to have t = 15.38% to outperform the best performance of chaining (the symmetric case). For any t beyond 30%, SP is systematically better than chaining whatever is p.

We also measure the relative benefits between SP and chaining. Figure 2.5(b) provides, for various values of the relative benefits, the associated curve of t as a function of p. We observe that the sensitivity of the relative benefit as a function of t decreases in p. The reason is that the number of customers 0 increases in p, which decreases the number of agents with two regular skills in chaining (i.e., decreases the cost sensitivity in t).

The main conclusion here is that SP can be better than chaining when the demand for skill 0 is important and/or when skill 0 is less costly than the other ones. The main idea is that as type 0 dominates, they profit in SP from a total pooling from all teams, while chaining, they do profit from a partial pooling from only two teams.

Asymmetry on the other Arrival Rates

The parameter p, that is defined on customers 0, is one way of measuring asymmetry in arrivals.

Here, we focus on the asymmetry between regular customer types, measured by

V = λ 1 λ 2 = λ 2 λ 3 = λ 3 λ 4 .
The simulation results for the cases V = 2 and 5 are shown in Table 2.3. The experiments for the case V = 1 reduces to those given in Table 2.2. 2.6 reveal that the performance of SP increases in V . An intuitive explanation is as follows. Remark that the team size s i = s(λ i ) is increasing and concave in λ i , for i = 1, ..., n. Applying then the Jensen inequality leads to

Table 2.3: Impact of V (µ i = µ 0 = 0.2, W * 0 = W * i = 0.2, ∑ 4 i=0 λ i = 8, i = 1, ..., 4, U = 1, p ′ = 20%, n = 4)
= µ 0 = 0.2, W * 0 = W * i = 0.2, ∑ 4 i=0 λ i = 8, i = 1, ..., 4, U = 1, p ′ = 20%, n = 4)
n ∑ i=1 s(λ i ) ≤ n • s ( ∑ n i=1 λ i n
)

. In this inequality, the left hand side corresponds to the overall staffing level for an arbitrary situation, i.e., with arbitrary values of λ i s. As for the right hand side, it gives the overall staffing level for a symmetric situation, i.e., all the λ i s are identical. We also observe from Table 2.3 that the performance of chaining is however relatively insensitive to V . Note that we change each time the configuration of chaining such that the large teams are close to each others in order to create more pooling effect. This is better than having small teams each of which connected to a large team.

Asymmetry in Service Rates

In this section, we focus on the comparison between chaining and SP with regard to the asymmetry in the customer service times. We first define the asymmetry only on customers 0, and second on all customer types.

Asymmetry on Customers 0

We measure the asymmetry on customers 0 by p ′ =

1 µ 0 ∑ n i=0 1 µ i
. The asymmetry here is defined by the difference between the value of the mean service time of customers 0 and that of the regular types. The results are shown in Table 2.4 and Figure 2.7(a).

Table 2.4: 2.4, we observe that the performance of both models chaining and SP improves in p ′ (from 0 until the symmetric case for p ′ = 25%). The reason is that for chaining we are approaching the symmetric case where it behaves well, and for SP we are profiting better from the pooling effect when all service times are statistically identical. However the performance of the two models deteriorates in p ′ (for p ′ above 25%), and no model performs well for a high asymmetry in service times. The explanation is related to a phenomenon referred to as the blocking effect. The blocking effect is the situation where the agents are excessively blocked by customers 0 (who are in need of large service times) which deteriorates the waiting time of the regular customers. This phenomenon is more apparent for single pooling since in the latter customers 0 have access to all teams, whereas in chaining they do only have access to two teams. We refer the reader to [START_REF] Tekin | Pooling Strategies for Call Center Agent Cross-Training[END_REF] for more details on how pooling could be counterproductive when service times are very different.

Impact of p ′ (λ i = λ 0 = 2, ∑ 4 i=0 1 µ i = 25, W 0 = W * i = 0.2, i = 1, ..., 4, p = 20%, U = V = 1, n = 4)
µ i = 25, W0 = W * i = 0.2, i = 1, ..., 4, p = 20%, U = V = 1, n = 4)
µ i = 25, W0 = W * i = 0.2, i = 1, ..., 4, p = 20%, U = V = 1, n = 4) Ϭй ϱй ϭϬй ϭϱй ϮϬй Ϯϱй Ϭ ϭ Ϯ ϯ ϰ ϱ ƚ Đ|Ɛ ŚĂŝŶŝŶŐ ^ŝŶŐůĞ WŽŽůŝŶŐ (c) Effect of variability in service times(λ0 = 2, λi = 1.5, µ0 = µi = 0.2, W0 = W * i = 0.2, i = 1, ..., 4, p = 25%, p ′ = 20%, U = V = 1, n = 4) ͲϯϬй ͲϮϬй ͲϭϬй Ϭй ϭϬй ϮϬй ϯϬй ϰϬй Ϭ ϭ Ϯ ϯ ϰ ϱ ͲϭϬй Ͳϱй Ϭй нϱй нϭϬй ƚ Đ|Ɛ ŚĂŝŶŝŶŐ ^ŝŶŐůĞ WŽŽůŝŶŐ ZĞůĂƚŝ|Ğ ďĞŶĞĨŝƚƐ ŽĨ ^W (d) Relative benefits (λ0 = 2, λi = 1.5, µ0 = µi = 0.2, W0 = W * i = 0.2, i = 1, ..., 4, p = 25%, p ′ = 20%, U = V = 1, n = 4)
Recall that this situation with a slow service rate for customers 0 is out of our context. In our context, customers 0 are in need of an easy skill, and are therefore likely to be served within a short duration. We also measure the relative benefits between SP and chaining as a function of p ′ (see Figure 2.7(a)). Similarly to the effect of p, we again observe that the sensitivity of the relative benefit as a function of t decreases in p ′ . The reason is that the service time duration of customers 0 increases in p ′ , which decreases the number of agents with two regular skills in chaining (i.e., decreases the cost sensitivity in t).

In what follows, we go further by defining the asymmetry on the variability of customers 0 service times. We choose to measure this variability by the coefficient of variation (ratio of standard deviation over expected value), denoted by cv s . We consider a log-normal distribution for the service times of customers 0 (inter-arrival times of all types, and service times of all regular types are Markovian). The choice of the log-normal distribution is based on the call center statistical analysis in [START_REF] Brown | Statistical Analysis of a Telephone Call Center: A Queueing-Science Perspective[END_REF]. The results are shown in Table 2.5 and Figure 2.7(c). We draw the same conclusions as those for service rates. Due to the blocking effect, both models do not behave well as the variability is increasing. Figure 2.7(d) reveals that the relative benefit as a function of t is not sensitive to the variation of cv s . To the contrary to the case for p and p ′ , the arrival and service rates of regular types do not vary here.

= µ 0 = 0.2, W 0 = W * i = 0.2, i = 1, ..., 4, p = 25%, p ′ = 20%, U = V = 1, ∑ 4 i=0 λ i = 8, n = 4)

Asymmetry on the Other Service Rates

We examine the impact of asymmetry by defining it on all service times. The service times can be now different from one regular customer to another. Recall that the ratio U is defined by

U = 1/µ 1 1/µ 2 = 1/µ 2 1/µ 3 = 1/µ 3 1/µ 4 .
We also consider cases with a high proportion of customers 0, p = 50%.

This can be seen as a worst case for SP, since the blocking effect is more apparent in such a case.

The simulation results are shown in

= 0.2, λ 0 = 4, λ i = 1, W 0 = W * i = 0.2, i = 1, ..., 4, ∑ 4 i=0 1 µ i = 25, p ′ = 20%, p = 50%, V = 1, n = 4)
µ i = 25, p ′ = 20%, p = 50%, V = 1, n = 4)
From the numerical results we observe that SP is preferred to chaining for a wide range of parameters. The performance of SP is quite insensitive to the asymmetry defined by U . The reason is that whatever is U , the agent teams in SP are divided to two types. One first type with two teams where customers 0 are served faster than regular customers (positive effect), and a second type with two teams where customer 0 are served slower than regular customers (negative effect of blocking). The performance of chaining is however decreasing in asymmetry. In chaining, each team receives two customer types with different service times, which creates a negative blocking effect in all teams and deteriorates as a consequence the performance. In general for both single pooling and chaining with U ̸ = 1, regular customers require different mean service times. We then have regular customers that are served faster than others. The slowly served ones block the teams in which they are routed to. This is more apparent in chaining because regular customers are routed to two teams (and to only one in SP). We also measure the relative benefits between SP and chaining. Figure 2.8(b) reveals that this benefit as a function of t is not sensitive to the variation of U . The reason is that although the service rates of regular types do vary, the total staffing level for the regular types do not.

Asymmetry in the Service Level Constraints

We define the asymmetry on the service level of customers 0, W * 0 . The results are shown in Table 2.7 and Figure 2.9(a).

Table 2.7: We observe as expected that SP behaves better than chaining in the case of a high asymmetry in the service levels. Chaining is requiring higher staffing levels than needed for some customer types. The agent teams are less correlated in SP than in chaining. This gives more flexibility under SP to adjust the size of the teams as required. However, the strong link between the chains in chaining forces the size of the teams to be adjusted with regard to the high requirement of some customer types while it is not needed for other types. As for the relative benefits between SP and chaining, we observe from Figure 2.9(b) that it is not sensitive to the variation of W * 0 . Since the parameters related to the regular types do not vary, the associated staffing levels do not change also.

Impact of W * 0 (λ 0 = 4, λ i = 1, µ i = µ 0 = 0.2 and W * i = 0.2 for i = 1, ..., 4, p = 50%, p ′ = 20%, U = V = 1, n = 4)
λ i = 1, µ i = µ 0 = 0.2 and W * i = 0.2 for i = 1, ..., 4, p = 50%, p ′ = 20%, U = V = 1, n = 4)

Asymmetry in Abandonments

We allow in this section customers to abandon. After entering the queue, a customer will wait a random length of time for service to begin. If service has not begun by this time she will abandon and be lost. Abandonment is an important feature in call centers. We first investigate the impact of abandonment on the performance of single pooling and chaining. We then investigate the effect of the asymmetry in the abandonment rate of customers 0. Recall that that γ i denotes the abandon rate of customers i, for i = 0, • • • n. In the experiments below, times before abandonment are assumed to be exponentially distributed. Note that with customer abandonment, new performance measures do appear for waiting times. Since the sojourn time in queue may end up with a start of service or an abandonment, we distinguish the conditional waiting time given service, that given abandonment, and the unconditional one. We focus here on the conditional waiting time given service.

Impact of Abandonment. We investigate the impact of abandonment on the performance of SP and chaining in various situations of asymmetries. We consider homogeneous abandonments for all customer types, γ i = γ for i = 0, • • • n. The results are shown in Figures 2.10(a)-2.10(d). Further results are also given in Tables A.6-A.9 in Section A.3 of the appendix. An important observation here is that the effect of the parameters asymmetry changes in the presence of abandonment. For example, to the contrary to the results with no abandonment, the performance of SP deteriorates in p, but improves in p ′ . The reason is that the abandonment of customers reduces the arrivals to service, which in turn reduces the asymmetry. This can be seen from Table 2.8, where the the probability to abandon of customers 0 increases in p. .11(d). Further results are also given in Tables A.10-A.13 in Section A.3 of the appendix. We again observe an important impact of the abandonment on the performance of SP and chaining. This impact mainly depend on how the abandonment affects the asymmetry. For example, we observe from Figure 2.11(a) that when regular customers have higher abandonment rates than customers 0, the asymmetry in terms of p is accentuated (which further improves SP performance). In the opposite case however, the asymmetry in p reduces because of the abandonment of customers 0.

(a) Impact of p (µi = µ0 = 0.2, ∑ 2 i=0 λi = 8, λi = λj, W * 0 = W * i = 0.2, for i, j = 1, ..., 4, p ′ = 20%, U = V = 1, n = 4) 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% p' t γ=0 γ=0.1 γ=0.2 Single Pooling Chaining (b) Impact of p ′ (λi = λ0 = 2, ∑ 4 i=0 1 µ i = 25, W0 = W * i = 0.2, i = 1, ..., 4, p = 20%, U = V = 1, n = 4) 0% 5% 10% 15% 20% 25% 30% 35% 40% 1 1.5 2 2.5 3 3.5 4 4.5 5 V t γ=0 γ=0.1 γ=0.2 Single Pooling Chaining (c) Impact of V (λ0 = 2, µ0 = µi = 0.2, ∑ 4 i=0 λi = 8, W0 = W * i = 0.2, i = 1, ..., 4, p = 25%, p ′ = 20%, U = 1, n = 4) 0% 5% 10% 15% 20% 25% 30% 35% 40% 1 1.5 2 2.5 3 3.5 4 4.5 5 U t γ=0 γ=0.1 γ=0.2 Single Pooling Chaining (d) Impact of U (µ0 = 0.2, λ0 = 4, λi = 1, W0 = W * i = 0.2, i = 1, ..., 4, ∑ 4 i=0 1 µ i = 25, p ′ = 20%, p = 50%, V = 1, n = 4)
= µ 0 = 0.2, ∑ 2 i=0 λ i = 8, λ i = λ j , W * 0 = W * i = 0.2, γ i = γ 0 = γ for i, j = 1, ..., 4, p ′ = 20%, U = V = 1, n = 4) γ = 0.1 γ = 0.
(c) Impact of V (λ0 = 2, W * 0 = W * i = 0.2, µi = µ0 = 0.2, i = 1, ..., 4, p = 25%, p ′ = 20%, U = 1, n = 4) Ϭй ϱй ϭϬй ϭϱй ϮϬй Ϯϱй ϯϬй ϯϱй ϰϬй ϰϱй ϱϬй ϭ Ϯ ϯ ϰ ϱ ɶϬхɶŝ ɶϬфɶŝ ƚ ŚĂŝŶŝŶŐ ^ŝŶŐůĞ WŽŽůŝŶŐ h (d) Impact of U (λ0 = 4, W * 0 = W * i = 0.2, µ0 = 0.2, i = 1, ..., 4, p = 50%, p ′ = 20%, V = 1, n = 4)
= µ 0 = 0.2, W * i = W * 0 = 0.2 for i = 1, ..., 4, p ′ = 20%, U = V = 1, n = 4) Small Call Center ( ∑ 4 i=0 λi = 1) Large Call Center ( ∑ 4 i=0 λi =

Impact of the Call Center Size

We focus in this section on the impact of the size of the call center on the comparison between the two models. [START_REF] Akşin | Characterizing the performance of process flexibility structures[END_REF] showed that a small call center benefits more from a flexible architecture than a larger one. From the simulation experiments conducted here, we confirm this conclusion. The results are shown in Table 2.9 and Figure 2.12. In Table 2.10 provides the achieved expected waiting times for the optimal staffing levels.

Table 2.10: Expected waiting times ( Because of the small teams, the lack of the pooling effect in small call centers makes the threshold values of t higher than those in large call centers. However in large call centers, the team sizes are quite large in the sense that we have a less need to the chains. This makes SP better than chaining even under the symmetric case of arrival rates. From Table 2.10 we observe that to the contrary to small call centers, the service level constraints are saturated for large call centers. Because of

µ i = µ 0 = 0.2, W * i = W * 0 = 0.2 for i = 1, ..., 4, p ′ = 20%, U = V = 1, n = 4)
µ i = µ 0 = 0.2, W * i = W * 0 = 0.2 for i = 1, ..., 4, p ′ = 20%, U = V = 1, n = 4)
the discrete nature of staffing levels, the impact of adding or removing an agent on performance is higher in small call centers. For the same reason, the staffing levels of the regular teams do not vary much in small call centers. This makes the relative benefits between SP and chaining not sensitive to the variation of p in small call centers, while the opposite is true for large call centers (see Figures 2.12(b) and 2.12(d)).

Impact of the Number of Skills

In this section, we investigate the effect of the number of skills (denoted by N = n + 1). For two cases with different number of skills, it is not possible to keep at the same time a constant workload on each team and a constant overall workload. We choose to separately treat each situation.

Constant Workload per Team. We consider identical service rates for all customer types. In the experiments below, the ratio

∑ n i=0 λ i N
is then hold constant. The results are presented in Table 2.11 and Figure 2.13(a). We observe that SP behaves much better than chaining as the number of skills increases. Figure 2.13(a) shows that for N = 10, the crossing value of t should be negative for high values of p (this means that SP is better in all cases). Single pooling behaves much better than chaining for the following two reasons. First as N increases, the flexibility in chaining decreases. A customer type in the chaining configuration has access to a fewer proportion of agent as N increases (the gap with the full flexible model increases). The second reason is related to the impact of the constant ratio

∑ n i=0 λ i N
, which increases the overall size of the call center as N increases. Having large call centers makes SP more efficient (see Section 2.5.5).

Table 2.11: Impact of the number of skills ( Constant Overall Workload. We again consider identical service rates for all customer types.

µ i = µ 0 = 0.2, W * 0 = W * i = 0.2, ∑ n i=0 λ i /N = 2, i = 1, ..., n, U = V = 1)
The summation ∑ n i=0 λ i is then hold constant. The results are presented in Table 2.12 and Figure 2.13(b). We distinguish two effects depending on p. For small values of p, the preference zone for SP reduces. The opposite is true for large values of p. The reason is related to the decreasing of the size of each team as N increases. Since we keep constant the overall workload, increasing the number of skills implies a lower demand per skill, which requires less agents per team. This makes the effect of pooling predominant. For the case of large p, the large number of customers 0 benefits from pooling under SP. For the case of small p, the system contains more regular customers, each of which benefits in chaining from the pooling of two adjacent teams.

Table 2.12: Impact of p, t and N on the staffing cost (

µ i = µ 0 = 0.2, W * 0 = W * i = 0.2, ∑ n i=0 λ i = 8, i = 1, ..., n, U = V = 1)

Mix of Asymmetry

In this section we mix the effects of more than a parameter at a time. We propose to interact the effects of p and p ′ , U and p, U and p ′ , U and V , and also all of them. The results are presented in Tables 2.13-2.17 and Figures 2.14(a)-2.14(d). From the numerical results, we observe that the individual effects are still present, but they may accumulate or make up for one another.

µi = µ0 = 0.2, W * 0 = W * i = 0.2, ∑ n i=0 λi/N = 2, i = 1, ..., 4, U = V = 1)
(b) Preference zone (µi = µ0 = 0.2, W * 0 = W * i = 0.2, ∑ n i=0 λi = 8, i = 1, ..., 4, U = V = 1)
Figure 2.13: Preference zone One important observation is that two asymmetries may lead to a bad performance for SP. For example SP behaves well in each one of the asymmetric situations (U = 2 and V = 1) and (U = 1 and V = 1/3) in isolation. However, it does not behave well for the mixed situation (U = 2 and V = 1/3). In such a situation, the customers types with large arrival rates are the faster to be served, and viceversa. Therefore, the different customer types workloads are likely to be symmetric. For the same reason, SP behaves well in the situation (U = 2 and V = 3) because the mix of asymmetries further accentuates the asymmetry in workloads.

Tables 2.14 and 2.15 reveal also that the most predominant effects are those of p (because of pooling) and p ′ (because of blocking). Various scenarios of mixed asymmetries are considered in Table 2.17. We find again that SP behaves well in large call centers (the first four scenarios).

Scenarios 3 and 7 are similar in terms of the values of p and p ′ (high values for the two parameters). This means that the effect of pooling and blocking are highly present in both scenarios.

An important observation here is that scenario 3 is the best among scenarios 1-4, while scenario 7

is the worst among scenarios 5-8. This gives an indication on the direct competition between the effects of p and p ′ . In large call centers, the pooling effect created by customers 0 is predominant over the blocking effect, and the opposite is true in small call centers.

Main Conclusions. In summary, the numerical analysis of this section confirms that single pooling Table 2.13: Impact of p and p ′ (µ i = µ j and λ i = λ j for i, j = 1, ..., 4, performs better than chaining for various cases of asymmetry in the system parameters. In the case of a predominance of customers 0 and/or an important asymmetry in the arrival rates of the regular types (captured by V ), SP is more robust than chaining even for small differences between the costs of a regular skill and that of skill 0. Because of the blocking effect, the performance of both chaining and SP deteriorates in the asymmetry defined by the service time duration of customers 0 relatively to that of regular customers. This is more apparent in single pooling because customers 0 have access to all teams, while in chaining they do only have access to two teams. We have also observed that SP is more robust than chaining against an increasing asymmetry between the service times of regular types. Since the teams under SP are less inter-dependent than under chaining, SP is again preferred in the case of an asymmetry between the objective service levels.

W * 0 = W * i = 0.2, ∑ 4 i=0 λ i = 8, ∑ 4 i=0 1 µ i = 25 i = 1, ..., 4, U = V = 1)
We therefore avoid over-staffing situations that may happen in chaining. Another important feature is that of abandonment, because it may affect the asymmetry of the parameters. Finally, all above conclusions are more apparent for the situations with a large number of skills, or for those with a large call center size.

i = λ j for i, j = 1, ..., 4, µ 0 = 0.2, W * 0 = W * i = 0.2, ∑ 4 i=0 λ i = 8, ∑ 4 i=0 1 µ i = 25 i = 1, ..., 4, p ′ = 20%, V = 1)

Concluding Remarks

We focused on a fundamental problem in the design and management of SBR call centers, for which it is important to choose an intelligent architecture. We considered the context of call centers with unbalanced workload, different service requirements, a predominant customer type and high costs of cross-training. With these asymmetry in the parameters, the well known existing architectures such as chaining lose their robustness. For those particular cases, we proposed a new call center architecture (single pooling) and demonstrated its efficiency. SP allows to balance the workload among the agents in a way that captures the benefits of pooling, without requiring every agent to process every type of call. The results of the comparison between SP and chaining have significant managerial implications. We showed that SP behaves well in most cases of asymmetry in the parameters. There might be then opportunities for managers of call centers to improve performance using the single pooling architecture.

= λ j = 1.5 for i, j = 1, ..., 4, λ 0 = 2, W * 0 = W * i = 0.2, ∑ 4 i=0 1 µ i = 25 i = 1, ..., 4, p = 25%, V = 1)
In a future research, it would be useful to extend the use of the fixed point algorithm to evaluate the performance measures of customers 0. Another interesting work is to generalize the functioning of single pooling in order to avoid the blocking effect in the case of long service times for customers 0.

Table 2.16: Impact of U and V (λ Optimal Email routing in a

0 = 0.8, µ 0 = 0.2, W * 0 = W * i = 0.2, ∑ 4 i=0 1 µ i = 25 and ∑ 4 i=0 λ i = 8, p = 10%, p ′ =
i, j = 1, ..., 4, W * 0 = W * i = 0.2, ∑ 4 i=0 λi = 8, ∑ 4 i=0 1 µ i = 25 i = 1, ..., 4, U = V = 1)
i, j = 1, ..., 4, µ0 = 0.2, W * 0 = W * i = 0.2, ∑ 4 i=0 λi = 8, ∑ 4 i=0 1 µ i = 25 i = 1, ..., 4, p ′ = 20%, V = 1)

Multi-Channel Call Center

Motivated by the call center practice, we consider a blended call center with calls arriving over time and an infinitely backlogged queue of emails. Calls have a non-preemptive priority over emails. The call service is characterized by three successive stages where the second one is a break, i.e., there is no required interaction between the customer and the agent for a non-negligible duration. This leads to a new opportunity to efficiently split the agent time between calls and emails.

We focus on the optimization of the email routing to agents. Our objective is to maximize the throughput of emails subject to a constraint on the call waiting time. We develop a general framework with two parameters for the email routing to agents. One parameter controls the routing between calls, and the other does the control inside a call. We then derive various structural results

with regard to the optimization problem and numerically illustrate them. Various guidelines to call center managers are provided. In particular, we prove for the optimal routing that all the time at least one of the two email routing parameters has an extreme value.

Introduction

Context and Motivation: Call centers are an important part of customers' service in many organizations. New technology-driven innovations are multiplying the opportunities to make more efficient use of an agent as she can handle different types of workflow, including inbound calls, outbound calls, email, chat, etc. However, several issues on the management of call center operations emerged also as a result of advanced technology. In this chapter, we consider a call center with two types of jobs, inbound calls and emails. We focus on how to efficiently share the agent time between the two types of jobs in order to improve the call center performance.

In practice, we encountered call center situations where inbound calls and emails are combined.

This is called blending. The key distinction of problems with blending comes from the fact that emails are less urgent and can be inventoried to some extent, relative to incoming calls. Therefore managers are likely to give a strict priority to calls over emails. An important question here is what should be the best way of routing of emails (or the non-urgent job) to agents, i.e., as a function of the systems parameters and the service level constraints (on calls and emails) when should we ask the agent to treat emails between the call conversations [START_REF] Bernett | Blended call center performance analysis[END_REF][START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF]Legros et al., 2013a). The email routing question is further important in the context of the call center applications we consider here. We encountered examples where a call conversation between an agent and a customer contains a natural break. We mean by this a time interval with no interaction between the agent and the customer. During the conversation, the agent asks the customer to do some necessary operations in her own (without the need of the agent availability).

After finishing those operations, the conversation between the two parties can start again. Inside an underway conversation, the agent is then free to do another task if needed. For an efficient use of the agent time, one would think about the routing of the less urgent jobs (emails) not only when the system is empty of calls, but also during call conversations. In practice, such a situation often occurs. For example, an agent in an internet hotline call center asks the customer to reboot her modem or her computer which may take some time where no interactions can take place. It is also often the case that a call center agent of an electricity supplier company asks the customer for the serial number of her electricity meter box. This box is usually located outside of the house and is locked, so, the customer needs some non-negligible time to get the required information. Another example is that of commercial call centers with a financial transaction during the call conversation.

After some time from the start of the call conversation, the customer is asked to do an online payment on a website before coming back to the same agent in order to finish the conversation.

The online paiement needs that the customer looks for her credit card, then she enters the credit card numbers, then she goes through the automated safety check with her bank (using SMS, etc.), which may take some minutes.

In the call center examples we encountered, the back office job could be a confirmation email of subscription or unsubscription, simple answers for various customer requests, etc. The answers to these emails usually consist on a set of preprepared text blocks that the agent should mix/adapt with the customer case. Some minutes are then enough for the agent to handle more than one email.

For such situations, it is natural that call center managers think about using the opportunity to route emails to an agent during the break of an undergoing call conversation, and not only when no calls are waiting in the queue. The main advantages are i) an efficient use of the agent time and therefore better call center performance, ii) also, agents become less bored because of the diversity of activities, and therefore they are kept from falling into a rut.

Main Contributions: In this chapter, we consider a call center with an infinite amount of outbound jobs (emails), and inbound jobs (calls) arriving over time for which a break is required in the middle of the call conversation. Given this type of call centers, we are interested in optimizing its functioning by controlling how the resource is shared between the two types of jobs. Calls are more important than emails in the sense that calls request a quasi-instantaneous answer (waiting time in the order of some minutes), however emails are more flexible and could be delayed for several tens of hours. An appropriate functioning is therefore that the agent works on calls as long as there is work to do for calls. The agent can then work on emails when she becomes free from calls, i.e., after a service completion when no calls are waiting in the queue, or during the call conversation break. We assume that calls have a non-preemptive priority over emails, which means that if a call is busy with an email (that has started after a service completion or during the break), the agent will finish first the email before turning to a new arrived call to the queue or a call that has accomplished the requested operations and wants to start again the conversation to finish her service. The non-preemption priority rule is coherent with the operations in practice and also to the call center literature [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF][START_REF] Deslauriers | Markov chain models of a telephone call center with call blending[END_REF]. It is not appropriate to stop the service of a low priority customer, and it is not efficient for the agent to stop the treatment of an email or to group emails for a simultaneous treatment. In Appendix in Section B.1 we prove that the simultaneous treatment is not an interesting opportunity in a call center.

We focus on the research question: when should the agent treat emails? Between calls, or inside a call conversation, or in both situations? Given the nature of the job types, a call center manager in practice would be interested in maximizing the number of treated emails while respecting some service level objective on the call waiting time [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF]. For calls, we are interested in the steady-state performance measures in terms of the expected waiting time, the probability that the waiting time is less than a given threshold, and the probability of delay. We do not consider the call waiting after the break because it is not perceived as badly as that before entering the first stage of service (less uncertainty for the customer because she has been already connected to an agent). For emails, we are interested in the steady-state performance in terms of the throughput of emails, i.e., the number of treated emails per unit of time.

The email routing problem considered in this chapter is a part of a collaboration between the authors and the French consulting company Interact-iv.com. In the small call center customers of Interact-iv.com, an efficient control of the agent time is important. For those call centers, Interactiv.com wants to implement, in the email dispatcher, an intelligent email routing algorithm adapted to the call system load.

Despite its prevalence, there are no papers in the call center literature addressing such a question. Most of the related papers only focus on the email routing between call conversations but not inside a call conversation. To answer this question, we develop a general framework with two parameters for the email routing to agents. One parameter controls the routing between calls, and the other does the control inside a call conversation. For the tractability of the analysis, we first focus on the single server case. We then discuss the extension of the results to the multi-server case.

For the single server modeling, we first evaluate the performance measures using a Markov chain analysis. Second, we propose an optimization method of the routing parameters for the problem of maximizing the email throughput under a constraint on the service level of the call waiting time.

As a function of the system parameters (the server utilization, the email service time, the severity of the call service level constraint, etc.), we derive various guidelines to managers. In particular, we prove for the optimal routing that all the time at least one of the two email routing parameters has an extreme value. As detailed later in this chapter, an extreme value means that the agent should do all the time emails inside a call (or between calls) or not at all. In other cases the parameters lead to randomized policies. We also solve our optimization problem by proposing 4 particular cases corresponding to the extreme values of the probabilistic parameters. We analytically derive the conditions under which one particular case would be preferred to another one. The interest from these particular cases is that they are easy to understand for agents and managers. Several numerical experiments are used to illustrate the analysis. To simplify the Markov chain analysis, we further propose an approximation method under the light-traffic regime.

The rest of the chapter is organized as follows. In Section 3.2 we review some of the related literature. In Section 3.3, we describe the blended call center modeling and the optimization problem. In Section 3.4, we develop a method based on the analysis of Markov chains in order to derive the performance measures of interest for calls and for emails. In Section 3.5, we focus on optimizing the email routing parameters. We also provide various numerical illustrations and discuss the results. In Section 3.6, we develop an approximation method for the system performance evaluation under the light-traffic regime. We also propose an approximation method to extend the results to the multi-server case. Finally in Section 3.7, we provide some concluding remarks and directions for future research.

Literature Review

There are three related streams of literature to this chapter. The first one deals with blended call centers. The second one is the Markov chain analysis for queueing systems with phase type service time distributions. The third one is related to the cognitive analysis, or in other words the ability for an agent to treat and switch between different job types.

The literature on blended call centers consists on developing performance evaluation and optimal blending policies. [START_REF] Deslauriers | Markov chain models of a telephone call center with call blending[END_REF] develop a Markov chain for the modeling of a Bell Canada blended call center with inbound and outbound calls. The performance measures of interest are the rate of outbound calls and the waiting time of inbound calls. Through simulation experiments they prove the efficiency of their Markov chain model to reflect reality. [START_REF] Brandt | A two-queue priority system with impatience and its application to a call center[END_REF] develop an approximation method to evaluate the performance of a call center model with impatient inbound calls and infinitely patient outbound calls of lower priority than the inbound traffic. [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF] consider a similar model to the one analyzed in this chapter, expect that the call service is done in a single stage without a possible break. The model consists on inbound and outbound jobs where the inbound jobs have a non-preemptive priority over the outbound ones. For the special case of identically distributed service times for the two jobs, they optimize the outbound jobs routing subject to a constraint on the expected waiting time of inbound jobs. Gans and Zhou (2003a) study a call center with two job types where one of the jobs is an infinitely backlogged queue. They develop a routing policy consisting on the reservation of servers in order to maximize the throughput on the jobs of the infinitely backlogged queue. Armony and Maglaras (2004a) analyze a similar model with a callback option for incoming customers. The customer behavior is captured through a probabilistic choice model. Other references include [START_REF] Bernett | Blended call center performance analysis[END_REF][START_REF] Keblis | Improving customer service operations at amazon.com[END_REF][START_REF] Pichitlamken | Modeling and simulation of a telephone call center[END_REF].

The analysis in this chapter is also related to the analysis of queueing systems with phase type service time distributions. We model the call service time through three successive exponentially distributed stages, where the second stage may also overlap with the service of one or several emails with an exponential time duration for each. The performance evaluation of such systems involves the stationary analysis of Markov chains and is usually addressed using numerical methods. We refer the reader to [START_REF] Kleinrock | Queueing Systems, Theory, volume I[END_REF] for simple models with Erlang service time distributions. For more complex systems, see [START_REF] Neuts | Explicit Steady-State Solutions to Some Elementary Queueing Models[END_REF]; [START_REF] Sze | A queueing model for telephone operator staffing[END_REF]; [START_REF] Bolotin | Telephone circuit holding time distributions[END_REF]; [START_REF] Brown | Statistical Analysis of a Telephone Call Center: A Queueing-Science Perspective[END_REF]; [START_REF] Guo | The Effects of Information on a Queue with Balking and Phase-Type Service Times[END_REF]. Our approach to derive the performance measures is based on first deriving the stationary system state probabilities for two-dimension and semi-infinite continuous time Markov chains. One may find in the literature three methods for solving such models. The first one is to truncate the state space, see for example [START_REF] Seelen | An Algorithm for Ph/Ph/c Queues[END_REF] and [START_REF] Keilson | Row-Continuous Finite Markov Chains: Structure and Algorithms[END_REF]. The second method is called spectral expansion [START_REF] Daigle | Queueing systems having phase-dependant arrival and service rates[END_REF][START_REF] Mitrani | Spectral Expansian Solution of a Class of Markov Models: Application and Comparison with the Matrix-Geometric Method[END_REF][START_REF] Choudhury | Numerical Solution of Mt/Gt/1 Queues[END_REF]. It is based on expressing the invariant vector of the process in terms of the eigenvalues and the eigenvectors of a matrix polynomial. The third one is the matrix-geometric method, see [START_REF] Neuts | Matrix-Geometric Solutions in Stochastic Models: an Algorithmic Approach[END_REF]. The approach relies on determining the minimal positive solution of a non-linear matrix equation. The invariant vector is then expressed in terms of powers of itself. In our analysis, we reduce the problem to solving cubic and quartic equations, for which we use the method of Cardan and Ferrari [START_REF] Gourdon | Les maths en tête : Algèbre[END_REF].

Finally, we briefly mention some studies on human multi-tasking, as it is the case for the agents in our setting. [START_REF] Gladstones | Division of Attention: The Single-Channel Hypothesis Revisited[END_REF] show that a simultaneous treatment of tasks is not efficient even with two easy tasks because of the possible interferences. In our models, we are not considering successive tasks in the sense that an agent can not talk to a customer and at the same time treats an email. More interestingly, [START_REF] Charron | Divided Representation of Concurrent Goals in the Huma Frontal Lobes[END_REF] studied the capacity of the frontal lobe to deal with different tasks by alternation (as here for calls and emails). They develop the notion of branching: capacity of the brain to remember information while doing something else. They show that the number of tasks done alternatively has to be limited to two to avoid loss of information. [START_REF] Dux | Training improves Multitasking Performance by Increasing the Speed of Information Precessing in Human Prefontal Cortex[END_REF] showed that training and experience can improve multi-tasking performance.

The risk from alternating between two tasks is the loss of efficiency because of switching times.

An important aspect to avoid inefficiency as pointed out by [START_REF] Dux | Training improves Multitasking Performance by Increasing the Speed of Information Precessing in Human Prefontal Cortex[END_REF] and [START_REF] Charron | Divided Representation of Concurrent Goals in the Huma Frontal Lobes[END_REF] is that the alternation should be at most between two tasks quite different in nature (like inbound and outbound jobs).

Problem Description and Modeling

For tractability and a better understanding of the results, we first focus on a single server queueing model. We then extend the analysis to the multi-server case in Section 3.6.2 using simulation.

We consider a single server queue with two types of jobs: calls and emails. The arrival process of customers is assumed to be Poisson with mean arrival rate λ. We assume to have an infinite amount of emails that are waiting to be treated in a dedicated first come, first served (FCFS) queue with an infinite capacity. One can think of a call center that stores a sufficiently large number of emails of a given day and handles them the next one.

Upon arrival, a call is immediately handled by the agent, if available. If not, the call waits for service in an infinite FCFS dedicated queue. Calls have a non-preemptive priority over emails. This means that the idle agent deals with a call first (the first in line). If the queue of calls is empty, this agent can handle an email (the first in line) from the queue of emails. Non-preemption priority is a natural assumption for our application. An agents in practice prefers to finish answering an underway email rather than starting it over later on. This is also preferred from an efficiency perspective. We assume in our models that there is no call abandonment or retrial.

As mentioned in Section 3.1, we consider call center applications where the communication between the agent and the customer includes a break (the customer does not need the agent availability). We model the service time of a call by 3 successive stages. The first stage is a conversation between the two parties. The second stage is the break, i.e., no interactions between the two parties. The third and final step is a again a conversation between the two parties. The service completion occurs as soon as the third stage finishes. We model each stage duration as an exponentially distributed random variable. The service rates of the first, second and third stages are denoted by µ 1 , µ 2 and µ 3 , respectively. This Markovian assumption, which is common in modeling in service operations, is reasonable for systems with high service time variability where service times are typically small but there are occasionally long service times. An agent handle an email within one single step without interruption. The time duration of an email treatment is random and assumed to be exponentially distributed with rate µ 0 .

We are interested in an efficient use of the agent time between calls and emails. More concretely,

we want to answer the question when should we treat emails for the following optimization problem

  
Maximize the throughput of emails subject to a service level constraint on the call waiting time in the queue.

(3.1)

To solve Problem (3.1), we propose a general model for the routing of emails to the agent.

It is refereed to as probabilistic model or Model PM and is described below. Recall that the call center has an infinite number of emails. Then defining the routing of emails consists of determining whether or not to start an email each time the agent becomes idle inside a call conversation, or after a call service completion with no waiting calls in the queue. Note that when the agent becomes idle during the second service stage of a currently served call, she cannot start to serve a new waiting call from the queue, if any. Such an overlap between two or more different call treatments would necessarily disturb the agent, which leads to errors and work inefficiency [START_REF] Gladstones | Division of Attention: The Single-Channel Hypothesis Revisited[END_REF].

Probabilistic Model (Model PM): We distinguish the two situations when the agent is available to handle emails between two call conversations, or inside a call conversation.

Between two calls: just after a call service completion (as soon as the third stage finishes) and no waiting calls are in the queue, the agent treats one or more emails with probability p (independently of any other event), or does not work on emails at all with probability 1 -p. In the latter case, the agent simply remains idle and waits for a new call arrival to handle it. In the former case (with probability p), she selects a first email to work on. After finishing the treatment of this email, there are two cases: either a new call has already arrived and it is now waiting in the queue, or the queue of calls is still empty. If a call has arrived, the agent handles that call. If not, she selects another email, and so on. At some point in time, a new call would arrive while the agent is working on an email. The agent will then handle the call as soon as she finishes the email treatment.

Inside a call : Just after the end of the first stage of an underway call service (regardless wether there are other waiting calls in the queue or not), the agent treats one or more emails with probability q (independently of any other event), or does not work on emails at all with probability 1 -q. In the latter case, the agent simply remain idle and waits for the currently served customer to finish her operations on her own (corresponding to the second call service stage, i.e., the agent break). As soon as the customer finishes on herself her second service stage, the agent starts the third and last service stage. In the former case (with probability q), she selects a first email to work on. After finishing the treatment of this email, there are two cases: either the currently served customer has already finished her second service stage, or not yet. If she does, the agent starts the third stage of the customer call service. If not, she selects another email, and so on. At some point in time, the Table 3.1: Particular cases of Model PM Model Description Model 1 p = q = 0, no treatment of emails Model 2 p = 1 and q = 0, systematic treatment of emails only between two calls Model 3 p = 0 and q = 1, systematic treatment of emails only during the break Model 4 p = q = 1, systematic treatment of emails between two calls and during the break currently served call would finish her second service while the agent is working on an email. The agent will then handle the call as soon as she finishes the the email treatment.

We further consider next 4 particular cases of Model PM as shown in Table 3.1. Although these models might appear to be too restrictive to solve Problem (3.1), we show later their merit in Section 3.5.2 when we focus on the optimization of p and q in Model PM. Moreover, they have the advantage of being easy to implement in practice, easy to understand by managers, and easy to follows by agents. Note that in Model 1, the throughput of emails is zero. The interest from Model 1 is in the extreme case of a very high workload of calls or a very restrictive constraint on the call waiting time.

Performance Analysis

In this section we provide an exact method to characterize the call waiting time in the queue and the email throughput for Model PM (Section 3.4.1) and its extreme cases (Section 3.4.2).

Our approach consists on using a Markov chain model to describe the system states and compute their steady-state probabilities. The computation of some of the steady-state probabilities involves the resolution of cubic (third degree) or quartic (fourth degree) equations for which we use the Cardan-Ferrari method.

Model PM

Let us define the random process {(x(t), y(t)), t ≥ 0} where x(t) and y(t) denote the state of the agent and the number of waiting calls in the queue at a given time t ≥ 0, respectively. We have y(t) ∈ {0, 1, 2, ...}, for t ≥ 0. The possible values of x(t) (corresponding to the possible states of the agent), for t ≥ 0, are -"Agent working on the first stage of a call service" denoted by x(t) = A, -"Idle agent that is waiting for the call to finish her second stage of service" denoted by x(t) = B, -"Agent working on an email while an underway call has already finished her second stage of service and is waiting for the agent to start her third stage of service" denoted by x(t) = B ′ , -"Agent working on the third stage of a call service" denoted by x(t) = C, -"Agent working on an email between two call conversations" denoted by x(t) = M , -"Agent idle between two call conversations" denoted by x(t) = 0.

Since call inter-arrival times, call service times in each stage, and email service times are exponentially distributed, {(x(t), y(t)), t ≥ 0} is a Markov chain (Figure 3.1).

M , 1 M , 2 M ,n+1 A , 0 B , 0 C , 0 C , 1 A , 1 B , 1 A , n B , n C , n A, n+1 B, n+1
C, n+1 and (M, n), respectively. We also define ρ i = λ µ i , for i ∈ {0, 1, 2, 3}. In Proposition 2, we give the probability of delay of a call (probability of waiting) denoted by P D and the throughput of emails denoted by T . Note that the stability condition of Model PM is λ <

3 pµ 3 µ 1 µ 2 qµ 0 µ B' , 0 B', 1 B' , n B', n+1 M , 0 M,n+2 0 , 0 λ (1 ) 3 p µ - λ (
q µ 0 + 1 µ 1 + 1 µ 2 + 1 µ 3 .
Proposition 2 For Model PM, we have

P D = 1 - 1 -p 1 + pρ 0 (1 -ρ 1 -ρ 2 -qρ 0 -ρ 3 ), T = µ 0 ( 1 + ρ 0 1 + pρ 0 p (1 -ρ 1 -ρ 2 -qρ 0 -ρ 3 ) + q(ρ 2 + ρ 0 )
) .

Proof. From the Markov chain of Model PM, we have

c 0 = ρ 3 (P 0 + m 0 ), c n = ρ 3 (a n-1 + b n-1 + b ′ n-1 + c n-1 + m n ), for n ≥ 1. Then ∞ ∑ n=0 c n = ρ 3 { P 0 + m 0 + ∞ ∑ n=0 ( a n + b n + b ′ n + c n ) + ∞ ∑ n=1 m n } . (3.2)
Since all system state probabilities sum up to 1, i.e.,

P 0 + ∑ ∞ n=0 (a n + b n + b ′ n + c n + m n ) = 1, Equation (3.2) becomes ∞ ∑ n=0 c n = ρ 3 . (3.3)
For the state (M, 0), we have pµ 3 c 0 = λm 0 , or equivalently c 0 = ρ 3 m 0 p . Therefore c 0 = ρ 3 P 0 1-p . We then may write

P 0 = 1 -p p m 0 . (3.4)
From the Markov chain, we also have

µ 2 ∑ ∞ n=0 b n = µ 3 ∑ ∞ n=0 c n = µ 0 ∑ ∞ n=0 b ′ n +(1-q)µ 2 ∑ ∞ n=0 b n = µ 1 ∑ ∞ n=0 a n . Using Equation (3.3), we then obtain ∞ ∑ n=0 b n = ρ 2 , ∞ ∑ n=0 b ′ n = qρ 0 , ∞ ∑ n=0 a n = ρ 1 . (3.5) For state (M, n), n ≥ 1, we have m n = ( ρ 0 1+ρ 0 ) n m 0 . Therefore ∞ ∑ i=0 m i = m 0 (1 + ρ 0 )
. Using now Equation (3.5) together with the normalization condition implies m 0 = p 1+pρ 0 (1-ρ 1 -ρ 2 -ρ 3 -qρ 0 ), and Equation (3.4) then becomes

P 0 = 1 -p 1 + pρ 0 (1 -ρ 1 -ρ 2 -qρ 0 -ρ 3 )).
A new call enters service immediately upon arrival, if and only if the system is in state (0, 0).

Since the call arrival process is Poison, we use the PASTA property to state that the steady-state probabilities seen by a new call arrival coincide with those seen at an arbitrary instant. Thus P D = 1 -P 0 , or

P D = 1 - 1 -p 1 + pρ 0 (1 -ρ 1 -ρ 2 -ρ 3 -qρ 0 ).
As for the email throughput, it is given by T = µ 0 (q

∑ ∞ i=0 b i + ∑ ∞ i=0 b ′ i + ∑ ∞ i=0 m i )
, which may be also written as

T = µ 0 ( 1+ρ 0 1+pρ 0 p (1 -ρ 1 -ρ 2 -ρ 3 -qρ 0 ) + q(ρ 2 + ρ 0 )
) . This finishes the proof of the proposition. ✷

Let us now define W , a random variable, as the call waiting time in the queue, and P(W < t)

as its cumulative distribution function (cdf) for t ≥ 0. Conditioning on a state seen by a new call arrival and averaging over all possibilities, we state using PASTA that

P(W < t) = P 0 • 1 + +∞ ∑ n=0 (P(W < t, (A, n)) • a n + P(W < t, (B, n)) • b n + P(W < t, (B ′ , n)) • b ′ n ) + P(W < t, (C, n)) • c n + P(W < t, (M, n)) • m n ). (3.6)
For n ≥ 0, the quantities P(W < t, (A, n)), P(W < t, (B, n)), P(W < t, (B ′ , n)), P(W < t, (C, n))

and P(W < t, (M, n)) are the cdf of the conditional call waiting times in the queue, given that a new arriving call finds the system in states (A, n), (B, n), (B ′ , n), (C, n) and (M, n), respectively.

In the Markov chain of Model PM, these conditional random variables correspond to first passage times to state (0, 0) starting from the system state upon a new call arrival. They are convolutions of independent exponential random variables with arbitrarily rates, not necessarily all equal (Erlang random variable) or all distinct. Using the results in [START_REF] Amari | Closed-form expressions for distribution of sum of exponential random variables[END_REF], we can explicitly derive the expressions of P(W < t, (A, n)), P(W < t, (B, n)), P(W < t, (B ′ , n)), P(W < t, (C, n))

and P(W < t, (M, n)), for n ≥ 0, as shown in Section B.2 of the Appendix.

It remains now to compute the probabilities a n , b n , b ′ n , c n and m n in n, for n ≥ 0. From the Markov chain of Model PM, we can write the following iterative equations

λX n-1 = AX n , (3.7)
for n ≥ 1, where

X n =            a n b n b ′ n c n m n            ,
for n ≥ 0 is the vector of probabilities to be computed and

A =            µ 1 -λ -λ -λ -λ -µ 1 λ + µ 2 0 0 0 0 -qµ 2 λ + µ 0 0 0 0 -(1 -q)µ 2 -µ 0 λ + µ 3 0 0 0 0 0 λ + µ 0           
.

The first step to solve Equation (3.7) is to find the eigenvalues of the matrix 1 λ A. These are solutions of the equation det( 1 λ A -yI) = 0 with y as variable. One obvious eigenvalue is 1 + 1 ρ 0

(see the last line of A), and the remaining ones are those of a 4 × 4 matrix (derived from 1 λ A by removing the last line and the last column) and they are solutions of the following quadric equation

σ 4 y 4 -(3σ 4 + σ 3 )y 3 + (3σ 4 + 2σ 3 + σ 2 )y 2 -(σ 4 + σ 3 + σ 2 + σ 1 )y + 1 + ρ 0 (1 -q) = 0, (3.8)
with y as variable,

σ 1 = ρ 0 + ρ 1 + ρ 2 + ρ 3 , σ 2 = ρ 0 ρ 1 + ρ 0 ρ 2 + ρ 0 ρ 3 + ρ 1 ρ 2 + ρ 1 ρ 3 + ρ 2 ρ 3 , σ 3 = ρ 0 ρ 1 ρ 2 + ρ 0 ρ 1 ρ 3 + ρ 0 ρ 2 ρ 3 + ρ 1 ρ 2 ρ 3 , and σ 4 = ρ 0 ρ 1 ρ 2 ρ 3 . Since the constant term 1 + ρ 0 (1 -q)
in Equation (3.8) is strictly positive, zero cannot be a solution of that equation. Then, 1 λ A is invertible. Therefore the eigenvalues of λA -1 are solutions of

(1 + ρ 0 (1 -q))x 4 -(σ 4 + σ 3 + σ 2 + σ 1 )x 3 + (3σ 4 + 2σ 3 + σ 2 )x 2 -(3σ 4 + σ 3 )x + σ 4 = 0, (3.9)
where x = 1 y . We solve the quadric Equation (3.9) using the Cardan-Ferrari method. In Section B.4 of the Appendix, we describe the details of this method.

The explicit expressions of the probability components of the vector X n , for n ≥ 0, can be derived, however they are too cumbersome for Model PM. We go further in providing their expressions for the extreme cases of Model PM in Section 3.4.2 and also using a light-traffic approximation in Section 3.6.1. In all cases, an exact numerical method is straightforward and easy to implement.

Numerical illustrations are shown later in Section 3.5.

Let us now compute the expected call waiting time in Model PM, denoted by E(W ). Consider first a model similar to Model PM except that emails can only be treated inside a call conversation.

We denote this model by Model PM'. With a little thought, one can see that the expected call waiting time in Model PM is that of Model PM' plus p µ 0 . This can be easily proven using the memoryless property of the email service duration and sample path arguments. The main idea of the proof is as follows. Consider each first call of the busy periods of Model PM'. The same calls arrive also at Model PM but not necessarily enter immediately service as in Model PM'. Each one of these calls in Model PM, will arrive either at a system that is empty of calls, or not. In the first case and with probability p, she will be delayed compared to Model PM' by the residual duration of an email treatment (exponential with rate µ 0 ). All the calls arriving after her (and seeing a system non-empty of calls) will be delayed by the same amount of time. In the second case, if the call arrives at a system non-empty of calls, then this means that the previous busy period of calls in Model PM has been delayed by an amount of time corresponding to the residual time of an email treatment. Then this call and all of those who arrive after her and see a system non-empty of calls will be delayed with the same amount of time, and so on and so forth.

Let us now compute the expected waiting time in Model PM', denoted by E(W ′ ). We use the Pollaczeck-Kinchin result for an M/G/1 queue. From [START_REF] Pollaczeck | Uber eine Aufgabe der Wahrscheinlichkeitstheorie[END_REF], we have

E(W ′ ) = ρ 2 (1+cv 2 ) 2λ(1-ρ) ,
where cv is the coefficient of variation of the service distribution (its standard deviation over its expected value) and ρ is the equivalent workload. Because of the possibility to do emails between calls, the random variable representing the service time duration, say S, can be written as

S = S 1 + S 2 + U S 0 + S 3
, where S i , a random variable, follows an exponential distribution with rate µ i , for i = 0, ..., 3, and U follows a binomial distribution with parameter q. We denote by E(Z) and V (Z) the expected (first moment) and variance of a given random variable Z. The first moment of S is given by

E(S) = 1 µ 1 + 1 µ 2 + q µ 0 + 1 µ 3 ,
and its variance can be written as (using the independence between S i and S j for i ̸ = j ∈ {0, ..., 3})

V (S) = V (S 1 ) + V (S 2 ) + V (U S 0 ) + V (S 3 ) = 1 µ 2 1 + 1 µ 2 2 + 2q -q 2 µ 2 0 + 1 µ 2 3 .
Then

cv 2 (S) = 1 µ 2 1 + 1 µ 2 2 + 2q-q 2 µ 2 0 + 1 µ 2 3 ( 1 µ 1 + 1 µ 2 + q µ 0 + 1 µ 3 ) 2 .
After some algebra, we obtain

E(W ′ ) = (ρ 1 + ρ 2 + ρ 3 ) 2 + ρ 2 1 + ρ 2 2 + ρ 2 3 + 2qρ 0 (ρ 0 + ρ 1 + ρ 2 + ρ 3 ) 2λ(1 -(ρ 1 + ρ 2 + qρ 0 + ρ 3 ))
, which leads to

E(W ) = p µ 0 + (ρ 1 + ρ 2 + ρ 3 ) 2 + ρ 2 1 + ρ 2 2 + ρ 2 3 + 2qρ 0 (ρ 0 + ρ 1 + ρ 2 + ρ 3 ) 2λ(1 -(ρ 1 + ρ 2 + qρ 0 + ρ 3 )) . (3.10)
This closes the performance measure analysis of Model PM.

Extreme Cases

We consider the 4 extreme cases of Model PM; Models 1,...,4. To derive the expressions of the email throughput, the call probability of delay, and the call expected waiting time, we simply apply the analysis of Section 3.4.1 and state the results as shown in Table 3.2.

One can derive the cdf of the call waiting time P(W < t) using Equation (3.6). For n ≥ 0,

T 0 µ 0 (1 -ρ 1 -ρ 2 -ρ 3 ) E(W ) (ρ 1 +ρ 2 +ρ 3 ) 2 (1+ ρ 2 1 +ρ 2 2 +ρ 2 3 (ρ 1 +ρ 2 +ρ 3 ) 2 ) 2λ(1-ρ 1 -ρ 2 -ρ 3 ) 1 µ 0 + E(W 1 ) P D ρ 1 + ρ 2 + ρ 3 1 Model 3 Model 4 T µ 0 (ρ 0 + ρ 2 ) µ 0 (1 -ρ 1 -ρ 3 ) E(W ) (ρ 1 +ρ 2 +ρ 3 +ρ 0 ) 2 (1+ ρ 2 1 +ρ 2 2 +ρ 2 3 +ρ 2 0 (ρ 1 +ρ 2 +ρ 3 +ρ 0 ) 2 ) 2λ(1-ρ 1 -ρ 2 -ρ 3 -ρ 0 ) 1 µ 0 + E(W 3 ) P D ρ 0 + ρ 1 + ρ 2 + ρ 3 1
the quantities P(W < t, (A, n)), P(W < t, (B, n)), P(W < t, (B ′ , n)), P(W < t, (C, n)) and P(W < t, (M, n)) can be derived using again the results of [START_REF] Amari | Closed-form expressions for distribution of sum of exponential random variables[END_REF] We then show that u n , for n ≥ 0, satisfies a recurrent cubic or quartic linear equation that we solve using the Cardan-Ferrari method.

Another method to compute P(W < t) is this case (service in three exponential stages) based on the Dunford decomposition is developed in Section B.3 of the Appendix.

Model 1: The Markov chain associated to Model 1 is given in Figure 3.2. For this model, we have

u n = a n + b n + c n , for n ≥ 0.
From the Markov chain of Model 1, we may write c n (λ + µ 3 ) = λc n-1 + µ 2 b n for n ≥ 1, and

c n = ρ 3 u n-1 for n ≥ 1. This leads to b n = ρ 2 ((1 + ρ 3 )u n-1 -ρ 3 u n-2 ) for n ≥ 2.
From the Markov chain, we also may write b n (λ + µ 2 ) = λb n-1 + µ 1 a n for n ≥ 1. Combining the last two equations yields to

a n = ρ 1 (((1 + ρ 3 )ρ 2 + ρ 3 + 1)u n-1 -(ρ 2 (1 + 2ρ 3 ) + ρ 3 )u n-2 + ρ 3 u n-3 ) for n ≥ 3. The 0 , 0 A , 0 B , 0 C , 0 C , 1 A , 1 B , 1 A , n B , n C , n A, n+1 B, n+1 C, n+1 3 µ 3 µ 3 µ 1 µ 1 µ 1 µ 1 µ 2 µ 2 µ 2 µ 2 µ λ λ λ λ λ λ λ 3 µ 3 µ Figure 3.2: Markov chain for Model 1 equation u n = a n + b n + c n for n ≥ 0 is then equivalent to u n+3 = ((1 + ρ 1 )(1 + ρ 2 )(1 + ρ 3 ) -1) u n+2 -(ρ 1 ρ 2 + ρ 1 ρ 3 + ρ 2 ρ 3 + 2ρ 1 ρ 2 ρ 3 ) u n+1 + ρ 1 ρ 2 ρ 3 u n ,
which leads to

u n+3 = (σ 1 + σ 2 + σ 3 ) u n+2 -(σ 2 + 2σ 3 ) u n+1 + σ 3 u n , (3.11)
where • If ∆ > 0, the cubic equation has one real solution denoted by x 1 and two complex solutions

σ 1 = ρ 1 + ρ 2 + ρ 3 , σ 2 = ρ 1 ρ 2 + ρ 1 ρ 3 + ρ 2 ρ 3 ,
x 2 and its conjugate. We denote by |x 2 | the module of x 2 and θ its argument. Since u n is real, it can be written as

u n = rx n 1 + s|x 2 | n cos(nθ) + t|x 2 | n sin(nθ),
for n ≥ 0, where the parameters r, s, t ∈ R. These parameters are easily computed using the initial conditions given by u 0 , u 1 and u 2 .

• If ∆ < 0, the cubic equation has three real solutions denoted by x 1 , x 2 and x 3 . Thus

u n = rx n 1 + sx n 2 + tx n 3 ,
for n ≥ 0, where again the parameters r, s and t are computed from the initial conditions.

• If ∆ = 0, the cubic equation has a simple solution denoted by x 1 and a double one denoted by x 2 . We then have

u n = rx n 1 + (sn + t)x n 2 ,
for n ≥ 0, and the real parameters r, s and t are again given by the initial conditions.

Model 2: The Markov chain associated to Model 2 is given in Figure 3.3. For this model, we have

u n = a n + b n + c n , for n ≥ 0.
Following the same reasoning as that for Model 1, we obtain

͕ Ŷ λ ͕ Ŷ ͕ Ŷнϭ 1 µ 1 µ 1 µ 2 µ λ 0 µ D ͕ Ϭ ͕ Ϭ ͕ Ϭ ͕ Ϭ ͕ ϭ ͕ ϭ ͕ ϭ 3 µ 3 µ 1 µ 2 µ 2 µ λ λ λ D ͕ ϭ D ͕ Ϯ 0 µ ͕ Ŷнϭ ͕ Ŷнϭ 3 µ 2 µ λ λ λ 0 µ 0 µ λ 3 µ D͕ ŶнϮ D ͕ Ŷнϭ 3 µ ͕ Ŷ ͕ Ŷ Figure 3.3: Markov chain for Model 2 u n+3 = [σ 1 + σ 2 + σ 3 ]u n+2 -[σ 2 + 2σ 3 ]u n+1 + σ 3 u n + K ( ρ 0 1 + ρ 0 ) n , (3.12)
with the constant K ∈ R. Since Equation (3.12) is similar to Equation (3.11) in the sense that the former has just an additional term proportional to

( ρ 0 1+ρ 0 ) n
, we again use the solutions of Equation

for n ≥ 4. From u n = a n + b n + b ′ n + c n we then state that u n+4 = (ρ 0 + ρ 1 + ρ 2 + ρ 3 + ρ 0 ρ 1 + ρ 0 ρ 2 + ρ 0 ρ 3 + ρ 1 ρ 2 + ρ 1 ρ 3 + ρ 2 ρ 3 + ρ 0 ρ 1 ρ 2 + ρ 0 ρ 1 ρ 3 + ρ 0 ρ 2 ρ 3 + ρ 1 ρ 2 ρ 3 + ρ 0 ρ 1 ρ 2 ρ 3 )u n+3 -(ρ 0 ρ 1 + ρ 0 ρ 2 + ρ 0 ρ 3 + ρ 1 ρ 2 + ρ 1 ρ 3 + ρ 2 ρ 3 + 3(ρ 0 ρ 1 ρ 2 + ρ 0 ρ 1 ρ 3 + ρ 0 ρ 2 ρ 3 + ρ 1 ρ 2 ρ 3 ) + 3ρ 0 ρ 1 ρ 2 ρ 3 )u n+2 + (ρ 0 ρ 1 ρ 2 + ρ 0 ρ 1 ρ 3 + ρ 0 ρ 2 ρ 3 + ρ 1 ρ 2 ρ 3 + 2ρ 0 ρ 1 ρ 2 ρ 3 )u n+1 -ρ 0 ρ 1 ρ 2 ρ 3 u n ,
for n ≥ 0, or equivalently

u n+4 = [σ 1 + σ 2 + σ 3 + σ 4 ]u n+3 -[σ 2 + 2σ 3 + 3σ 4 ]u n+2 + [σ 3 + 3σ 4 ]u n+1 -σ 4 u n , (3.13)
for n ≥ 0, where

σ 1 = ρ 0 + ρ 1 + ρ 2 + ρ 3 , σ 2 = ρ 0 ρ 1 + ρ 0 ρ 2 + ρ 0 ρ 3 + ρ 1 ρ 2 + ρ 1 ρ 3 + ρ 2 ρ 3 , σ 3 = ρ 0 ρ 1 ρ 2 + ρ 0 ρ 1 ρ 3 + ρ 0 ρ 2 ρ 3 + ρ 1 ρ 2 ρ 3 and σ 4 = ρ 0 ρ 1 ρ 2 ρ 3 .
In order to compute u n we then solve the associated quadric equation using the Cardan-Ferrari method and the initial conditions u 0 , u 1 , u 2 and u 3 .

Model 4: The Markov chain associated to Model 4 is given in Figure 3.5. For this model, we have

u n = a n + b n + b ′ n + c n , for n ≥ 0. M , 1 M , 2 M ,n+1 λ λ A , 0 B , 0 C , 0 C , 1 A , 1 B , 1 A , n B , n C , n A, n+1 B, n+1 C, n+1
Following the same reasoning as that for Model 3, we obtain

u n+4 = [σ 1 + σ 2 + σ 3 + σ 4 ] u n+3 -[σ 2 + 2σ 3 + 3σ 4 ]u n+2 + [σ 3 + 3σ 4 ]u n+1 -σ 4 u n + K ( ρ 0 1 + ρ 0 ) n , (3.14)
with the constant K ∈ R. Since Equation (3.14) is similar to Equation (3.13) in the sense that the former has just an additional term proportional to

( ρ 0 1+ρ 0 ) n
, we again use the solution of Equation (3.13) to easily obtain those of Equation (3.14). This closes the characterization of the performance measures for the special cases Models 1,...,4.

Comparison Analysis and Insights

We start in Section 3.5.1 by a comparison analysis between the extreme cases Models 1,...,4. The comparison is based on the optimization problem (3.1). We derive various structural results and properties for this comparison. In particular, we investigate the impact of the mean arrival rate intensity of calls on the comparison between Models 1,...,4. One could think of a call center manager that adjusts the job routing schema as a function of the call arriving workload over the day. In Section 3.5.2 we focus on the general case Model PM. We prove that the optimization of the parameters of Model PM lead to extreme situations in the sense of a systematic email treatment of emails either between calls or inside a call conversation, which gives an interest in practice for Models 1,...,4.

Comparison between the Extreme Cases

We first compare between Models 1,...,4 based on their performance in terms of the email throughput, denoted by T 1 ,..., T 4 , respectively. It is obvious that Model 4 is the best and Model 1 is the worst (no emails at all). Let us now compare between Models 2 and 3. From Table 3.2 we have

T 2 = µ 0 (1 -ρ 1 -ρ 2 -ρ 3 ) and T 3 = µ 0 (ρ 0 + ρ 2 ). Thus T 3 > T 2 is equivalent to λ > 1 1 µ + 1 µ 2
,

where 1 µ = 1 µ 0 + 1 µ 1 + 1 µ 2 + 1 µ 3 . Since the stability condition for Model 3 is λ < µ, Model 3 is better than Model 2 if 1 1 µ + 1 µ 2 < λ < µ. (3.15)
Denoting the left term in Inequality (3.15) by R, the condition under which T 3 > T 2 is then

R = 1 1 µ 0 + 1 µ 1 + 1 µ 3 + 2 µ 2 < λ < µ. (3.16)
From Inequality (3.16), we first see that treating emails only inside a call conversation (Model 3) becomes better that treating them only between calls (Model 2) is likely the case for high arrival workloads (in such a case, busy period durations are reduced). We also see that ∂R ∂µ 2 > 0 for µ 2 > 0, ∂R ∂µ 0 > 0 for µ 0 > 0, ∂R ∂µ 1 > 0 for µ 1 > 0, and ∂R ∂µ 3 > 0 for µ 3 > 0. This means that decreasing the expected duration of the call service second stage (1/µ 2 ) relative to the expected durations of the other call service stages or the email service duration (1/µ 1 , 1/µ 3 and 1/µ 0 ) increases the range of arrival workloads where it is preferred to use Model 2 instead of Model 3. In other word, there is no sufficient time to treat emails inside the call conversation.

Comparison with a constraint on E(W )

As a function of the mean call arrival rate, we want to answer the question when should we treat emails (which model among Models 1 to 4 should a manger choose?) for the following problem

   Maximize T subject to E(W ) ≤ w * , (3.17)
where w * is the service level for the expected waiting time, w * > 0. Let W i , a random variable, denote the expected call waiting time in Model i, i = 1, ..., 4. It is clear that for some periods of a working day with a very high call arrival rate λ, the manager is likely to choose Model 1 (no emails), and for other periods with a very low λ, she is likely to choose Model 4 (emails between calls and inside a call). However for intermediate values of λ, the optimal choice is not clear. This is what we analytically analyze in what follows.

Under the condition of stability of Model i, E(W i ) is continuous and strictly increasing in λ (see

Table 3.1), for i = 1, ..., 4. The constraint E(W i ) ≤ w * is then equivalent to λ ≤ λ i , for i = 1, ..., 4,
where

λ 1 = 2w * 2w * ( 3 ∑ i=1 1 µ i ) + ( 3 ∑ i=1 1 µ i ) 2 + 3 ∑ i=1 1 µ 2 i , (3.18) λ 2 = 2 ( w * -1 µ 0 ) 2 ( w * -1 µ 0 ) ( 3 ∑ i=1 1 µ i ) + ( 3 ∑ i=1 1 µ i ) 2 + 3 ∑ i=1 1 µ 2 i , λ 3 = 2w * 2w * ( 3 ∑ i=0 1 µ i ) + ( 3 ∑ i=0 1 µ i ) 2 + 3 ∑ i=0 1 µ 2 i , λ 4 = 2 ( w * -1 µ 0 ) 2 ( w * -1 µ 0 ) ( 3 ∑ i=0 1 µ i ) + ( 3 ∑ i=0 1 µ i ) 2 + 3 ∑ i=0 1 µ 2 i .
For a given λ and under the condition of stability of Model i (i = 1, ..., 4), the choice of Model i happens if λ ≤ λ i and T i = max j∈{1,...,4}, λ≤λ j (T j ). When λ ≤ λ 4 , the choice is obviously for Model 4. When λ ≤ λ 1 and λ > λ i for i = 2, 3, 4 the only possibility is Model 1. Proposition 3 provides the conditions under which an optimal choice of Model 2 or Model 3 may happen.

Proposition 3 The following holds:

1. For λ < 1 µ 1 + 1 µ 2 + 1 µ 3 , it exists some values of λ for which Model 2 is optimal if and only if

λ 2 > 0. 2. For λ < 1 µ 0 + 1 µ 1 + 1 µ 2 + 1 µ 3 , it exists some values of λ for which Model 3 is optimal if and only if                      R ≤ λ 3 or λ 2 < λ 3 .
3. We have λ 2 < λ 3 if and only if 1 µ 0 < w * < w * , where

w * = 1 2 4 µ 2 0 + ( 1 µ 3 + 1 µ 2 + 1 µ 1 ) 4 µ 0 + 5 3 ∑ i=1 1 µ 2 i + 6 3 ∑ i,j=1;i̸ =j 1 µ i µ j - ( 1 µ 3 + 1 µ 2 + 1 µ 1
) .

Proof. See Section C.2 of the appendix. ✷

Using Equation (3.18), the condition in the first statement of Proposition 3 may be rewritten

as                        w * > 1 µ 0 or w * < 1 µ 0 - ( 3 ∑ i=1 1 µ i ) 2 + 3 ∑ i=1 1 µ 2 i 2 ( 3 ∑ i=1 1 µ i
) .

(3.19)

The second inequality in Relation (3.19) implies w * < 1 µ 0 . Since at least the expected waiting time in Model 2 is strictly higher than 1 µ 0 (any new call has at least to wait for the residual time of an email treatment), this second inequality is impossible. Roughly speaking, the condition for the optimality of Model 2 (for some values of λ) holds when the service level on the call waiting is higher than the expected email service time.

In what follows, we numerically illustrate the analysis above. For various system parameters, Figure 3.6 gives the optimal model choice as a function of the mean arrival rate of calls, λ. An intuitive reasoning of a manager would choose the ordering Model 4 (emails between calls and inside a call), then 2 (emails only between calls), then 3 (emails only inside a call), then 1 (no emails) as λ increases. The ordering Model 2 then Model 3 is not always appropriate, and some situations may require to consider some counterintuitive ordering. For instance, Model 3 is better than Model 2 for small values of λ if R ≤ λ 4 and λ 3 < λ 2 , see Figure 3.6(c). In other words, this happens when the constraint on E(W ) is not too restrictive and when the expected second stage service duration is long. Another more surprising ordering, as λ increases, is Model 2, then Model 3, then again Model 2 (see Figure 3.6(d)) which happens for system parameters such that λ 4 < R < λ 3 < λ 2 .

Comparison with a constraint on P (W < AW T )

In the constraint of Problem (3.1), we want that the probability that a call waits less than a given threshold, defined as AW T is at least a given service level, defined as SL, i.e., P(W < AW T ) ≥ SL.

Note that a special case of this constraint is that on P D , the call probability of waiting. The expressions involved in the analysis of P (W < AW T ) are quite complicate to allow an analytical comparison between the models as we have done for a constraint on E(W ). We have then conducted a numerical comparison analysis (not totaly illustrated here). The main qualitative conclusions are similar to those for the case of a constraint on E(W ). As λ increases, it is not always true as one would intuitively expect that a manager should choose first Model 2 and then at some point of λ she shifts to Model 3 (Figure 3.7(a)). The optimal choice may change with the system parameters and we may have the ordering Model 3 then Model 2 (Figure 3.7(b)).

Optimization of Model PM

In this section we focus on the general case, Model PM. We are interested in the optimization of the parameters p and q in Model PM for Problem (3.1). Concretely, we want to find the optimal routing parameters of Model PM that allows the manager to maximize the email throughput while respecting a call service level constraint. Recall that the stability condition of Model PM is

λ < q µ 0 + 1 µ 1 + 1 µ 2 + 1 µ 3 .
The expression of the email throughput T for Model PM is given in Proposition 2. It is straightforward to prove that for p, q ∈ [0, 1] the maximum of T (best situation) is reached for p = q = 1. The proof is then omitted. Also, the expected call waiting time of Model PM (given in Equation (3.10)) is maximized (worst) for p = q = 1. Therefore in order to solve Problem (3.1), one would be interested analyzing the sensitivity of T with respect to p and q. In Lemma 1 we prove a sensitivity result for T . The result will be used later in our analysis.

Lemma 1 We have ∂T ∂p > ∂T ∂q if and only if 0 < ρ 0 < ρ 0 , where

ρ 0 = √ (p -q -ρ 1 -ρ 2 -ρ 3 ) 2 -4(p 2 -p -q)(1 -ρ 1 -ρ 2 -ρ 3 ) -q + p -(ρ 1 + ρ 2 + ρ 3 ) 2(q -p 2 + p) .
Proof. See Section B.6 of the appendix.

✷

In what follows we address the question: starting from p = q = 1, in which direction should we decrease T ? Should we decrease p or q first? For p = q = 1, we have

ρ 0 = 1 2 ( √ (ρ 1 + ρ 2 + ρ 3 ) 2 + 4(1 -ρ 1 -ρ 2 -ρ 3 ) -(ρ 1 + ρ 2 + ρ 3 )
) . Let us now prove (for p = q = 1) that ρ 0 > ρ 0 . From the one hand, proving ρ 0 > ρ 0 is equivalent

to proving √ (ρ 1 + ρ 2 + ρ 3 ) 2 + 4(1 -ρ 1 -ρ 2 -ρ 3 ) > 2ρ 0 + (ρ 1 + ρ 2 + ρ 3 ) or equivalently ρ 2 0 + ρ 0 (ρ 1 + ρ 2 + ρ 3 ) -(1 -(ρ 1 + ρ 2 + ρ 3 )) < 0 or also (ρ 0 + 1)(ρ 0 -(1 -(ρ 1 + ρ 2 + ρ 3 )) < 0.
From the other hand, we have ρ 0 + 1 > 0. Also, the stability condition of Model 4 (Model PM with

p = q = 1) is ρ 0 + ρ 1 + ρ 2 + ρ 3 < 1. Then ρ 0 < 1 -(ρ 1 + ρ 2 + ρ 3 ).
As a conclusion the inequality ρ 0 > ρ 0 is true, for p = q = 1. Using Lemma 1, this means that starting from p = q = 1, we have ∂T ∂p > ∂T ∂q > 0. As a consequence, when we need to modify the values of p and q in order to decrease the expected call waiting time (the constraint in Problem (3.1)), the maximum of T is guarantied by first decreasing q (the email throughput is less sensitive to the variation of q than that of p). The question now is: which direction to use next? In other words when p = 1 and some value of q such that 0 < q < 1, is it possible that it is better to decrease p instead of q? The answer is no and the proof is as follows. For p = 1, let us try to find a value of q for which we have ρ 0 ≤ ρ 0 . This is equivalent to

√ (1-q-ρ 1 -ρ 2 -ρ 3 ) 2 +4q(1-ρ 1 -ρ 2 -ρ 3 )-q+1-ρ 1 -ρ 2 -ρ 3 2q ≤ ρ 0 . Thus, q 2 ρ 2 0 + qρ 0 -(1 -ρ 1 -ρ 2 -ρ 3 )(1 + ρ 0 ) > 0.
This trinomial in q has two real solutions;

q 1 = - 1+ √ 4ρ 0 +5-4(ρ 1 +ρ 2 +ρ 3 )(ρ 0 +1) 2ρ 0 and q 2 = -1+ √ 4ρ 0 +5-4(ρ 1 +ρ 2 +ρ 3 )(ρ 0 +1) 2ρ 0
. Obviously q 1 < 0. We also have q 2 > 1 because: proving q 2 -1 > 0 is equivalent to proving ρ 2 0 + (ρ 1 + ρ 2 + ρ 3 )ρ 0 + 1 > 0.

The discriminant of this latter trinomial in ρ 0 is negative and it is equal to (ρ 1 + ρ 2 + ρ 3 ) 2 -4. So q 2 > 1 for any ρ 0 > 0. Therefore it is impossible to find a value of q between 0 and 1 for which 0 < ∂T ∂p < ∂T ∂q . In conclusion starting from p = q = 1, when we need to change the values of p and q, the best direction to maximize T is to first decrease q until q = 0 and only then start to decrease p from p = 1.

Consider now Problem (3.1) with a constraint on E(W ). From the one hand, the constraint

E(W ) ≤ w * implies p µ 0 + (ρ 1 + ρ 2 + ρ 3 ) 2 + ρ 2 1 + ρ 2 2 + ρ 2 3 + 2qρ 0 (ρ 0 + ρ 1 + ρ 2 + ρ 3 ) 2λ(1 -(ρ 1 + ρ 2 + qρ 0 + ρ 3 )) ≤ w * , for p, q ∈ [0, 1], or equivalently q ≤ 2λ(1 -ρ 1 -ρ 2 -ρ 3 )(w * -p/µ 0 ) -(ρ 1 + ρ 2 + ρ 3 ) 2 -ρ 2 1 -ρ 2 2 -ρ 2 3 2ρ 0 (ρ 0 + ρ 1 + ρ 2 + ρ 3 + λ(w * -p/µ 0 )) , (3.20)
for p, q ∈ [0, 1]. From the other hand, the condition in Lemma 1, 0 < ρ 0 < ρ 0 , is equivalent to In what follows we prove an interesting result on the optimal values of p and q. Consider for example Figure 3.8(a) and assume that the agent is in a situation such that (p, q) belongs to Zone 1 or 2. Then the constraint on E(W ) is respected, but T can be improved. Using Lemma 1, we should increase p first (q first) for Zone 1 (Zone 2). From Figure 3.8(a), we also see that we should decrease p first (q first) for Zone 3 (Zone 4). It is clear that the optimal couple (p, q) will be on the curve of f . Moreover, we prove in Theorem 1 that the optimal point is such that p ∈ {0, 1} or q ∈ {0, 1}.

q < 1 -(ρ 1 + ρ 2 + ρ 3 )(1 + ρ 0 ) ρ 0 (1 + ρ 0 ) + 1 -ρ 0 1 + ρ 0 p + ρ 0 1 + ρ 0 p 2 , ( 3
Theorem 1 For p, q ∈ [0, 1], the optimal values of p and q of the optimization problem

   Maximize T subject to E(W ) ≤ w * , (3.22)
are always extreme values (0 or 1) for at least p or q.

Proof. We want to maximize the email throughput T (p, q) while respecting a constraint on the expected call waiting time (E(W )(p, q) ≤ w * ). We use the method of Lagrange multipliers to find the optimal point (p, q). Let us denote the Lagrange multiplier by α (α is real). Then α and the extremum (p, q) of our optimization problem are solutions of the set of the 3 equations D(T + α(W -w * )) = 0, where D is the differential applicator in α, p and q. These 3 equations are

∂(T + α(W -w * )) ∂p = µ 0 (1 -ρ 1 -ρ 2 -qρ 0 -ρ 3 )(1 + ρ 0 ) (1 + pρ 0 ) 2 + α 1 µ 0 = 0, ∂(T + α(W -w * )) ∂q = µ 0 (1 -p)ρ 0 1 + pρ 0 + α 1 2λ ρ 0 (2(ρ 0 + ρ 1 + ρ 2 + ρ 3 )(1 -(ρ 1 + ρ 2 + ρ 3 )) + (ρ 1 + ρ 2 + ρ 3 ) 2 + ρ 2 1 + ρ 2 2 + ρ 2 3 (1 -(ρ 1 + ρ 2 + ρ 3 + qρ 0 )) 2 = 0, ∂(T + α(W -w * )) ∂α = p µ 0 + (ρ 1 + ρ 2 + ρ 3 ) 2 + ρ 2 1 + ρ 2 2 + ρ 2 3 + 2qρ 0 (ρ 0 + ρ 1 + ρ 2 + ρ 3 ) 2λ(1 -(ρ 1 + ρ 2 + qρ 0 + ρ 3 )) -w * = 0.
Solving this set of 3 equations leads to two couples of solutions (p 1 , q 1 ) and (p 2 , q 2 ). The expressions of these solutions are too long (see for example Section B.7 of the appendix for the expression of q 2 ), but we show for any case of system parameters ρ i (i = 0, ..., 3) under the condition of stability that all the values of p 1 , q 1 , p 2 and q 2 are out of the interval [0, 1]. We then deduce for the optimal values of p and q that at least one them has an extremum value (0 or 1). This finishes the proof of the theorem. ✷ .9 provides a numerical illustration of Theorem 1. We observe as a function of the system parameters that at least one of the routing parameters is either 0 or 1. This gives the merit to the study of the extreme cases Models 1,...,4. Moreover they are easy to implement and easy to understand for both managers and agents. Note that the same observation holds also for Problem (3.1) with a constraint on P (W < AW T ). This is based on several numerical examples (a rigorous proof as that for E(W ) is too complex to develop). An illustration is given in Figure 3.10.

Approximations

In this section, we develop two approximation methods. In Section 3.6.1 we focus on a lighttraffic approximation for Model PM in the case of a single server. In Section 3.6.2, we propose an approximation method to extend our email routing optimization problem to the multi-server case.

Light-Traffic Approximation for Model PM

Although the analysis in Section 3.4.1 provides an exact method to derive the system steady-state probabilities of Model PM, the expressions involved are cumbersome. In this section, we develop an approximate analysis that allows to obtain simpler expressions for those quantities under the light-traffic regime.

Under the light-traffic regime, we have ρ i << 1 and also ρ i 1+ρ i ≃ ρ i , for i = 0, ..., 3. Equation (3.9) can be then rewritten as

x 4 -σ 1 x 3 + σ 2 x 2 -σ 3 x + σ 4 = 0, or equivalently (x -ρ 0 )(x -ρ 1 )(x -ρ 2 )(x -ρ 3 ) = 0.
Also, the matrix 1 λ A becomes

1 λ A =            1/ρ 1 0 0 0 0 -1/ρ 1 1/ρ 2 0 0 0 0 -q/ρ 2 1/ρ 0 0 0 0 -(1 -q)/ρ 2 -1/ρ 0 1/ρ 3 0 0 0 0 0 1/ρ 0           
, which easily leads to

λA -1 =            ρ 1 0 0 0 0 ρ 2 ρ 2 0 0 0 qρ 0 qρ 0 ρ 0 0 0 ρ 3 ρ 3 ρ 3 ρ 3 0 0 0 0 0 ρ 0           
.

The computation of the steady-state probabilities simplifies as follows. We have

P 0 = 1-p, m 0 = p, a 0 = ρ 1 , b 0 = ρ 2 , b ′ 0 = qρ 0 and c 0 = ρ 3 . Using Equation (3.7), we obtain a n = ρ n+1 1 , for n ≥ 0. Also, b n = ρ 2 a n-1 + ρ 2 b n-1 for n ≥ 1. Therefore, if µ 1 ̸ = µ 2 , we have b n = ρ 2 n ∑ k=0 ρ k 1 ρ n-k 2 = ρ 2 ρ n+1 2 -ρ n+1 1 ρ 2 -ρ 1
for n ≥ 0, and if µ 1 = µ 2 , we have b n = (n + 1)ρ n+1 1 for n ≥ 0. From Equation (3.7), we may write b

′ n = qρ 0 a n-1 + qρ 0 b n-1 + ρ 0 b ′ n-1 , for n ≥ 1. Thus if µ 1 ̸ = µ 2 , we have b ′ n = qρ 0 ρ n+1 2 -ρ n+1 1 ρ 2 -ρ 1 + ρ 0 b ′ n-1 , for n ≥ 1. Therefore if µ i ̸ = µ j for i ̸ = j ∈ {0, 1, 2} we have b ′ n = qρ 0 ρ 2 (ρ 2 -ρ 1 )(ρ 2 -ρ 0 ) (ρ n+1 2 -ρ n+1 0 ) - qρ 0 ρ 1 (ρ 2 -ρ 1 )(ρ 1 -ρ 0 ) (ρ n+1 1 -ρ n+1 0 ), or equivalently b ′ n = qρ 0 ( ρ n+2 0 (ρ 0 -ρ 1 )(ρ 0 -ρ 2 ) + ρ n+2 1 (ρ 1 -ρ 0 )(ρ 1 -ρ 2 ) + ρ n+2 2 (ρ 2 -ρ 0 )(ρ 2 -ρ 1 )
) ,

for n ≥ 0. If µ i = µ j and µ i ̸ = µ k for distinctive i, j, k ∈ {0, 1, 2}, then b ′ n = qρ 0 ( (n + 2)ρ n+1 i ρ i -ρ k + ρ n+2 k (ρ k -ρ i ) 2 ) , for n ≥ 0. If µ 1 = µ 2 = µ 0 , we have b ′ n = qρ n+1 0 (n + 1)(n + 2) 2 ,
for n ≥ 0. From Equation (3.7), we also may write

c n = ρ 3 (a n-1 + b n-1 + b ′ n-1 + c n-1 ), for n ≥ 1.
Thus, if µ i ̸ = µ j for i ̸ = j ∈ {0, 1, 2, 3}, we have

c n =ρ 3 ( ρ n+2 0 (ρ 0 -ρ 0 (1 -q)) (ρ 0 -ρ 1 )(ρ 0 -ρ 2 )(ρ 0 -ρ 3 ) + ρ n+2 1 (ρ 1 -ρ 0 (1 -q)) (ρ 1 -ρ 0 )(ρ 1 -ρ 2 )(ρ 1 -ρ 3 ) + ρ n+2 2 (ρ 2 -ρ 0 (1 -q)) (ρ 2 -ρ 0 )(ρ 2 -ρ 1 )(ρ 2 -ρ 3 ) + ρ n+2 3 (ρ 3 -ρ 0 (1 -q)) (ρ 3 -ρ 0 )(ρ 3 -ρ 1 )(ρ 3 -ρ 2 )
),

for n ≥ 0. The expression of c n , for n ≥ 0, can be also easily derived in closed form in the other remaining cases for µ i and µ j for i ̸ = j ∈ {0, 1, 2, 3}. Finally, we deduce from Equation (3.7) that m n = pρ n 0 , for n ≥ 0. Because of the light-traffic approximation in the above computations, the system steady-state probabilities do not sum up to one. We then normalize them by dividing each one of them by the summation of all the steady-state probabilities. In Table 3.3 we evaluate the light-traffic approximation. The comparison of the approximate results with those from the exact numerical analysis show that the approximation is appropriate under the light-traffic regime. Note that the probabilities for n ≥ 4 are too close to zero for both the exact and the approximation methods.

Multi-Server Case

In this section, we focus on the email routing problem for the multi-server case. The modeling is the same as described in Section 3.3, expect that instead of one server, there are s identical, parallel servers. As previously, we consider a call center manager that wants to optimize the email routing as a function of the system parameters (Problem (3.1)). In other words, we want to either optimize p and q for Model PM, or give the ordering of the extreme cases Models 1 to 4 (easier to use in Table 3.3: Validation of the light-traffic approximation 0.0451% 0.0048% 0,1559% 0,0790% 0.0000% 0.0000% b 0.0607% 0.0145% 0,1114% 0,0691% 0,0001% 0.0000% b ′ 0.2374% 0.1378% 0,0648% 0,0432% 0.0000% 0.0000% c 0.4380% 0.3414% 0,1042% 0,0758% 0,0001% 0.0001% m 2 0.0606% 0.0755% 0,0459% 0,0493% 0,0012% 0.0012% a 3 0.0071% 0.0002% 0,0195% 0,0079% 0.0000% 0.0000% b 0.0091% 0.0008% 0,0146% 0,0074% 0.0000% 0.0000% b ′ 0.0382% 0.0153% 0,0087% 0,0048% 0.0000% 0.0000% c 0.0787% 0.0499% 0,0146% 0,0089% 0.0000% 0.0000% m 3 0.0055% 0.0076% 0,0011% 0,0012% 0.0000% 0.0000%

λ µ =28%, λ = 0.2, λ µ =20.83%, λ = 0.1, µ 1 = 1, λ µ =2.00%, λ = 0.01 µ 1 = µ 2 = 5, µ 3 = µ 0 = 2, µ 2 = 2, µ 3 = 3, µ 0 = 4, µ i = 2 for i = 0, • • • , 3 p = 10%, q = 80% p = q = 1 p = q = 50%
= 2, µ 1 = µ 3 = 1, µ 2 = 3 s = 10, µ 0 = µ 1 = µ 2 = µ 3 = 2 Constraint
-1.2 0 -1.1 0 -1.3 0 -1.2 P (W < 1) ≥ 0.8 Model 2 1.2 -1.9 1.1 -1.7 1.3 -2.0 1.2 -1.8 Model 1 1.9 -2.3 1.7 -2.1 2.0 -2.5 1.8 -2.4
practice and also efficient).

An exact analysis as that done for the single server case is too complex. We propose an approximation based on the single server results. It consists on replacing the s servers by a single super server. The service rates become sµ 0 (for emails), sµ 1 (first stage of call service), sµ 2 (second stage of call service), sµ 3 (third stage of call service). We have conducted an extensive simulation study in order to assess the quality of this approximation. Some of the comparison results between approximation and simulation are given in Table 3.4.

In Table 3.4 we give as a function of the interval of the call arrival rate value the ordering of Models 1 to 4 with respect to the optimization problem. The intervals are given using the single server approximation and also using a combined simulation and optimization of the multi-server system. The same intervals hold also when considering Model PM. We observe from Table 3.4 that the approximation provides an appropriate solution for the email routing optimization.

Conclusions

We considered a blended call center with calls and emails. The call service is characterized by successive stages where one of them is a break for the agent in the sense that inside the conversation there is no required interaction during a non-negligible time between the two parties. We addressed an important question in the call center practice: how should managers make use of this opportunity in order to better improve performance? We focused on the optimization of the email routing given that calls have a non-preemptive priority over emails. Our objective was to maximize the throughput of emails subject to a constraint on the call waiting time.

We developed a general framework (Model PM) with two probabilistic parameters for the email routing to agents. One parameter controls the routing between calls, and the other does the control inside a call conversation. We have also considered 4 particular cases corresponding to the extreme values of the probabilistic parameters. For these routing models, we have derived various structural results with regard to the optimization problem. We have also numerically illustrated and discussed the theoretical results in order to provide guidelines to call center managers. In particular, we proved for the optimal routing that all the time at least one of the two email routing parameters has an extreme value.

There are several avenues for future research. It would be interesting to extend the structural results to the multi-server case. It would also be useful but challenging to extend the analysis to cases with an additional channel, in particular the chat which is increasingly used in call centers.

Using the chat channel, an agent may handle many customers at the same time, which represent an additional opportunity to efficiently use the agent time.

Chapter 4

Adaptive Threshold Policies for

Multi-Channel Call Centers

In the context of multi-channel call centers with inbound and outbound calls, we consider a threshold policy on the reservation of agents for the inbound calls. We study a general non-stationary model where the call arrival follows a non-homogeneous Poisson process. The optimization problem consists on maximizing the throughput of outbound calls under a constraint on the waiting time of inbound calls. We propose an efficient adaptive threshold policy easy to implement in practice in the Automatic Call Distribution (ACD). This scheduling policy is evaluated through a comparison with the optimal performance measures found in the case of a constant stationary arrival rate, and also a comparison with other intuitive adaptive threshold policies in the general non-stationary case.

Introduction

Call centers require a very accurate match of demand and supply. The delay of the treatment of a call, its waiting time, is usually not allowed to exceed 20 seconds. Thus a very accurate prediction of the demand is required. However, this can rarely be obtained, because of the volatility of call arrival patterns. Therefore there is often a mismatch between demand and the scheduled supply, consisting of rostered call center employees (usually called agents). Moreover, even if the demand is accurately forecasted, a considerable overcapacity should be scheduled to be able to deal with the random Poisson fluctuations of the demand. Usually queueing models are used to quantify this overcapacity, most often Erlang C.

To prevent idle overcapacity and to limit the necessity to have extremely accurate forecast inbound calls are sometimes mixed with other types of customer contacts which have a less strict allowable delay, such as emails or outbound calls. This is called (call) blending. The amount of capacity assigned to the other channels is supposed to adapt to the number of inbound calls, giving at the same time a good service level for the inbound calls and a good occupancy of the call center agents.

Because of the strict waiting time requirement on inbound calls it is best to give them priority over the other channels. To maximize agent productivity it would be best to assign an email to every idle agent where there are no inbound calls in queue. This would lead to a 100% productivity.

There are two objections against this policy: a 100% occupancy is not sustainable over a longer period because of agent fatigue, and it would lead to long waiting times for inbound calls because never an agent is waiting for an inbound call to arrive. For these reasons a more sophisticated assignment policy is required.

In [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF] and [START_REF] Gans | A call-routing problem with service-level constraints[END_REF] it is shown that the optimal assignment policy is of the following form: outbound calls should only be scheduled when there are no waiting inbound calls and when the number of idle agents exceeds a certain threshold. Thus the problem of controlling our blended call center reduces to determining the right threshold level. This threshold however depends on all the parameters of the system. But these parameters, especially the arrival rate, are often hard to determine. This calls for a policy in which the threshold is adapted to the current situation without using explicitly the parameters of the system. In this chapter such adaptive policies are studied, both for systems with a constant (but unknown) arrival rate and for the more realistic situation of a fluctuating arrival rate. The parameter that is used to update the threshold is the service level up to that moment, a number which is always available in call centers.

The overall objective is to reach a certain service level by the end of the day, while maximizing the number of emails or outbound calls that are done.

We discuss the relevant literature. There is a rich literature on planning and scheduling in call centers, see Gans et al. (2003); Akşin et al. (2007). However, few papers focus on blending. The general context of multi-channel call centers is described in Koole ( 2013), Chapter 7. show in such a problem, which is close to ours, that a threshold policy outperforms a common routing policy used in call centers (that routes to the agent that has been idle the longest).

Milner and Olsen ( 2008) consider a call center with contract and non contract customers. They explore the use to give priority to contract customers only in off peaks. They show that this choice is a good one under classical assumptions (steady-states) and the same performance measure as ours and present also examples when it is not. This result is important since we found an insight arguing that the service level for inbound calls has to be very strictly respected during off peaks.

This chapter is organized as follows. Section 4.2 presents our model. Sections 4.3 and 4.4 contain our results, first for a constant arrival rate in Section 4.3 and then in Section 4.4 with a fluctuating arrival rate. We end with a short conclusion.

Model

We consider a call center modeled as a multi-server queuing system with two types of jobs, foreground jobs (inbound calls) and background jobs (outbound calls, emails, etc.). We use the terms foreground jobs and calls (background jobs and emails) interchangeably. The arrival process of calls is assumed to be a non-homogeneous Poisson process with rate λ(t), for t ≥ 0. Calls arrive at a dedicated first come, first served (FCFS) queue with infinite capacity. There is an infinite supply of background jobs, waiting for treatment in a dedicated FCFS queue. There are s identical, parallel servers (agents in call center parlance). Each agent can handle both types of jobs. We assume that the service times of foreground and background jobs are exponentially distributed with rates µ and µ 0 , respectively. Neither abandonment nor retrials are modeled.

Foreground jobs are more important than background ones in the sense that the former request a quasi-instantaneous answer (waiting time in the order of seconds or minutes), while the latter are more flexible and could be delayed for several (tens of) hours. The objective of the call center manager over a working day is to maximize the email throughput while satisfying a constraint on the call waiting time in the queue.

Since the model is transient, we can not define the waiting time of an arbitrary customer as a unique random variable. There is a random number of served customers during the working period S, if S > 0 we define the random variable for the waiting time of customer n for n ∈ {1, . . . , S}, W n . We want that the expected proportion of calls that wait less than a predefined threshold τ is at least equal to α, i.e., S -1 E S ∑ n=1

I{W n ≤ τ } ≥ α, for τ ≥ 0 and 0 ≤ α ≤ 1. Note that we do not consider arriving customers at the end of the working period which can not be served.

We then aim to find the best routing rules in terms of efficiency for the considered problem and easiness of implementation in call center software. We assume that preemption of jobs in service is not allowed. This is a quite natural assumption. An agent usually prefers to finish answering an underway email rather than starting it over later on. This is also preferred from an efficiency perspective. Evidently, when the background jobs are outbound calls, then it is neither acceptable to preempt.

For a similar model but with a constant arrival rate and equal service requirements for the two job types, [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF] prove that the optimal policy is a threshold policy with the priority given to calls (some servers reserved for calls). Their result is mainly based on the fact that it is optimal to handle calls as long as the queue of calls is not empty. For our general modeling, the analysis is more complicated. Even for a constant arrival rate but different service requirements, the optimal policy is hard to obtain, and might not be useful in practice (for software implementation for example). For simplicity and usefulness of the results in practice, we then restrict ourselves to the case of threshold policies. Moreover, [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF] numerically show, for more general cases, that the appealing threshold policies are good approximations of the optimal ones.

More concretely, the functioning of the call center under a threshold policy is as follows. Let us denote the threshold by u, 0 ≤ u ≤ s. Upon arrival, a call is immediately handled by an available agent, if any. If not, the call waits in the queue. When an agent becomes idle, she handles the call at the head of the queue with calls, if any. If not, the agent may either handle an email at the head of the email queue, or she remains idle. If the number of idle agents (excluding her) is at least s -u, then the agent in question handles an email. Otherwise, she remains idle. In other words, there are s -u agents that are reserved for calls, u is the number of agents that are always working.

In this chapter, we propose an adaptive threshold policy which adjusts the threshold as a function of the non-stationary arrival process of calls. We divide the working day into N identical intervals, each with length θ. The total working duration in a day is L, L = N θ. At the beginning of each interval i (i = 1, ..., N ), we define the threshold u i , 0 ≤ u i ≤ s, under which the job routing policy works during interval i. Let T denote the expected throughput of emails over the whole day, i.e., the ratio between the number of treated emails and L. Let also SL be the proportion, for the whole day, of calls that have waited less than τ ,

SL = S -1 E S ∑ n=1
I{W n ≤ τ }. In summary, our optimization problem can be formulated as

   Maximize T subject to SL ≥ α, (4.1)
where the decision variables are u i with 0 ≤ u i ≤ s, for i = 1, ..., N . It is clear that the best case for calls is such that u i = 0 for all i, which means that no email is treated and SL is maximized (case of an M(t)/M/s with only calls). We therefore assume from now on that the parameters λ(t)

for t ≥ 0, µ and s are such that SL ≥ α for u i = 0 (i = 1, ..., N ).

Constant Arrival Rate

We consider a basic case with a constant arrival rate, λ(t) = λ for t ≥ 0 and a constant threshold, u i = u for i = 1, ..., N and 0 ≤ u ≤ s. The purpose of the analysis in this section is to understand the behavior of the performance measures as a function of the threshold in order to build an efficient method for the threshold adaptation rule (u i for i = 1, ..., N ) in the case of a non-constant arrival rate. In Section 4.3.1 we propose a method to compute the performance measures, then in Section 4.3.2 we use them to provide a useful insight to construct our adaptive policy.

Performance Measures

In Section 4.3.1 we provide close form formula of the performance measures in the case of equal service rates and study the form of these measures as a function of the threshold. Then in Section 4.3.1 we propose a numerical method to compute the performance measures in the case of unequal service rates. Since we consider a stationary model we can define a unique random variable for the waiting time of an arbitrary customer W , and denote by P (W < τ ) the probability that an arbitrary customer waits less than τ (τ > 0).

Equal Service Rates

We consider the case µ = µ 0 . First, we compute the performance measures of interest for calls and emails for a given constant reservation threshold, denoted by u, 0 ≤ u ≤ s. We then develop some structural results that will be used in Section 4.3.2.

Let us define the stochastic process {x(t), t ≥ 0}, where x(t) ∈ N is the number of jobs in service plus the number of jobs in the queue of calls. Since µ = µ 0 , we need not distinguish between the two job types in service. The process {x(t), t ≥ 0} is a birth-death process. The transition rate from state x to state x -1 is min{x, s}µ, for x ≥ 1, and that from state x to state x + 1 is λ, for x ≥ 0. Under the stability condition λ sµ , we denote by p x the steady-state probability to be in state x ∈ N, and by a the ratio λ µ . In Theorem 2, we give the expression of the email throughput, T (s, u, a), and that of the probability that the call waiting time is less than τ , SL = P (W < τ ).

Theorem 2 For 0 ≤ u ≤ s, we have

T (s, u, a) = µ ( s-u ∑ k=0 a k u! (u + k)! + a s-u u! s! a s -a ) -1 ( u + s-u ∑ k=1 a k u! (u + k -1)! + a s-u+1 u! (s -1)!(s -a) ) -λ, (4.2) P (W < τ ) = 1 -C(s, u, a)e -τ (sµ-λ) , (4.3) with C(s, u, a) = a s-u u! s!(1 -a/s) ( s-u ∑ k=0 a k u! (u + k)! + a s-u u! s! a s -a ) -1 . (4.4) Proof. For 0 ≤ x < u, we have p x = 0. For 0 ≤ k ≤ s -u, we have p u+k = a k u! (u+k)! p u . For k ≥ 0,
we have p s+k = a k s k p s . Since all probabilities sum up to one, we obtain

p u = ( s-u ∑ k=0 a k u! (u + k)! + a s-u u! s! a s -a ) -1 . (4.5)
The email throughput can be seen as the overall throughput (of calls and emails) minus the call throughput. Thus

T (s, u, a) = s-u ∑ k=0 (u + k)µp u+k + sµ ∞ ∑ k=1 p s+k -λ.
After some algebra, we deduce that

T (s, u, a) = µp u ( u + s-u ∑ k=1 a k u! (u + k -1)! + a s-u+1 u! (s -1)!(s -a) ) -λ.
Note that the lower bound of T (s, u, a) is T (s, 0, a) = 0, which corresponds to the case when all servers are reserved to calls. As for the upper bound, it is T (s, s, a) = sµ -λ, which corresponds to the case of no server reservation for calls (the infinite amount of emails leads to sµ as a total throughput for the two job types).

The call service level, P (W > τ ), is obtained using the PASTA property. We have

P (W > τ ) = ∞ ∑ n=0 p s+n P (W > τ |x = n + s)
, where P (W > τ |x = s + n) is the conditional probability that the waiting time of a new call exceeds τ , given that it finds all servers busy and n calls waiting ahead in the queue, n ≥ 0. It is easy to see that this conditional waiting time follows an Erlang distribution with n + 1 stages and a rate of sµ per stage. Then,

P (W > τ |x = s + n) = n ∑ k=0 e -sµτ (sµτ ) k k! , which leads to P (W > τ ) = ∞ ∑ n=0 p s a n s n n ∑ k=0 e -sµτ (sµτ ) k k! = lim n→∞ ( p s e -sµτ n ∑ k=0 ∞ ∑ n=k (sµτ ) k k! ( a s ) n ) .
Observing that

∞ ∑ n=k ( a s ) n = ( a s ) k 1 1-a/s implies P (W > τ ) = C(s, u, a)e -τ (sµ-λ) , (4.6) with C(s, u, a) = p s 1 -a/s = a s-u u! s!(1 -a/s) ( s-u ∑ k=0 a k u! (u + k)! + a s-u u! s! a s -a ) -1 . (4.7)
Note that the upper bound of C(s, u, a) is C(s, s, a) = 1 (no server reservation for calls, then, any arriving call has to wait for service), and its lower bound is C(s, 0, a) =

a s s!(1-a/s) ( s ∑ k=0 a k (k)! + a s s!(1-a/s) ) -1
(all servers are reserved to calls, which corresponds for calls to the case of an M/M/s queue with no emails). This finishes the proof of the theorem. ✷

In Section C.1 of the appendix, we prove some monotonicity results of the performance measures in the threshold. We prove that the email throughput is strictly increasing and neither convex nor concave in u, for 0 ≤ u ≤ s. However the end of the email throughput, for s -2 ≤ u ≤ s and s ≥ 2, is concave in u. The inbound service level P (W < τ ) is strictly decreasing in u, for 0 ≤ u ≤ s.

We prove that it is concave in u, for the cases a < 1 and also for a ≥ u + 1 (u < s). An extensive numerical study for the remaining cases shows that the concavity still holds.

Unequal Service Rates

In this section we focus on the performance evaluation (email throughput and call waiting time distribution) for the case of unequal service rates, µ ̸ = µ 0 . In contrast to the case of equal service rates, the performance expressions are here too cumbersome to allow the development of useful structural results. The results of this section are however still useful for the numerical experiments in Section 4.3.2 in order to build the insights on the threshold policy for the more general case with a non-constant call arrival rate.

As in [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF], our approach consists on using a Markov chain analysis to derive the steady-state probabilities of the system, from which the performance measures are characterized thereafter. To simplify the presentation, we focus on the particular case u = s. The analysis for the case u = 0 is obvious, and that of the remaining cases, 0 < u < s, is done similarly to the case u = s. It simply adds a finite number of additional equations but does not impact the general form of the steady-state probabilities. Consider the stochastic process {(x(t), y(t)), t ≥ 0}, where

x(t) is the number of waiting calls in the queue and y(t) is the number of emails being in service, (x, y) ∈ N 2 . This process is a Markov chain. For x ≥ 0 and y ≥ 0, the transition rate from (x, y)

to (x + 1, y) is λ. For x ≥ 1 and y ≥ 0, the transition rate from (x, y) to (x -1, y) is (s -y)µ. For

x ≥ 1 and y ≥ 1 the transition rate from (x, y) to (x -1, y -1) is yµ 0 . For y ≥ 0, the transition rate from (0, y) to (0, y + 1) is (s -y)µ. Under the stability condition λ sµ < 1, we denote by p x,y the steady-state probability that the system is in state (x, y).

For y = s and x > 0, we have p x,s (λ

+ sµ 0 ) = λp x-1,s , then p x,s = ( λ λ+sµ 0 ) x p 0,s . Using λp 0,s = µp 0,s-1 , we deduce that p x,s = µ λ ( λ λ+sµ 0 ) x p 0,
( λ + µ + (s -1)µ 0 + √ (λ + µ + (s -1)µ 0 ) 2 -4λµ
)

and z 2 = 1 2µ ( λ + µ + (s -1)µ 0 - √ (λ + µ + (s -1)µ 0 ) 2 -4λµ
) . We deduce, for x ≥ 0, that

p x,s-1 = β 1 (z 1 ) x + β 2 (z 2 ) x + β 3 ( λ λ + sµ 0 ) x ,
with β 3 = -µs λ+(µ 0 -µ)s p 0,s-1 . From the boundaries x = 0 and x = 1, we obtain p 0,s-

1 = β 1 + β 2 - µs λ+(µ 0 -µ)s p 0,s-1 and p 1,s-1 = β 1 z 1 + β 2 z 2 - µs λ+(µ 0 -µ)s ( λ λ+sµ 0
) p 0,s-1 , respectively, which implies

β 1 = 2µp 0,s-2 (λ + s(µ 0 -µ)) + p 0,s-1 (z 2 µµ 0 s -λµ 0 s + λµs + z 2 λµ -λµ -λ 2 (z 2 -z 1 ) µ ((µ 0 -µ)s + λ) , β 2 = 2µp 0,s-2 (λ + s(µ 0 -µ)) + p 0,s-1 (z 1 µµ 0 s -λµ 0 s + λµs + z 1 λµ -λµ -λ 2 (z 2 -z 1 ) µ ((µ 0 -µ)s + λ) .
For y = k, 0 ≤ k < s, and x > 0 we have

(λ + (s -k)µ + kµ 0 )p x,k = λp x-1,k + (s -k)µp x+1,k + (k + 1)µ 0 p x+1,k+1 . (4.8)
The homogeneous equation associated to Equation (4.8) is

(s -k)µz 2 -(λ + (s -k)µ + kµ 0 )z + λ = 0,
with z as a variable, for z ∈ C. It has two solutions denoted by z 1,k and z 2,k and are given by

z 1,k = 1 2(s -k)µ ( λ + (s -k)µ + kµ 0 - √ (λ + (s -k)µ + kµ 0 ) 2 -4(s -k)λµ ) , z 2,k = 1 2(s -k)µ ( λ + (s -k)µ + kµ 0 + √ (λ + (s -k)µ + kµ 0 ) 2 -4(s -k)λµ ) ,
for 0 ≤ k < s. Because of the last term in the right hand side of Equation (4.8), one may write, for 0 ≤ k ≤ s and x > 0,

p x,k = s ∑ i=k A i,k z x 1,i + B i,k z x 2,i ,
with z 1,s and z 2,s defined as z 1,s = λ λ+sµ 0 and z 2,s = 0, respectively, and

A i,k , B i,k ∈ R for 0 ≤ k < s
and k ≤ i ≤ s. Using Equation (4.8), we can prove that

A i,k+1 = A i,k -(s -k)µz 2 1,i + (λ + (s -k)µ + kµ 0 )z 1,i -λ (k + 1)µ 0 z 2 1,i
, for 0 ≤ k < s and k < i ≤ s. Similarly, we have

B i,k+1 = B i,k -(s -k)µz 2 2,i + (λ + (s -k)µ + kµ 0 )z 2,i -λ (k + 1)µ 0 z 2 2,i , for 0 ≤ k < s and k < i ≤ s. For i = k,
(λ + (s -k)µ)p 0,k = (s -k)µp 1,k + (k + 1)µ 0 p 1,k+1 + (s -k -1)p 0,k-1 ,
for 0 ≤ k < s. We also have z 1,0 = 1 and z 2,0 = λ sµ . Since all probabilities sum up to one, we obtain A 0,0 = 0. In conclusion through the above analysis, we have characterized all steady-state probabilities, p i,k , for i ≥ 0 and 0 ≤ k ≤ s. The email throughput T (λ, µ, µ 0 , s) may be written as

T (λ, µ, µ 0 , s) = µ 0 s ∑ k=1 ∞ ∑ i=0 kp i,k , or equivalently T (λ, µ, µ 0 , s) = λ + s ∑ k=1 kµ 0 p 0,k ,
for s ≥ 1. As for the call waiting performance, it is given by

P (W > τ ) = s ∑ k=1 ∞ ∑ i=0 p i,k P (W > τ |(x, y) = (i, k)),
where P (W > τ |(x, y) = (i, k)) is the conditional probability that the waiting time of a new call exceeds τ , given that it finds i emails in service, s -i calls in service, and k calls waiting ahead in the queue, for 0 ≤ i ≤ s and k ≥ 0. The computation of P (W > τ |(x, y) = (i, k)), for 0 ≤ i ≤ s and k ≥ 0, is as follows. For k = 0 and 0 ≤ i ≤ s, the new call has to wait for a service completion of one of the i emails, or one of the s-i calls, so, P (W > τ |(x, y) = (i, 0)) = e -τ (iµ 0 +(s-i)µ) . For k = 1 and 0 ≤ i ≤ s, the probability that the next service completion is that of an email is iµ 0 iµ 0 +(s-i)µ . Thus, the waiting time of the new call follows a hypoexponential distribution consisting of the summation of two exponential random variables with rates iµ 0 + (s -i)µ and max(0, i -1)µ 0 + min(s, s -i + 1)µ with probability

iµ 0 iµ 0 +(s-i)µ
, and it follows an Erlang distribution with 2 phases and iµ 0 + (s -i)µ

as a rate per stage with probability 1 -iµ 0 iµ 0 +(s-i)µ . This leads to

P (W > τ |(x, y) = (i, 1)) = iµ 0 iµ 0 + (s -i)µ × ((i -1)µ 0 + (s -i)µ)e -τ (iµ 0 +(s-i)µ) -(iµ 0 + (s -i)µ)e -τ ((i-1)µ 0 +(s-i)µ) µ -µ 0 + (s -i)µ iµ 0 + (s -i)µ e -τ (iµ 0 +(s-i+1)µ) (1 + τ (iµ 0 + (s -i)µ)),
for 0 ≤ i ≤ s. One can continue in the same way to derive all the conditional waiting time probabilities for k > 1, which finishes the characterization of the performance measures (email throughput and call waiting time distribution) in the case of unequal service rates.

Construction of the Adaptative Threshold Policy

In this section, we use the previous results to find an insight on how we should adapt the threshold as a function of the intensity of the call arrivals. The objective is to maximize the throughput of emails while reaching the constraint on the call waiting times for the whole day. We find that during the periods with low demand, the need of having a good service level is more important than during the periods with high demand. On the basis of this observation, we build a method for adapting the threshold. We then evaluate this method by comparing it with the optimal threshold policy.

Numerical Observations

For a given time interval long enough to reach the stationary regime, one can use the results of Section 4.3.1 to obtain the optimal threshold, denoted by u * , for Problem (4.1). Consider now a working day with two time intervals, each with a different call arrival rate, and on each of which the stationary regime is reached. We want to find the optimal couple of thresholds that answers our optimization problem, where the call service level constraint is for the whole day. We denote the first (second) time interval by I 1 (I 2 ) and its arrival rate by λ 1 (λ 2 ). Without loss of generality, we consider cases where λ 1 ≤ λ 2 .

In Table 4.1, we consider various scenarios of arrival rates, service rates, and relative durations between the two time intervals. We give the optimal threshold of each interval in isolation (i.e. the highest threshold which verifies the service level constraint). They are denoted by u * 1 and u * 2 for I 1 and I 2 , respectively. We also evaluate the couple of thresholds which answers Problem (4.1) on the set of the two intervals. This couple is found by an exhaustive test of all the possible values for the couple (u 1 ,u 2 ). We denote by (u 1 , u 2 ) * this optimal couple. Remark that for this couple, Problem (4.1) does not have to be answered on each interval but on the set of the two intervals. Finally, we give the performance measures for each interval and for the set of the two intervals for the couple

(u 1 , u 2 ) * .
We observe that u * 1 (respectively u * 2 ) is always higher or equal to u 1 (respectively lower or equal to u 2 ) in the optimal couple (u 1 , u 2 ) * . An interesting insight here is that, while respecting the global call service level, we should strictly respect the service level during the interval with a small arrival rate (I 1 ), and violate the constraint when the arrival rate is high (I 2 ). The reason is related to the sensitivity of the service level to an increase of the threshold. When the workload increases the sensitivity of the service level for a given threshold (∆SL(u) = SL(u+1)-SL(u) for 0 ≤ u < s)

first increases and then decreases (except for u = s -1, the sensitivity only decreases). This can be seen in Figure 4.1. Proposition 4 proves that there is less waste for the service level, when increasing the threshold in situations for which the sensitivity of the service level is decreasing in the workload.

Proposition 4 If λ 1 < λ 2 and ∆SL(u) is decreasing in the workload then

|∆SL λ 1 (u * 1 )| ≥ |∆SL λ 2 (u * 2 )|. Proof. If u * 1 = u * 2 then increasing u is less sensitive in SL l 2 than in SL l 1 since the sensitivity of SL is decreasing in the workload. The other case is u * 1 > u * 2 because a 1 < a 2 . Since SL is decreasing and concave in u, SL l 1 is more sensitive to the increasing of u from u * 1 than from u * 2 . From u * 2 ,
SL l 1 is more sensitive to the increasing of u than SL l 2 , then SL l 2 is less sensitive to the increasing of u from u * 2 than SL l 1 would be from u * 1 . ✷

Yet the opposite seems to be more intuitive, since it would be hard to compensate a very bad service level during an interval with a high number of calls. Counterexamples can be found when

λ 1 << λ 2 . For example for λ 1 = 0.1, λ 2 = 1, µ 0 = µ = 0.2, s = 10, I 1 = I 2 = 50%, α = 80%, we
have u * 1 = 9, u * 2 = 8 and (u 1 , u 2 ) * = (10, 8).

Our Adaptive Threshold Policy (ATP)

We propose for Problem (4.1) an adaptive threshold policy which adjusts the threshold as a function of the call workload. The idea of the policy comes from the numerical observations in Section 4.3.2. As mentioned in Section 4.2, the threshold is reevaluated at the beginning of each interval i (i = 1, ..., N ). The threshold associated to interval i is denoted by u i . The global service level for the whole day (all N intervals) is denoted by SL, and the global one from interval 1 to interval i is denoted by SL i , for i = 1, ..., N .

If SL i is higher (lower) than α at the beginning of an interval i (i = 2, ..., N ) then the policy increases (decreases) the threshold. To change the threshold, we use a real parameter denoted by c i (i = 1, ..., N ). The threshold u i is defined as the closest integer to c i , for i = 1, ..., N . Note that the parameter c i is chosen to be real in order to smooth the change in the threshold u i . We start with u 1 = c 1 = s. For i ≥ 2, if we need to increase the threshold, then we take

c i = c i-1 + 1 -c i-1 /s. If not, then c i = c i-1 -c i-1 /s.
This policy is refereed to as ATP.

As the workload of calls decreases, ATP increases the threshold with a decreasing speed. This decreasing speed allows a slow increase in the threshold and then gives advantage to calls, which is coherent with the insight of Section 4.3.2. The opposite is also true and coherent with the insight.

The advantages of ATP is that it is simple, easy to implement, and at the same time efficient as we will show later.

Evaluation of the Adaptative Threshold Policy

In this section, we evaluate the quality of the ATP policy by comparing it with the optimal one.

First we provide the optimal threshold policy. Because of the discrete nature of the threshold, one may intuitively see that the threshold should vary between two or more values. The reason is that we need to satisfy exactly the constraint on calls in Problem (4.1) in order to maximize the email throughput. Both for cases µ 0 = µ and µ 0 ̸ = µ, Theorem 3 provides a weak condition under which the optimal policy is a randomization of the threshold between two values.

Theorem 3 Consider 0 ≤ u 1 , u 2 ≤ s such that SL(u 1 ) ≤ α ≤ SL(u 2 ). If there exists γ ∈ R for which (u 1 , u 2 ) ∈ arg max u T (u) + γSL(u), then randomizing between u 1 and u 2 is optimal.

Proof. Let p ∈ [0, 1] be the parameter of randomization between u 1 and u 2 . Assume that we can find a couple (u 3 , u 4 ) ̸ = (u 1 , u 2 ) and a parameter of randomization q ∈ [0, 1] such that the constraint on calls is also saturated and SL(u 3 ) ≤ α ≤ SL(u 4 ). We have pT

(u 1 ) + (1 -p)T (u 2 ) + γpSL(u 1 )+γ(1-p)SL(u 2 ) ≥ qT (u 3 )+(1-q)T (u 4 )+γqSL(u 3 )+γ(1-q)SL(u 4 ). Since γpSL(u 1 )+ γ(1 -p)SL(u 2 ) = γqSL(u 3 ) + γ(1 -q)SL(u 4 ) = γα, we deduce that pT (u 1 ) + (1 -p)T (u 2 ) ≥ qT (u 3 ) + (1 -q)T (u 4 ).
Then the couple (u 1 , u 2 ) is optimal, which completes the proof. ✷

The randomization between two thresholds allows for the constraint on calls to be met exactly.

For our system with constant parameters, we believe that the randomization is between two suc-Table 4.2: Comparison under steady-states assumption (θ=15min) Optimal c Optimal T ATP T Difference Scenario 1 (λ = 4, µ = µ 0 = 0.2, s = 28) 25.49 1.39 1.37 1.46% Scenario 2 (λ = 0.02, µ = µ 0 = 0.2, s = 1) 0.13 0.02 0.02 0.00% Scenario 3 (λ = 18, µ = µ 0 = 0.2, s = 100) 93.91 1.65 1.58 4.43% Scenario 4 (λ = 4, µ = 0.27, µ 0 = 0.15, s = 28) 26.63 1.89 1.89 0.00% Scenario 5 (λ = 4, µ = 0.17, µ 0 = 1, s = 28) 23.21 2.00 1.79 11.73% cessive thresholds. Since the throughput is neither convex nor concave it is difficult to rigorously prove this result. However, if we denote by u * (0 ≤ u * ≤ s) the highest threshold that verifies SL(u * ) > α, we numerically checked that with γ = -T (u * +1)-T (u * ) SL(u * +1)-SL(u * ) (for 0 ≤ u * < s), the expression T (u)+γ ×SL(u) is strictly increasing from u = 0 to u = u * , strictly decreasing from u = u * +1 to u = s and T (u * ) + γSL(u * ) = T (u * + 1) + γSL(u * + 1). Then for all the considered numerical situations the optimal policy is a randomization between two adjacent values when 0 ≤ u * < s.

When u * = s, the optimal policy is to keep the threshold constant and equal to s.

In Table 4.2, we propose 5 representative scenarios with constant arrival rates and compare the optimal throughput with the one found with under ATP. A comprehensive numerical study can be found in Section C.2 of the appendix. Although the ATP method is not optimal, the difference with the optimum is quite small. This shows the advantage of ATP in the case of constant arrival rates. Recall that our main purpose in this chapter is the analysis of the case with a fluctuating arrival rate. In the next section, we consider the case of a fluctuating arrival rate and evaluate the performance of ATP through a comparison with other intuitive methods.

Non-Constant Arrival Rates

In Section 4.4.1 we compare ATP with methods that use constant step sizes. Then in Section 4.4.2

we analyze the impact of the parameters on the choice of the method. In Section 4.4.3 we propose some other intuitive adaptive methods.

We consider cases where the length of the working day equals eight hours (L = 8h) and a frequent possibility of reevaluating the real threshold c, at the beginning of each time interval with length θ = 1 min, 5 min or 15 min. We use simulation to obtain the performance measures.

For each scenario, we run n replications. We then introduce a measure of the bias after the n simulations, denoted by r n and calculated as

r n = n ∑ k=1 Max(α-SL k ,0) n
, where SL k is the expected service level of simulation k (1 ≤ k ≤ n). Since the value of r n should be as small as possible, we introduce a coefficient A which would be the aversion of the call center manager to the risk and introduce an utility indicator denoted by U n and given by T n -A × r n , where T n is the average throughput after n simulations. The confidence intervals are a safe way to evaluate the required number of equivalent simulations, n. The confidence interval for a proportion p and a risk of 5% is

( p -1.96 √ p(1-p) n , p + 1.96 √ p(1-p) n
) in which n is the number of terms used to calculate the proportion p. If we want a precision of one decimal we need 2 × 1.96

√ 0.8(1-0.8) n < 0.001 then n >2
458 624. In order to have safe results we run each simulation 3 000 000 times.

Comparison with Constant Step Methods

We propose different scenarios to compare ATP with constant step size methods. We denote by h the step size (0 < h ≤ 1). When we need to increase (respectively decrease) the real threshold c i after i intervals (1 ≤ i < N ) under the case SL i > α (respectively SL i < α) we add h to c i (respectively we add -h to c i ). In each scenario we use an aversion of risk equal to 100 and initialize the system with c 0 = u 0 = s. In some scenarios the number of agents varies over the day.

When the number of agents decreases, we could be in a situation in which c > s, i.e., the number of busy agents becomes higher than the new value for s. To avoid such a situation, we force in the simulation the change of c to the new smaller value of s. Any undertaken task by a removed agent is lost. In all scenarios the constraint on calls is such that the proportion of calls that wait less than 30 seconds is at least 80%, τ = 30s and α =80%. We consider the following scenarios:

• Scenario 1: λ = 4, µ = µ 0 = 0.2, s = 28 and N = 480 (θ = 1 min),

• Scenario 2: λ = 4, µ = µ 0 = 0.2, s = 28 and N = 32 (θ = 15 min),

• Scenario 3: λ = 4, µ = 0.27, µ 0 = 0.15, s = 28 and N = 480,

• Scenario 4: λ = 4, µ = 0.17, µ 0 = 1, s = 28 and N = 480,

• Scenario 5: λ linearly decreasing from 5 to 3 , µ = µ 0 = 0.2, s = 34 if λ > 4.5, s = 28 if 4.5 > λ > 3.5, s = 23 in the remaining cases, and N = 480,

• Scenario 6: λ linearly increasing from 3 to 5 , µ = µ 0 = 0.2, s = 34 if λ > 4.5, s = 28 if 4.5 > λ > 3.5, s = 23 in the remaining cases, and N = 480,

• Scenario 7: During the first quarter of the period λ is linearly increasing from 1 to 5, during the second quarter λ is linearly decreasing from 5 to 3, during the third quarter λ is linearly increasing from 3 to 5 and during the last quarter λ is linearly decreasing from 5 to 1, µ = µ 0 = 0.2, s = 34 if λ > 4.5, s = 28 if 4.5 > λ > 3.5, s = 23 in the remaining cases, and

N = 480,
• Scenario 8: The period T is divided into 10 sub-periods and the value of λ alternates between the values 5 and 0.5, i.e., it is 5 in the first sub-period, 0.5 in the second one, again 5 in the third one, and so on, µ = µ 0 = 0.2, s = 28 and N = 480.

The results are presented in Table 4.3. We consider values of h equal to 0.1, 0.2, 0.5 and 1.

We observe that ATP performs better or at least similarly to the constant step methods with an aversion of risk equal to 100.

In Figures 4.2(a), 4.2(b) and 4.2(c), we present the evolution of the threshold, the proportion of customers that wait less than 30 seconds and the email throughput as a function of time in one simulation of scenario 2. This is an illustration that could help to understand why ATP is efficient.

With a small value of h (h = 0.2), the initialization has an important impact on the evolution of the threshold. At the beginning with u 0 = c 0 = s = 28, there is a need to decrease the threshold. A small value of h does not allow to do this decreasing quickly enough. Then there is a need to keep on decreasing the threshold in order to have a chance to reach the service level on calls over the whole day. On the other hand a high value of h (h = 1) goes with a fluctuation of the threshold, with sometimes bad call service levels and other times bad email throughput. Note that the higher is h, the faster the service level converges its target. In what follows we go further in analyzing the impact of the main parameters on the choice of h.

Impact of the Parameters

In this section, we analyze the impact of the parameters in the choice of a constant value for h. This might permit opportunities for providing new methods for adapting h and also help to understand the performance of ATP.

Impact of the Number of Intervals, N : The comparison between scenarios 1 and 2 in Table 4.3 indicates that there is a link between h and N . In scenario 2 with only 32 intervals, a small 1,34 1,33 1,31 1,30 SL 80,03% 80,01% 80,03% 80,04% value of h does not allow to reach the call service level constraint. In scenario 1, a large number of intervals and a high value for h lead to an important fluctuation from u = 0 to u = s. For a given value of h there should be an optimal number of intervals. In Table 4.4 we present an example of optimal values of N for a given h. We observe that increasing N decreases in h.

An advantage of ATP is its ability of adaptation to the number of intervals. A high number of intervals could lead to a high probability of reaching extreme and inefficient states (u = 0 and u = s usually). Thanks to the slowing in the speed of the threshold reduction when c is small and the slowing in the threshold increasing when c is high, this is unlikely to happen. ATP also provides a high capacity of reaction when the threshold is too low or too high, which is important in the case of a small number of intervals.

Impact of the Aversion of Risk A:

The choice of the method for adapting the threshold depends on the risk aversion of the manager. In our simulations we considered that A = 100 and showed the efficiency of ATP. This value provides a balance between the performance (T ) and the risk (r). One further question would be about the choice of the manager if her aversion to the risk is different. If we consider the extreme case of an infinite aversion to the risk (A = ∞), the choice will be made for the smallest value of r. We observe that the choice will still be for ATP rather than the other methods even more than with A = 100. This is an important advantage of ATP; it is a safe method (i.e., the probability the reach the service level constraint at the end of the working period is high). Since the threshold c is usually closer to s than to 0 because of the concavity in the service level we usually have a higher speed in decreasing the threshold than in increasing it, which is safe and explains the small values for r. On the other hand if the manager has no risk aversion (A = 0) then the choice will be made for the highest average throughput (T ).

ATP is then not the best one but it still provides results close to the best ones in Table 4.3.

Impact of the Arrival Process: From Table 4.3 we observe that the choice of a constant step h is influenced by the volatility of the arrival process (comparison between scenarios 1 and 8). When the arrival rate is constant the optimal solution may be an oscillation between a threshold u and a threshold u+1. Then, in order to reach this solution, we need small values of h in order not to leave too often those two states. When the arrival rate has a high volatility, there is a need to change the threshold very often in order to reach a good value of u in a short time. This explains why h = 0.2 is better than h = 0.5 in scenario 1 and the opposite is true in scenario 8. The comparison between scenarios 5 and 6 is also interesting; in scenario 5 the arrival rate decreases over the day and in scenario 6, it increases. We observe that small values of h are better in scenario 6 because of the initial conditions; u 0 = c 0 = s. In scenario 5 those initial conditions with a small value of h lead to a difficulty in reaching a good threshold. We observe that ATP performs well in various scenarios of volatility in the arrival process thanks to its capacity of creating slow or fast changes in the threshold.

Impact of the Emails Service Rate: Consider scenarios 3 and 4. We observe that ATP performs better when the emails are served slower (scenario 3) than when they are served faster (scenario 4) than the calls. Because of the concavity of the service level, the threshold is usually closer to s than to 0. This implies a higher speed in decreasing than in increasing the threshold.

When the emails are served faster than the calls, the need to increase c is more important because an email does not occupy an agent for a long period of time but with our method this increasing might be too slow. However, we notice that this case has less meaning for our study since the problem of reserving agents is interesting in the case of long service times for background jobs (relatively to calls).

Comparison with other Intuitive Methods

In this section, we compare ATP with other intuitive adaptive methods. We propose the following ones based on the reevaluation of the step h i after each intervals i (i = 1, ..., N ).

Method 1: The first intuitive idea is to propose a decision based on the distance from the achieved service level and the target after each interval. The intuition is that the need to change the threshold increases with this distance. We initialize with h 0 = 0, c 0 = u 0 = s and SL 0 = 100%.

After each interval we reevaluate h i according to the relation:

h i+1 = Min {1, |SL i -α|} , for i ∈ {0, ..., N -1}.
Method 2: Method 2 is a variation of Method 1. We propose a decision based on the cumulative distance with the service level target, α. The intuition is that the need to change the threshold not only increases with the distance to the target service level, but also increases with the time spent above or under this target. More precisely, we initialize by h 0 = 0, c 0 = u 0 = s and SL -1 = SL 0 = 100%. After each interval we reevaluate h i according to the relation:

h i+1 = Min {1, h i + |SL i -α|} × I (SL i -α)(SL i-1 -α)>0 , for i ∈ {0, ..., N -1}.
Method 3: We propose the same evaluation of h i as in Method 2 but instead of using the service level SL i of the last i intervals (i = 1, • • • , N ), we use the service level measured only on the last interval i (i = 1, • • • , N ). This method is made to correct a too important weight that could be given to the past in the previous method.

Methods 4a and 4b: Methods 4a and 4b are not really intuitive. The idea behind them is the question of when the strongest decisions in the change of the threshold should be taken. If we choose the strongest changes in the threshold at the beginning of the period we could quickly reach the service level constraint (Method 4a). If we choose the strongest changes at the end of the period we could maximize the email throughput at the beginning and do an efficient correction at the end of the working period to reach the service level constraint (Method 4b). More precisely, in Method 4a we propose after i intervals to choose

h i = 1 - i N ,
and in Method 4b we choose

h i = i N ,
for i = 1, • • • , N .
We compare the proposed methods in Table 4.5 with the constant step sizes methods and ATP under scenarios 1 and 2. We observe that those methods are not as good as ATP and even sometimes not as good as the constant step size methods. Methods 4a and 4b are not efficient for a simple reason; the choices in changing the threshold mainly depend on the demand and not on the closeness to the end of the working day. We observe on other simulations that Method 4a is efficient when the variability in the demand is high at the beginning of the working period and the opposite is true for Method 4b. Although Methods 1 and 2 are the most intuitive, we observe that they are not efficient. The weight of the past is too heavy and entails extreme choices in the threshold (which are often inefficient) so as to compensate the past values. Method 3 is often more efficient in terms of the email throughput, however it converges very slowly. We observe on other simulations that Method 3 could be a good proposition only if a working day is long enough (at least 1000 hours). An intermediate solution between Methods 2 and 3 would be to propose a decision in the changes of the threshold based on the average of the service levels measured on all past intervals weighted by coefficients which are increasing with the proximity to the last interval.

Many solutions can be proposed in that direction but none of them seems to be efficient for a representative number of scenarios.

Conclusion and Future Research

We considered call centers with inbound calls and an infinite supply of non-preemtable outbound jobs. We proposed a scheduling policy, refereed to as ATP, where the objective is to do as much outbound as possible while satisfying a service level constraint on the call waiting time. In the reallife call center context with fluctuating call arrival rate, the assignment policy for outbound adapts itself to the current service level. We showed the efficiency of ATP by comparing its performance with those of other policies. One of the main advantage of ATP is its ability to quickly react when an important change in the arrival process happens and also its ability to avoid inefficient states when the arrival rate remains constant.

Future research on this subject may follow two directions. First, a theoretical modeling for the adaptive blending might be useful to better understand ATP. This is now hindered by the fact that no theory seems to exist on this type of control problems. One of the difficulties to build a Markov chain is the non-exponentiality of the decision interval length defined in the ACD. Another difficulty is the lack of transient results for the performance measures of queueing call center models.

Second, the complexity of a real-life call center has been partly avoided in our study. Features such as abandonments, retrials, different types of the inbound calls, switching times between different tasks, and the finite number of back office tasks, are important but including them considerably complicate the analysis.

Introduction

In the context of highly congested call centers, the use of an alternative service channel can be proposed to customers so as to balance workload and avoid excessive abandonments. Such an alternative channel could be a Web site, an e-mail service, a proposition to the customer to call at a less busy time, a possibility of routing to another team of agents (see [START_REF] Manitz | The impact of a waiting-time threshold in overflow systems with impatient customers[END_REF]) or to leave a number and be called back later. We focus on this last alternative. Armony and Maglaras (2004a,b) propose a model in which customers are given a choice of whether to wait online for their call to be answered or to leave a number and be called back within a specified time.

They show that this callback scheme allows the system to increase its performance.

Combining inbound and outbound calls leads to a call center with blended operations. The key distinction of problems with blending comes from the fact that callbacks have less urgency relative to inbound calls. The call blending problem has led to research on performance evaluation (see [START_REF] Bernett | Blended call center performance analysis[END_REF], [START_REF] Pichitlamken | Modeling and simulation of a telephone call center[END_REF] and [START_REF] Deslauriers | Markov chain models of a telephone call center with call blending[END_REF]) and analysis of blending policies (see Gans et al. (2003), [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF], Armony and Maglaras (2004a) and Legros et al. (2013a)). Because of the strict waiting time requirement on inbound calls, it is best to give them priority over the outbound calls. To maximize agent productivity it would be best to assign a call (inbound or outbound) to every idle agent when it is possible. The objection against this policy is that it could lead to long waiting times for future arrivals from inbound calls.

In [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF] and [START_REF] Gans | A call-routing problem with service-level constraints[END_REF], it is shown under certain assumptions that the optimal assignment policy is of the following form: outbound calls should only be scheduled when there are no waiting inbound calls and when the number of idle agents exceeds a certain threshold.

In the case of a callback option, this policy has to be completed. The blending models in the literature usually assume an infinite number of outbound tasks. With a callback option this number is finite and depends on all the parameters of the system. Thus this number should be considered as another threshold for the choice of the threshold in the agents reservation. Otherwise, it could lead to long waiting times for outbound calls.

In this chapter, we consider a call center modeling with a single customer type and a callback option. We develop a method based on Markov chains to evaluate its performance measures. We also provide an efficient routing policy based on the number of reserved agents for the inbound calls and the number of outbound calls. The overall objective is to reach a certain long term service level for the inbound calls, while minimizing the expected waiting time of the outbound ones. Few of our key findings related to the performance measures are highlighted next. We prove that the expected waiting time for the outbound calls is decreasing in the proportion of outbound calls and in the threshold of reserved agents. We also show that sensitivity of the performance measures is higher on the threshold of reserved agents than on a limit in the number of outbound calls. These findings provide an intuitive justification of the optimal curve relating the two thresholds.

The rest of the chapter is organized as follows. In Section 5.2 we present the modeling of a call center with a callback option. In Section 5.3 we evaluate the performance measures and study the impact of the main parameters. In Section 5.4, we numerically study the general form of the optimal curve. We finally give some concluding remarks and highlight future research.

Model

We consider a call center modeled as a multi-server queueing system with two types of jobs, inbound and outbound calls. The arrival process of inbound calls is assumed to be a homogeneous Poisson process with rate λ. Inbound calls arrive at a dedicated first come, first served (FCFS) queue with infinite capacity (queue 1). There are s identical, parallel servers (agents in call center parlance).

Each agent can handle both types of jobs. We assume that the service times of inbound or outbound calls are exponentially distributed with the same rate µ. Neither abandonment nor retrials are modeled.

Inbound calls are more important than outbound ones in the sense that the former request a

quasi-instantaneous answer (waiting time in the order of seconds or minutes), while the latter are more flexible and could be delayed for several hours. The possibility of the callback option is given in this purpose. The objective of the call center manager is to minimize the expected waiting time for the outbound calls while satisfying a strict constraint on the inbound calls expected waiting time.

The functioning we consider for the callback option and the related customers behavior is inspired by Armony and Maglaras (2004a) and described as follow. This option is proposed to a newly arrival only when her expected waiting time is too long (i.e. too many waiting calls in queue 1). We define a limit k (k ∈ N) in the number of waiting calls in queue 1. Upon her arrival, a customer can find three situations. If at least one agent is available then she is routed to one of the idle agents. Otherwise if the number of waiting calls in queue 1 is strictly lower than k, the callback option is not proposed and she waits in queue 1. If this number is higher than or equal to k the callback option is proposed. The customer reaction is assumed to be probabilistic. Directly upon arrival, she decides to accept the callback with a probability q or she decides to wait in queue 1 with a probability 1 -q (q ∈ [0, 1]). Note that the callback option is only proposed upon arrival.

If a callback is asked, the call arrives at a dedicated FCFS queue with infinite capacity (queue 2).

The calls in queue 1 (inbound calls) have a non preemptive priority over the calls in queue 2 (outbound calls). As mentioned in Chapter 4, [START_REF] Bhulai | A Queueing Model for Call Blending in Call Centers[END_REF] prove for a similar model that the optimal policy is a threshold policy with the priority given to inbound calls (some servers reserved for inbound calls). More concretely, the functioning of the call center under a threshold policy is as follows. Let us denote the threshold by c, 1 ≤ c ≤ s (Note that the case c = 0 can not be considered for stability reasons). When an agent becomes idle, she handles the call at the head of queue 1, if any. If not, the agent may either handle a call at the head of queue 2 if any, or she remains idle. If the number of idle agents (excluding her) is at least s -c, then the agent in question handles an outbound call (from queue 2). Otherwise, she remains idle. In other words, there are s -c agents that are reserved for inbound calls.

In this chapter, we propose a threshold policy which adjusts the threshold as a function of the number of outbound calls in queue 2, denoted by y (y ∈ N). We want to minimize the expected waiting time of the outbound calls, denoted by E(W 2 ) and respect a constraint of service level, denoted by w * 1 on the expected waiting time of the inbound calls, denoted by E(W 1 ). We are also interested for the analysis on the performance in terms of proportion of customers who asks for a callback, denoted by π. In summary, our optimization problem can be formulated as

   Minimize E(W 2 ) subject to E(W 1 ) ≤ w * 1 , (5.1)
where the decision variable is c(y). We provide a relation between the threshold c and the number of calls in queue 2 so as to build the optimal curve that answer Problem (5.1).

Performance Measures Results in the Case c(y) = c

In this section we simplify the model by assuming that c(y) = c in order to evaluate the main performance measures (E(W 1 ), E(W 2 ) and π). We first analyze the underlying Markov chain.

Then we determine closed form expressions for E(W 1 ) and π and a numerical method to exactly evaluate E(W 2 ). Second we examine the behavior of the performance measures as a function of k, q and c. Finally, we reconsider the modeling by adding a limit in the capacity of the second queue.

Let us define the random process {(x(t), y(t)), t ≥ 0} where x(t) and y(t) denote the number of calls in queue 1 or in service and the number of outbound calls at a given time t ≥ 0, respectively. We have x(t), y(t) ∈ {0, 1, 2, ...}, for t ≥ 0. Since call inter-arrival and service times are exponentially distributed, {(x(t), y(t)), t ≥ 0} is a Markov chain. We denote by p x,y the steady-state probability to be in state (x, y) (x, y ∈ N), and by a the ratio λ µ . Because of the priority for inbound calls we can find y > 0 if and only if x ≥ c. For 0 ≤ x < s + k and y ≥ 0 the transition rate from state (x, y) to state (x + 1, y) is λ. For x ≥ s + k and y ≥ 0 the transition rate from state (x, y) to state (x + 1, y) is (1 -q)λ. For x ≥ 1, x ̸ = c and y ≥ 0 the transition rate from state (x, y) to state (x -1, y) is min(x, s)µ. For x ≥ s + k and y ≥ 0 the transition rate from state (x, y) to state (x, y + 1) is qλ. For x = c and y ≥ 1 the transition rate from state (c, y) to state (c, y -1) is cµ.

The stability condition is a s < 1. We also have

P s+x = a x s x P s ,
for 0 ≤ x ≤ k and

P s+k+x = ((1 -q)a) x s x P s+k , for x ≥ 0. With cµ(P c -p c,0 ) = λq ∞ ∑ x=0
P s+k+x , we obtain

P c = a c c! p 0,0 1 -q c a s-c c! s! a k s k 1 1- a(1-q) s
.

Since the overall sum of the probabilities equals to one, we have

p 0,0 =    c-1 ∑ x=0 a x x! + a c c! 1 -q c a s-c c! s! a k s k 1 1- a(1-q) s ( s-c-1 ∑ x=0 a x c! (c + x)! + a s-c c! s! k-1 ∑ x=0 a x s x + a s-c c! s! a k s k 1 1 -a(1-q) s )    -1
.

(5.2)

Thus we have a close form formula for p x,0 for 0 ≤ x ≤ c. For x ≥ 0, we have (1 -q)λp s+k+x,0 + sµp s+k+x+2,0 = (λ+sµ)p s+k+x+1,0 . The associated homogeneous equation is sµx 2 -(λ+sµ)x+(1q)λ = 0. The solutions of this equation are

x 1 = λ+sµ+ √ (λ+sµ) 2 -4sµ(1-q)λ 2sµ = 1 2 ( 1 + a s + √ (1 + a s ) 2 -4(1-q)a s ) and x 2 = 1 2 ( 1 + a s - √ (1 + a s ) 2 -4(1-q)a s ) . Then p s+k+x,0 = α 1 x x 1 + α 2 x x 2 . We have p s+k,0 = α 1 + α 2 and because cµp c,1 = qλ ∞ ∑ x=0
p s+k+x,0 we also have cµp c,1 = qλ

( α 1 1-x 1 + α 2 1-x 2
) . Thus we define α 1 and α 2 as a function of p s+k,0 and p c,1 . Since we have cµp c,1 + (c + 1)µp c+1,0 = λp c,0 , p c+1,0 , λp c+x,0 + (c + x + 2)µp c+x+2,0 = (λ + (c + x + 1)µ)p c+x+1,0 for 0 ≤ x ≤ sc -2 and λp s+x,0 + sµp s+x+2,0 = (λ + sµ)p s+x+1,0 for 0 ≤ x ≤ s -2, we evaluate p c+x,0 as a function of p c,1 for x > 0. Finally with the expression of p s+k+x,0 , we can compute p c,1 as a function of p s+k,0 . The others probabilities p s+k+x,y for x ≥ 0 and y > 0 are solutions of (1 -q)λp s+k+x,y + sµp s+k+x+2,y + qλp s+k+x+1,y-1 = (λ + sµ)p s+k+x+1,y . Consequently, we have Proportion of called back. We have

p s+k+x,y = Q 1,y (x)x x 1 + Q 2,y (x)x x 2 in which Q 1,y (
π = q ∞ ∑ x=0 P s+k+x = qP c a s-c c! s! a k s k 1 1 -a(1-q) s = qP s+k 1 -a(1-q) s .
Expected waiting time for the first queue. We have λE(

W 1 ) = ∞ ∑ x=0 xP s+x . Then λE(W 1 ) = k-1 ∑ x=0 xP s+x + ∞ ∑ x=0 (x + k)P s+k+x = P s k-1 ∑ x=0 x a x s x + P s+k    k ( 1 -a(1-q) s ) + a(1-q) s ( 1 -a(1-q) s ) 2    .
Expected waiting time for the second queue. For c = s, the system is work conserving.

Since the overall system is an M/M/s queue, we have a closed form expression for the expected waiting time of the overall customers; E(W ). We evaluate the expected waiting time for the second queue with the relation;

πE(W 2 ) + (1 -π)E(W 1 ) = E(W ).
When c < s the system is not work-conserving. We can evaluate E(W 2 ) with the steady states 5.1 illustrate the evolution of the performance measures as a function of k, q and c. As a function of q, we observe that E(W 1 ) and E(W 2 ) are decreasing and convex and π is increasing and concave. As a function of k, E(W 1 ) and E(W 2 ) are increasing and concave and π is decreasing and convex. As a function of c, E(W 1 ) is increasing and convex, E(W 2 ) is decreasing and concave and π is increasing and convex. We observe a negative correlation between π and E(W 2 ). It means that the smallest the proportion of called backs is, the more they will have to wait. When the proportion of called back is small, the proportion of inbound calls is important. Those inbound calls benefit from reservation and priority. Thus they impact more negatively the performance of the outbound calls. Increasing q is similar to decreasing k, an increasing in the proportion of outbound calls and a decreasing in the expected waiting time of both inbound and outbound calls. We note for E(W 1 ) and E(W 2 ) that the convexity is more important in q than the concavity in k. Recall that q only depends on the behavior of the customers. This implies that a little increasing in q from a low value will impact more than from a high value. The parameter k is defined by the manager of the call center but can not be too high since the manager usually wants to avoid abandonment. The only parameter of control is c. We observe that E(W 1 ) and π are increasing in c. This is not surprising. As c increases the reservation for inbound calls decreases which implies an increasing of E(W 1 ). Thus the proportion of callbacks also increases and E(W 2 ) decreases. Control on the number of outbound tasks. In what follows we limit the number of outbound waiting tasks to K (limited capacity for queue 2), the parameter K could be another parameter of control for the manager. In other words the callback option is not proposed when the number of outbound calls is higher or equal than K. We propose a method to answer the optimization problem by choosing the optimal couple (c,K). The Markov chain and the evaluation of the performance measures are similar to the case with infinite capacity in queue 2. The difference is that the variable y is limited to values between 0 and K. Table 5.2 illustrates the numerical results. The intuitive observation in Table 5.2 is that the increasing of K induces a decreasing of E(W 1 ), an increasing of E(W 2 ) and an increasing of π because increasing K allows more space for the callbacks. We find values of K for which E(W 1 ) and E(W 2 ) are both decreasing in c (example with K = 1). Since the number of outbound calls is very limited, the reservation does not impact much, an increasing of c increases the overall performance of the system and both E(W 1 ) and E(W 2 ) benefit from the increasing of c. Thus the best solution when K is small is to have c = s. The limit in this model is the risk of abandon when the arrival rate increases.

General Case

In this section we consider the general case, using simulation experiments. The problem remains the same, but we do the optimization in finding the optimal curve, c(y). In the impact of the number of outbound calls in the choice of the threshold c. Second, we propose an intuitive construction of the optimal curve c(y) in some simple cases. Finally, we present the problem of the construction of the optimal curve c(y) via dynamic programming.

We examine the impact of the number of outbound calls (y) on the choice of the threshold (c) by considering a simple curve y(c). We define a limit in the number of outbound calls denoted by y * . We reserve s -c agents for inbound calls only if we have at most y * outbound tasks, otherwise we do not reserve any agent (c = s). Note that with y * = ∞ we have the same model as in Section 5.3 and with y * = 0 we have no reservation (c = s). Although this method is not optimal, it helps to understand the impact of a decision made on the number of outbound calls. In Table 5. This result is important, it means that for a given service level w * 1 for the inbound calls we should choose the highest possible value for c and adapt y * . For example if w * 1 = 0.83, we find in Table We start the construction with c = s, then we choose to decrease c to s -1 and start with y * = 1 and increase y * so as to decrease E(W 1 ). If for any values of y * (even y * = ∞) we have

E(W 1 ) > w *
1 , hence we choose to have y * = ∞ when c = s -1 and therefore we decrease c to s -2.

If we find a value of y * for which E(W 1 ) > w * 1 and for y * + 1 we stop there and choose this value of y * and c = s -1. We continue on the same way until we reach c = 1.

Example: Curve c(y) for w * 1 = 0.78 with λ = 1, µ = 0.2, s = 6, k = 3 and q = 80%. We have Using dynamic programming we define the problem of the search of the optimal curve. We denote by c 1 and c 2 the cost of a waiting customer in queue 1 and queue 2, respectively. We denote by V n (x, y) the expected costs over n steps. We have

V ′ n+1 (x, y) = c 1 (x -s) + + c 2 y (5.3) + λ λ + sµ
(1 x<k+s V n (x + 1, y) + 1 x≥k+s (qV n (x, y + 1) + (1 -q)V n (x + 1, y))

+ µ λ + sµ (min(x, s)V n (x -1, y) + (s -min(x, s))V n (x, y)) ,
and

V n+1 (x, y) =    V ′ n+1 (x, y) if x ≥ s min(V ′ n+1 (x, y), V ′ n+1 (x + 1, y -1)) if x < s
(5.4)

The long-term average optimal actions are a solution of the optimal equation (in vector notation)

T V = g + V .
Another way of obtaining them is through value iteration, by recursively defining V n+1 = T V n , for arbitrary V 0 . In the expression we see that it is optimal to schedule a type 2 task only if V (x + 1, y -1) < V (x, y). With this relation we can numerically build the optimal curve y(c).

In summary, we numerically find an efficient curve y(c) that answers Problem (5.1). This curve is defined by two parameters: c * (0 < c * ≤ s) and y * (y * ≥ 0). For c < c * , we do not do any reservation and y(c) = 0. For c = c * , y(c * ) = y * . For c > c * , y(c) = +∞. The main reason is that the sensitivity of the performance measures is more important in c than on the number of outbound calls. We still have to compare this intuitive result with the optimal curve that can be found through a dynamic programming approach.

Conclusions and Future Research

In this chapter we considered simple call center modeling with a callback option. We assumed a threshold policy on the reservation of agents for inbound calls and derived the main performance measures. The optimization problem is defined as minimizing the expected waiting time of the outbound calls while respecting a service level constraint on the inbound ones. We answered this problem by building a curve representing the relation between the threshold and the number of outbound tasks.

For future research, it might be interesting to evaluate via closed form expressions, the performance measures on the optimal curve. Moreover, the complexity of the customer behavior has been partly avoided. It is also interesting to include the feature of customer abandonment for the outbound calls. Another extension is to consider different service time distribution for inbound and outbound calls.

Chapter 6

Conclusions and Perspectives

In this chapter, we give general concluding remarks and present directions for future research. For further details, we refer the reader to the concluding sections of the previous chapters.

Conclusions

This thesis focused on operations management issues for SBR and multi-channel call centers. We investigated an important problem in the design and management of SBR call centers. The purpose of this study was to propose an intelligent architecture. We considered the context of call centers with unbalanced parameters. Under most cases of asymmetry, the well known existing architectures such as chaining lose their robustness. Thus we proposed a new call center architecture (single pooling) and demonstrated its efficiency. We showed that SP behaves well in most cases of asymmetry in the parameters. There are opportunities for managers of call centers to improve performance using this new architecture.

Next, we considered a blended call center with calls and emails. The call service is characterized by successive stages where one of them is a break for the agent. We focused on the optimization of the email routing given that calls have a non-preemptive priority over emails. Our objective was to maximize the throughput of emails subject to a constraint on the call waiting time. We developed a general framework with two probabilistic parameters for the email routing to agents.

One parameter controls the routing between calls, and the other does the control inside a call conversation. We have also considered 4 particular cases corresponding to the extreme values of the probabilistic parameters. For these routing models, we have derived various structural results and discussed the theoretical results in order to provide guidelines to call center managers. In particular, we proved for the optimal routing that all the time at least one of the two email routing parameters has an extreme value.

Next, we considered call centers with inbound calls and an infinite supply of non-preemptable outbound jobs. We proposed a scheduling policy, refereed to as ATP, where the objective is to do as much outbound as possible while satisfying a service level constraint on the call waiting time. In the real-life call center context with fluctuating call arrival rate, the assignment policy for outbound adapts itself to the current service level. We showed the efficiency of ATP by comparing its performance with those of other policies. One of the main advantage of ATP is its ability to quickly react when an important change in the arrival process happens and its ability to avoid inefficient states when the arrival rate remains constant.

Finally, we considered a call center modeling with a callback option. We assumed a threshold policy on the reservation of agents for inbound calls and derived the main performance measures.

The optimization problem is defined as minimizing the expected waiting time of the outbound calls while respecting a service level constraint on the inbound ones. We answered this problem by building the optimal curve representing the relation between the threshold and the number of outbound tasks.

Future Research

As detailed in the concluding sections of the previous chapters, several interesting areas of future research arise. In what follows, we point out some of these research directions.

For the architectures of chaining as well as for single pooling, there is a need to provide and analyse a much general model. Another interesting work is to find an architecture which resists better than chaining and single pooling to blocking, especially in small call centers.

In the modeling of Chapter 3, it would be interesting to extend the structural results to the multi-server case. It would also be useful but challenging to extend the analysis to cases with an additional channel, in particular the chat which is increasingly used in call centers. Using the chat channel, an agent may handle many customers at the same time, which represent an additional opportunity to efficiently use the agent time.

Future research on the subject of adapting blending in Chapter 4 may follow two directions.

First, further analytical analysis of the adaptive blending might be useful to better understand ATP and maybe to propose a better adaptive policy. Second, the complexity of call centers has been partly avoided in our study. The customers abandonment, the callback option, the different types of the inbound calls, the switching times between different tasks, the tiredness of the agent, a finite number of back office tasks could interfere in our results and might introduce other interesting performance measures in the optimization problem.

For future research on the callback option in Chapter 5, it may be useful to derive the optimal performance measures. Moreover, again, the complexity of the customer behavior has been partly avoided. The customers abandonment could occur even when the callback option is proposed. The service of the outbound calls could be done with a part of already informed customers which induces very short service times duration.

length up to the point that ensures the accuracy objective. This implies that the simulation length may vary from one set of parameters to another. The total number of generated calls varies and is in the order of tens of millions. The delay to run a simulation also varies and is in the order of several minutes.

A.1.1 Single Pooling

In what follows we present three staffing heuristics, and then select the best one. The heuristics consist on adaptations of greedy and local search algorithms.

Algorithm 1: Decreasing Greedy in Team 0

Without customers 0, single pooling is simply an FD model. A first idea of staffing is then to use a decreasing greedy algorithm as follows. We start such that we have a collection of n + 1 independent M/M/s queues. In each team i, the number of agents is the minimum required one to reach W i ≤ W * i = 0.2, for i = 0, 1, ..., n. In each iteration, we decrement the number of agents in team 0 by one, and evaluate all W i , for i = 0, 1, ..., n. We stop the algorithm when all the service levels are no longer reached for the first time. We then consider the results of the before last iteration. Table A.1 presents the results of the decreasing greedy algorithm in a single pooling model with 3 customer types, and compare it with the those of FF and FD models. The staffing level in team i is denoted by s i , i = 0, 1, 2. Algorithm 2: Increasing Greedy

= µ 0 = 0.2, W * i = W * 0 = 0.2, i = 1, 2)
Another idea is to proceed by introducing customers 0 in the system step by step. We start from an FD model with no customers 0, and we define the staffing level in team i such that W i ≤ 0.2, i = 1, ..., n (team 0 is being empty). In each iteration, we increase λ 0 by a given small step value (we have chosen in the experiments a sufficiently small step of λ 0 /100). If W 0 ≥ W * 0 , we add one agent in team 0. If W 0 ≤ W * 0 and W i > W * i for some customer types, we then add an agent in the team with the highest W i . We stop the algorithm once λ 0 reaches its value and the constraints W i ≤ W * i are all satisfied, for i = 0, 1, ..., n. Algorithm 3: Increasing Greedy with No Agents in Team 0

= µ 0 = 0.2, W * i = W * 0 = 0.2, i = 1, 2)
The algorithm is identical to the previous one, expect that we force team 0 to be empty. Table A .3 presents the simulated results for this algorithm.

From the results of all algorithms, we observe that the decreasing greedy algorithm (algorithm 1) is the worst. The reason is that it is not possible to increase or decrease the number of agents in a regular team i (i = 1, ..., n). Many effective configurations could not then be reached under this algorithm. The other two algorithms are equivalent in terms of the total number of agents in our simulation experiments. We have chosen to use algorithm 2 in the experiments of Section 5 of the main paper. We go further in order to check the quality of algorithm 2. In Table A.4, we provide optimization results using algorithm 2 and also using other configurations with one agent in less. We observe that these other configurations do not allow to satisfy all the constraints W i ≤ W * i = 0.2, for i = 0, 1, ..., n, which proves the efficiency of algorithm 2.

= µ 0 = 0.2, W * i = W * 0 = 0.2, i = 1, 2) s1 s2 s0 W1 W2 W0 Constraints λ1 =

A.1.2 Chaining

We also use a greedy algorithm to optimize the staffing of chaining. The simulation results reveal that increasing and decreasing greedy algorithms are efficient and very similar if we start the optimization heuristic with a good initialization of the team sizes. We choose to use a decreasing greedy algorithm since it is faster than an increasing greedy one (no need to increase the λ 0 with

= W * i = 0.2, ∑ 4 i=0 si = 48, i = 1, ..., 4, p = 25%, U = V = 1) Ϭ͕ϱϬ Ϭ͕ϱϱ Ϭ͕ϲϬ Ϭ͕ϲϱ Ϭ͕ϳϬ Ϭ͕ϳϱ Ϭ͕ϴϬ Ϭ͕ϴϱ Ϭ͕ϵϬ Ϭ͕ϵϱ ϭ͕ϬϬ ϭ Ϯ ϯ ϰ ϱ ĐŚĂŝŶŝŶŐ ^W s t (c) Impact of V (λ0 = 2, W * 0 = W * i = 0.2, ∑ 4 i=0 si = 44, i = 1, ..., 4, p = 25%, p ′ = 20%, U = 1) Ϭ͕ϭ Ϭ͕ϭϱ Ϭ͕Ϯ Ϭ͕Ϯϱ Ϭ͕ϯ Ϭ͕ϯϱ Ϭ͕ϰ ϭ Ϯ ϯ ϰ ϱ ĐŚĂŝŶŝŶŐ ^W h t (d) Impact of U (λ0 = 2, λi = 1.5, W * 0 = W * i = 0.2, ∑ 4 i=0 si = 48, µ0 = 0.2, i = 1, ..., 4, p = 25%, p ′ = 20%, V = 1)

A.3 Impact of Abandonment

The experiments of Tables A.6-A.9 are associated to Figures 2.10 We consider the performance measures of the average waiting time (denoted by W q ), the average service time (denoted by W s ) and the overall time spent in the call center (W q + W s ) for a task.

In the case s = 1 we can find formulas for the performance measures. The Markov chain in this case is the one of a M/M/1 queue. Only the performance measures are evaluated differently. If we denote by p x (x ∈ N) the probability to have x tasks in the call center, we have p x = a x (1 -a)

with a = λ/µ. Thus with Little law and a capacity of simultaneity of k, we have

W q = 1 λ +∞ ∑ x=0 xp x+k = 1 λ (1 -a)a k+1 1 (1 -a) 2 = 1 λ a k+1 1 -a , and
W s = 1 λ ( k-1 ∑ x=0 xp x + k k-1 ∑ x=0 p x+k ) = 1 λ a(1 -a k ) 1 -a .
We note that W q + W s = 1 λ a 1-a is the overall spent time in a classical M/M/1 queue. Thus is we increase the level of simultaneity k, we decrease the expected waiting time but we increase the service time. There is no benefit in the overall spent time in the call center.

In Table B.1 we present numerical results for the evolution of the performance measures in function of the level of simultaneity. We use simulation for s > 1. When s > 1 we still observe a decreasing of the average waiting time and an increasing of the service time in function of k.

Although the overall spent time increases with k, then the call center is less efficient when the level of simultaneity increases. The reason is the possibility of inefficient states in the system when s > 1; for example an agent with two tasks and one idle agent at the same time.

These observations reduce the possible interest for a simultaneous treatment. Note that the conclusions are based on the assumption that there is no increasing in productivity when an agent treats tasks simultaneously. That is why we also consider another simple modeling similar as the previous one. The only difference is when an agent treats n tasks simultaneously, the service rate is µ (1+t) n-1 for 1 ≤ n ≤ k and t ≥ 0. This model is more compatible with human behavior. In this model the productivity is better for a simultaneous treatment than for a successive one when the number of simultaneous tasks is smaller than a limit (determined with t). After this limit the productivity is worse. Roughly speaking the agent is overwhelmed. In Table B.2 we present the simulated results for different values of t. We observe an opportunity in working simultaneously on different tasks when t is small. For t = 0.1, the average waiting time decreases until k = 10 and the overall spent time also decreases until k = 2. For t = 0.5 or t = 1 there is no benefit in simultaneity for the overall spent time and neither for the average waiting time when t = 1. This implies that simultaneity in the treatments is an opportunity only if working simultaneously on more than one task induces a huge increase in productivity (t small). Moreover it is necessary to control the number of possible tasks done simultaneously (k) to avoid important decrease in the performance measures.

B.2 Computation of P (W < t)

This section is related to Sections 2.5.1 and 2.5.2. We give details on the derivation of the quantities P(W < t, (A, n)), P(W < t, (B, n)), P(W < t, (B ′ , n)), P(W < t, (C, n)) and P(W < t, (M, n)), Table B.2: Performance Measures in function of the level of simultaneity (λ = 1,µ = 0.2, s = 6) Consider a distinct sums of exponential random variables with rates µ 1 , µ 2 , ..., µ a and a number of terms in each sum equals to r 1 , r 2 , ..., r a , respectively. Let F (t), for t ≥ 0, denote the cdf of the summation of all the random variables. From Amari and Misra (1997), we have

t = 0.1 t = 0.5 t = 1 k W q W s W q + W s W q W s W q + W s W q W s W q + W s 1 2
F (t) = 1 -   a ∏ j=1 µ r j j   a ∑ k=1 r k ∑ l=1 Ψ k,l (-µ k )t r k -1 exp(-µ k t) (r k -l)!(l -1)! , (B.1) for t ≥ 0, with Ψ k,l (x) = - ∂ l-1 ∂x l-1   a ∏ j=1,j̸ =k (µ j + x) -r j   ,
for x ∈ R. For example for Model 1, we obtain

P(W < t, (A, n -1)) = 1 -(µ 1 µ 2 µ 3 ) n 3 ∑ k=1 n ∑ l=1 Ψ k,l (-µ k )t n-1 exp(-µ k t) (n -l)!(l -1)! ,
for n ≥ 1. Using the Leibnitz formulae, we have

Ψ k,l (x) = (-1) l (l -1)! l-1 ∑ i=0 ( n i )( n l -1 -i ) (µ ′ k + x) -(n+i) (µ ′′ k + x) -(n+l-1-i) , with µ ′ k , µ ′′ k ∈ {µ 1 , µ 2 , µ 3 \ µ k } and µ ′ k ̸ = µ ′′ k . We then deduce, for n ≥ 1, that P(W < t, (A, n -1)) = 1 -(µ 1 µ 2 µ 3 ) n n ∑ l=1 l-1 ∑ i=0 ( (-1) l ( n i )( n l-1-i ) (µ 3 -µ 1 ) -(n+i) (µ 2 -µ 1 ) -(n+l-1-i) t n-1 exp(-µ 1 t) (n -l)! + (-1) l ( n i )( n l-1-i ) (µ 3 -µ 2 ) -(n+i) (µ 1 -µ 2 ) -(n+l-1-i) t n-1 exp(-µ 2 t) (n -l)! + (-1) l ( n i )( n l-1-i ) (µ 1 -µ 3 ) -(n+i) (µ 2 -µ 3 ) -(n+l-1-i) t n-1 exp(-µ 3 t) (n -l)!
) .

For Model PM, let us take for example the case of a new call that arrives to a system with n (n ≥ 1) calls (excluding the new call): n -1 calls are waiting in the queue, and assume that the one in service is being in the first stage of service. The waiting time in the queue of the new call is the time it takes to clean the system from the n customers ahead of her. Recall that in Model PM, an agent works on emails during the second stage of a call service with probability q. Then,

( n k ) q k (1 -q) n-k , for 0 ≤ k ≤ n and n ≥ 1.
The cdf of the waiting time of our new call corresponds to the quantity P(W < t, (A, n -1))

(n ≥ 1). Using Equation B.1, we obtain

P (A,n-1) (W < t) = 1 -(µ 1 µ 2 µ 3 ) n n ∑ k=0 ( n k ) q k (1 -q) n-k µ k 0 × ( 3 ∑ i=1 n ∑ l=1 Ψ i,l (-µ i )t n-1 exp(-µ i t) (n -l)!(l -1)! + k ∑ l=1 Ψ 0,l (-µ 0 )t k-1 exp(-µ 0 t) (k -l)!(l -1)! ) ,
for n ≥ 1, with only cases when at least we find i, j = 1, 2, 3 which verify µ i ̸ = µ j . We present here the most general case is when µ i ̸ = µ j for i, j ∈ {1, 2, 3}. The other cases can be easily derived from this method. The sub-generator matrix has three different eigenvalues: -µ 1 , -µ 2 and -µ 3 . Each of these eigenvalues is associated to an eigenspace of dimension 1. The eigenvector associated with

-µ 1 is v 1 = (1, 0, . . . , 0). The eigenvector associated with -µ 2 is v 2 = (v 2,1 , v 2,2 , . . . , v 2,3n ) for which v 2,i = ( -µ 2 -µ 1 µ 1 ) i-1 for i ∈ {1, . . . , n + 1} and v 2,i = 0 else. The eigenvector associated with -µ 3 is v 3 = (v 3,1 , v 3,2 , . . . , v 3,3n ) for which v 3,i = ( -µ 3 -µ 1 µ 1 ) i-1 for i ∈ {1, . . . , n + 1}, v 3,i = ( -µ 3 -µ 1 µ 1 ) n ( -µ 3 -µ 2 µ 2 ) i-(n+1)
for i ∈ {n + 1, . . . , 2n + 1} and v 2,i = 0 else. If we denote by e 1 , e 2 , ..., e 3n the usual basis B. In the new basis B ′ in which we change e 1 with v 1 , e n+1 with v 2 and e 2n+1 with v 3 and keep the other vector. We denote by P the transition matrix from basis B to basis B ′ . We can write the sub generator matrix in this new basis as

     E 1 O O O E 2 O O O E 3      , (B.2) in which E i (i = 1, 2, 3) is a sub-matrix of dimension n × n equals to                -µ i µ i 0 . . . 0 0 0 -µ i µ i . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . -µ i µ i 0 0 0 . . . 0 -µ i µ i 0 0 . . . 0 0 -µ i                , (B.3)
and O is as sub-matrix of dimension n × n only composed with 0. We can decompose E i into a sum of a diagonal matrix (denoted by D i ) and a nilpotent one (denoted by N i ). We have

D i =                -µ i 0 0 . . . 0 0 0 -µ i 0 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . -µ i 0 0 0 0 . . . 0 -µ i 0 0 0 . . . 0 0 -m i                , (B.4)
and

N i =                0 µ i 0 . . . 0 0 0 0 µ i . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 µ i 0 0 0 . . . 0 0 µ i 0 0 . . . 0 0 0                . (B.5)
We denote by D the matrix

     D 1 O O O D 2 O O O D 3      , (B.6)
and by N the matrix and for i = 1, 2, 3

     N 1 O O O N 2 O O O N 3      . ( B
e xD i =                e -xµ i 0 0 . . . 0
0 . . . 0 0 e -xµ i                , (
e xN i =                1 xµ i x 2
2 µ 2 i . . . x n-2 (n-2)! µ n-2 i x n-1

(n-1)! µ n-1 i 0 1 xµ 1 . . . Then for i=1,2,3 we have We need to evaluate e xM = P -1 e x(N +D) P . In other words we need to write the matrix e x(N +D) in the basis B. In fact the only differences between the basis B and the basis B ′ are the n + 1 and 2n + 1 vectors. We denote by f the linear application associated to the matrix e x(N +D) . The target is to calculate 3 )z + c -ab 3 + 2a 3 27 = 0. Define now p and q as p = b -a 2 3 and q = c -ab 3 + 2a 3 27 .

e x(N i +D i ) =                e -
Then, Equation (B.22) becomes z 3 + pz + q = 0. (B.23)

The discriminant denoted by ∆ of this equation is ∆ = q 2 + 4 27 p 3 . Let us define u, v and j by u = ( 1 2 (-q + √ ∆)) 1/3 , v = ( 1 2 (-q -√ ∆)) 1/3 and j = e i 2π 3 . Three situations are possible for ∆:

∆ > 0, ∆ < 0 or ∆ = 0.

• If ∆ > 0, Equation (B.23) has one real solution z 1 = u + v and two complex solutions z 2 = ju + jv and the conjugate of z 2 denoted by z 2 .

• If ∆ < 0, then v = u and Equation (B.23) has three real solutions: z 1 = u + u, z 2 = ju + ju and z 3 = j 2 u + j 2 u.

• If ∆ = 0, Equation (B.23) has two real solutions: a simple one z 1 and a double one z 2 .

Solution of a cubic equation: Consider the general cubic equation written as . This is equivalent to (1 -ρ 1 -ρ 2 -ρ 3 -qρ 0 ) (1 + ρ 0 ) -ρ 0 (1 -p) (1 + pρ 0 ) > 0, or also (p 2 -p -q)ρ 2 0 + (p -q -ρ 1 -ρ 2 -ρ 3 )ρ 0 + 1 -ρ 1 -ρ 2 -ρ 3 > 0.

(B.29)

The discriminant for this inequality is ∆ = (p-q -ρ 1 -ρ 2 -ρ 3 ) 2 -4(p 2 -p-q)(1-ρ 1 -ρ 2 -ρ 3 ) > 0.

Equation (B.29) has then the two following solutions:

√ (p -q -ρ 1 -ρ 2 -ρ 3 ) 2 -4(p 2 -p -q)(1 -ρ 1 -ρ 2 -ρ 3 ) -q + p -(ρ 1 + ρ 2 + ρ 3 ) 2(q -p 2 + p)
and -√ (p -q -ρ 1 -ρ 2 -ρ 3 ) 2 -4(p 2 -p -q)(1 -ρ 1 -ρ 2 -ρ 3 ) + q -p + ρ 1 + ρ 2 + ρ 3 2(q -p 2 + p) .

Since the first solution is positive (denoted by ρ 0 ) and the second one is negative, ρ 0 ∈ [0; ρ * 0 ), which finishes the proof of the lemma. ✷ for 0 ≤ u ≤ s. Thus T (s, u, a) = µ ( a + 1 gu ) -λ, with

g u = 1 u + a u(u + 1) + • • • + a s-u-1 u(u + 1)(u + 2) • • • (s -1)
+ a s-u u(u + 1)(u + 2) • • • (s -1)(s -a) , for 0 < u ≤ s (and T (s, 0, a) = 0). We may write for 0 < u < s

g u+1 -g u = ( 1 u + 1 - 1 u
) + ( a (u + 1)(u + 2) -a u(u + 1)

) + • • • (C.1) + ( a s-u-1 (u + 1)(u + 2) • • • (s -1)(s -a) - a s-u-1 u(u + 1) • • • (s -1) ) + ( - a s-u u(u + 1) • • • (s -1)(s -a)
)

.

Since each term of the summation in the right hand of Equation (C.1) is strictly negative, g u+1 < g u for 0 < u < s. Then, g u is strictly decreasing in u for 0 < u ≤ s. We also have T (s, 1, a) > 0 =

T (s, 0, a). This implies that T (s, u, a) is strictly increasing in u, for 0 ≤ u ≤ s. (u+1+k)! . This implies that f u < f u+1 , for 0 ≤ u < s. Then, P (W < τ ) is strictly decreasing in u, for 0 ≤ u ≤ s.

We next focus on the proof of convexity of f u in u (for s ≥ 2). We prove that f u +f u+2 -2f u+1 > 0 for 0 ≤ u ≤ s -2. Since

f u + f u+2 -2f u+1 = 2f u f u+1 f u+2 ( 1 2 1 f u+1 ( 1 fu + 1 f u+2 ) -1 fu 1 f u+2 )
, it suffices to show that 1 2 1 f u+1 ( 1

( 1 fu + 1 f u+2 ) -1 fu 1 f u+2 > 0, for 0 ≤ u ≤ s -2. We have 1 2 1 f u+1 ( 1 f u + 1 f u+2 ) - 1 f u 1 f u+2 =
f u + 1 f u+2 ) - 1 f u 1 f u+2 = 1 2 a s!(s -a) a u-s ( - 1 u! + a (u + 1)! ) + a 2u-2s 2 ( a u!(u + 1)! + ( a (u + 1)! ) 2 - 1 u! s-u ∑ k=2 a k (u + k)! ) ,
for 0 ≤ u ≤ s -2. We thus need to show that a s!(s -a)

( -1 + a u + 1 )
+ a u-s ( a (u + 1)! + a 2 (u + 1)!(u + 1) -

s-u ∑ k=2 a k (u + k)! ) > 0. (C.2)
for 0 ≤ u ≤ s -2.

Case 1: a ≥ u + 1. We have The result then follows.

Case 2: a < 1. Consider the two last terms in the right hand side of Equation (C.2). We may write a 2 (u + 1)!(u + 1) -

( s-u ∑ k=2 a k (u + k)! ) (C.3) = a 2 (u + 1)! ( 1 u + 1 - 1 u + 2 - a (u + 2)(u + 3) -• • • - a s-u-2 (u + 2)(u + 3)(u + 4) • • • s ) > a 2 (u + 1)! ( 1 u + 1 - 1 u + 2 - 1 (u + 2)(u + 3) -• • • - 1 (u + 2)(u + 3)(u + 4) • • • s ) .
for 0 ≤ u ≤ s -2. We next consider the last line of Equation (C.3) and prove the following statement:

P (k) : 1 u + 1 - 1 u + 2 - 1 (u + 2)(u + 3) -• • • - 1 (u + 2)(u + 3)(u + 4) • • • (u + k) = Q k (u) (u + 1)(u + 2)(u + 3)(u + 4) • • • (u + k) ,
for k ≥ 2, where Q k (u) is a polynomial in u, for 0 ≤ u ≤ s -2, with all coefficients higher than or equal to 1. For k = 2, we have 1 u+1 -1 u+2 = 1 (u+1)(u+2) . So, Q 1 (u) = 1 and P (1) is true. Assume now that P (k) is true for k ≥ 2. We have

1 u + 1 - 1 u + 2 - 1 (u + 2)(u + 3) -• • • - 1 (u + 2)(u + 3)(u + 4) • • • (u + k) (C.4) - 1 (u + 2)(u + 3)(u + 4) -• • • - 1 (u + 2)(u + 3)(u + 4) • • • (u + k + 1) = Q k (u) (u + 1)(u + 2)(u + 3)(u + 4) • • • (u + k) - 1 (u + 2)(u + 3)(u + 4) • • • (u + k + 1) = (u + k + 1)Q k (u) -(u + 1) (u + 1)(u + 2)(u + 3)(u + 4) • • • (u + k)(u + k + 1)
.

Then Q k+1 (u) = (u + k + 1)Q k (u) -(u + 1). Using in addition that k ≥ 2 and the fact that all the coefficients of Q k (u) are higher than or equal to one, we deduce that P (k + 1) is also true. As a conclusion, the statement P (k) is true for all k ≥ 2. Using the statement P (k), we state that the last line of Equation (C.3) is strictly positive. Then a 2 (u + 1)!(u + 1) - + a u-s 1 (u + 1)! > 1 s!(s -1)

( s-u ∑ k=2 a k (u + k)! ) > 0
( -1 + a u + 1 )
+ a u-s 1 (u + 1)! = 1 s!(s -1)

a u + 1 + a u-s (u + 1)! - 1 s!(s -1)
.

for 0 ≤ u ≤ s -2. Since u ≤ s -2 and 0 ≤ a < 1, a u-s > 1 and (u + 1)! < s!(s -1). Therefore a u-s (u + 1)! -1 s!(s -1) > 0. We then obtain the convexity result of f u for 0 ≤ u ≤ s -2 and a < 1.

This completes the proof of the proposition. ✷

From an extensive numerical study, we observe for both cases (equal or unequal service rates) that P (W < τ ) is also concave in u for 1 ≤ a < u + 1. Also, the results of Proposition 5 still hold

C.2 Evaluation of The ATP

In Table C.1 we propose scenarios to compare the throughput found with the ATP method and the optimal throughput. We consider various situations by changing the service rates, the workload on calls and the size of the call center. abandonment,5,6,8,17,30,[39][40][41]47,60,94,117,119,120,126,132,135,152 ACD,2,7,91,117 approximations,6,[23][24][25]28,29,56,57,67,85,87,89,95 arrival process,20,23,59,65,94,113,114,117,120,134,146 ATP,[106][107][108][109][110][111][112][113][114]116,117 back office tasks,54,117,135,97 blending,4,53,57,92,93,117,119,135 blocking,[35][36][37]46,47,49,135 callbacks,5,7,58,[118][119][120][121][122][123]126,127,132,134,135 chaining,3,6,10,12,14,15,[18][19][20]31,32 distribution Erlang,58,66,152 exponential,58,60,63,120,122,146 hypoexponential,152 dynamic programming,129 expected waiting time,7,9,22,27,29,42,58,68,69,77,78,118,120 FCFS,14,19,59,94,120,121,152 idle agent,19,21,59,60,92,95,96,119,121,147 inbound,4,6,7,53,54,[57][58][59]94,99,[117][118][119][120]122,126,129,132,134,135 Markov chain,6,20,21,23,[56][57][58][62][63][64]66,70,[72][73][74]93,100,117,120,122,123,127,128,[1][2][3]5,6,91,93,1,3,6,8,19,52,55,57,90,95,117 optimal curve,7,122 optimization problem,52,56,57,60,75,83,85,89,90,118,127,132,134,135 outbound,4,7,53,54,57,59,[91][92][93][94][95][117][118][119][120][121][122]126,127,129,132,134,135

Index

Résumé

Les centres d'appels connaissent un grand succès depuis leur introduction dans les entreprises de service. Ils sont le principal point de contact avec les clients, et une composante essentielle de la majorité des entreprises. L'émergence des centres d'appels à grande échelle a suscité de nombreuses problématiques de management. Dans cette thèse, nous considérons des problématiques de management orientées sur les centres multi-canaux et multi-compétences. L'objectif de notre travail est de trouver des résultats qualitatifs et quantitatifs utiles pour le management.

Dans la première partie, nous considérons les architectures de centres multi-compétences à flexibilité limitée.

Abstract

Call centers have been introduced with great success by many service-oriented companies. They become the main point of contact with the customer, and an integral part of the majority of corporations. The large-scale emergence of call centers has created a fertile source of management issues. In this PhD thesis, we focus on various operations management issues of multi-skill and multichannel call centers. The objective of our work is to derive, both qualitative and quantitative, results for practical management.

In the first part, we focus on architectures with limited flexibility for multi-skill call centers. The context is that of call centers with asymmetric parameters: unbalanced workload, different service requirements, a predominant customer type, unbalanced abandonments and high costs of crosstraining. The most knowing architectures with limited flexibility such as chaining fail against such asymmetry. We propose a new architecture referred to as single pooling with only two skills per agent and we demonstrate its efficiency under various situations of asymmetry.

In the second part, we focus on routing problems in multi-channel call centers. In the first study, we consider a blended call center with calls arriving over time and an infinitely backlogged queue of emails. The call service is characterized by three successive stages where the second one is a break. We focus on the optimization of the email routing to agents. The objective is to maximize the throughput of emails subject to a constraint on the call waiting time. Various guidelines to call center managers are provided. In particular, we prove for the optimal routing that all the time at least one of the two email routing parameters has an extreme value.

In the second study, we examine a threshold policy on the reservation of agents for the inbound calls. We study a general non-stationary model where the call arrival follows a non-homogeneous Poisson process. The optimization problem consists on maximizing the throughput of outbound tasks under a constraint on the waiting time of inbound calls. We propose an efficient adaptive threshold policy easy to implement. This scheduling policy is evaluated through a comparison with the optimal performance measures found in the case of a constant stationary arrival rate, and also a comparison with other intuitive adaptive threshold policies in the general non-stationary case.

In the third study, we consider a call center model with a call back option, which allows to transform an inbound call into an outbound one. The optimization problem consists on minimizing the expected waiting time of the outbound calls while respecting a service level constraint on the inbound ones. We propose a routing policy with two thresholds, one on the reservation of the agents for inbound calls, and another on the number of waiting outbound calls. A curve relating the two thresholds is determined.

Keywords

  Figure 2.1: Call center configurations

Figure 2

 2 Figure 2.4: Impact of p on SP performance (n = 4,µ i = µ 0 = 0.2, ∑ 2 i=0 λ i = 8, λ i = λ j , s i = 12 for i, j = 1,..., 4)

Figure

  Figure 2.7: Preference zone

Figure 2 . 8 :

 28 Figure 2.8: Preference zone (µ 0 = 0.2, λ 0 = 4, λ i = 1 for i = 1, ..., 4, ∑ 4 i=0 1

Figure 2 . 9 :

 29 Figure 2.9: Preference zone (λ 0 = 4,λ i = 1, µ i = µ 0 = 0.2 and W * i = 0.2 for i = 1, ..., 4, p = 50%, p ′ = 20%, U = V = 1, n = 4)

Figure 2 .

 2 Figure 2.10: Impact of abandonment

  Impact of p (µi = µ0 = 0.2, W * 0 = W * i = 0.2, ∑ 4 i=0 λi = 8, i = 1, ..., 4, p ′ = 20%, U = V = 1, n = 4) 25, i = 1, ..., 4, p = 50%, U = V = 1, n = 4)

Figure 2 .

 2 Figure 2.11: Impact of the asymmetry in abandonment

Figure 2

 2 Figure 2.12: Preference zone (µ i = µ 0 = 0.2, W * i = W * 0 = 0.2 for i = 1, ..., 4, p ′ = 20%, U = V = 1, n = 4)

  Impact of U and p (λi = λj for

  Figure 2.14: Preference zone

Figure 3 . 1 :

 31 Figure 3.1: Markov chain for Model PM

  as shown in Section B.2 of the appendix. Fortunately, the computation of the probabilities a n , b n , b ′ n , c n and m n , for n ≥ 0, simplifies further. In what follows, we avoid the matrix analysis (as developed in Section 3.4.1) by providing for each model the expression of each one of the latter probabilities as a function of u n , where u n = a n + b n + b ′ n + c n , for n ≥ 0.

  and σ 3 = ρ 1 ρ 2 ρ 3 . The cubic equation associated to the recurrent Equation (3.11) is x 3 = (σ 1 + σ 2 + σ 3 ) x 2 -(σ 2 + 2σ 3 ) x + σ 3 , with x as variable. It remains now to apply the Cardan-Ferrari method in order compute u n , for n ≥ 0. The type of the solutions depends on the discriminant ∆ of the cubic equation (see Section B.4 of the Appendix).

  Figure 3.6: Comparison between Models 1 to 4 with a constraint on E(W )

  Figure 3.7: Comparison between Models 1 to 4 with the constraint on P(W < 1)) ≥ α *

  .21) for p, q ∈ [0, 1]. Let us denote the right hand sides of Inequalities (3.20) and (3.21) by the functions in p ∈ [0, 1] f (p) and g(p), respectively. Illustrations of these functions are given in Figure 3.8.

  Figure 3.8: Behavior of f (p) and g(p)

Figure 3 . 9 :

 39 Figure3.9: Optimal p and q with w * = 1

  Figure3.10: Optimal p and q with P(W < 1) ≥ 80%)

  [START_REF] Deslauriers | Markov chain models of a telephone call center with call blending[END_REF] extend the earlier mentioned papers by having different types of agents. outbound calls are served only by blend agents, whereas inbound calls can be served by either inbound-only or blend agents. Inbound callers may balk or abandon. They evaluate several performance measures of interest, including the rate of outbound calls and the proportion of inbound calls waiting more than some fixed number of seconds. A collection of continuous-time Markov chain (CTMC) models which capture many real-world characteristics while maintaining parsimony that results in fast computation are presented. They discuss and explore the tradeoffs between model fidelity and efficacy and compare different CTMC models with a realistic simulation model of a Bell Canada call center.[START_REF] Armony | Fair Dynamic Routing in Large-Scale Heterogeneous-Server Systems[END_REF] present an optimization problem; the objective is to minimize steadystate expected customer wait time subject to a fairness constraint on the workload division. They

  s-1 for x ≥ 0. For y = s -1 and x > 0, we may write (λ+µ+(s-1)µ 0 )p x,s-1 = λp x-1,s-1 +µp x+1,s-1 +sµ 0 p x+1,s . The associated homogeneous equation in the variable z is µz 2 -(λ + µ + (s -1)µ 0 )z + λ = 0, for z ∈ C. It has two solutions denoted by z 1 and z 2 and are given by z 1 = 1 2µ

  we can easily derive A k,k and B k,k as a function of p 0,k and p 1,k , for 0 ≤ k < s. The relation between p 1,k and p 0,k is given by

Figure 4 . 1 :

 41 Figure 4.1: Evolution of the Sensitivity of the Service Level function of the Threshold and the Workload (s = 10, τ = 30 seconds, µ = µ 0 = 0.2)

Figure 4 . 2 :

 42 Figure 4.2: Evolution of the threshold, the service level and the throughput (Scenario 2)

FigureFigure 5

 5 Figure 5.1: Markov chain for the callback option, c(y) = c

  x) and Q 2,y (x) are polynomials in the variable x with a degree equals to y. We evaluate the first coefficients of those polynomials with the ones of Q 1,y-1 (x) and Q 2,y-2 (x). The constant term in Q 1,y (x) and in Q 1,y (x) is found as a function of p c,y+1 with the relation cµp c,y+1 = qλ ∞ ∑ x=0 p s+k+x,y . Then we compute p c,y+1 using the relations λp c+x,y + (c + x + 2)µp c+x+2,y = (λ + (c + x + 1)µ)p c+x+1,y for 0 ≤ x ≤ s -c -2 and λp s+x,y + sµp s+x+2,y = (λ + sµ)p s+x+1,y for 0 ≤ x ≤ s -2.

Figure 5

 5 Figure5.2: E(W 1 ), E(W 2 ) and π as a function of q with λ = 1, s = 6, c = 6 and k = 3

Figure 5

 5 Figure 5.3: E(W 1 ), E(W 2 ) and π as a function of k with λ = 1, s = 6, c = 6 and q = 80%

Figure 5

 5 Figure 5.4: E(W 1 ) and π as a function of c with λ = 1, s = 6, k = 3 and q = 80%

  3 we present simulation results for different values of c and y * . The number of states for which it is possible to have waiting outbound calls increases in y * . Thus we observe that for a given value of c, E(W 1 ) decreases and E(W 2 ) increases in y * . Since the number of counterproductive states (with y > 0 and x < s) also increases in y * , we observe that E(W ) (expected average waiting time of both inbound and outbound calls) increases in y * . When c increases the performance measures are less sensitive to the increasing of y * . In particular, the extreme case y * = 0 is equivalent to c = s.

  Figure 5.5: Markov chain with a general curve c(y)

  the following simulated results. For c = 5 and y * = ∞ we have E(W 1 ) = 0.81, then we need to decrease c to c = 4, then we have for y * = 7, E(W 1 ) = 0.7801 and E(W 2 ) = 43.15 and for y * = 8, E(W 1 ) = 0.7777 and E(W 2 ) = 44.57. Hence the best solution to reach w * 1 = 0.78 is to have y * = 0 for c < 4, y * = 8 when c = 4 and y * = ∞ when c > 4. With this method, the best couple would

  Impact of p (µi = µ0 = 0.2, W * 0 = W * i = 0.2, ∑ 4 i=0 si = 44, i = 1, ..., 4, p ′ = 20%, U = V = 1) Impact of p ′ (λ0 = 2, λi = 1.5, W * 0

Figure A. 1 :

 1 Figure A.1: Performance measures of SP and chaining

  (a)-2.10(d) of Chapter 2, respectively. The experiments of Tables A.10-A.13 are associated to Figures 2.11(a)-2.11(d) of Chapter 2, respectively.

  the waiting time of our new call can be represented by Figure B.1. Each branch in B.1 is a possible scenario with a given summation of exponential random variables (with rates inside the circles). The subscript n or k in Figure B.1 indicates that the corresponding random variable is summed n or k times, respectively. A scenario with n exponential stages with rate µ 1 , n exponential stages with rate µ 2 , n exponential stages with rate µ 3 and k exponential stages of rate µ 0 , occurs with probability

. 7 )

 7 We have M = P -1 (N + D)P . Then e xM = P -1 e x(N +D) P and because N and D are commuting (Dunford decomposition) we have e x(N +D) = e x(N +D) = e xD × e xN . Because xD is a

  B.9)for i = 1, 2, 3. Since N n i = O for i = 1, 2, 3, we have N n = O. For k ≥ n, N k = O. We can write e xN = +∞ ∑ k=0 (xN ) k k! .The difficulty is to evaluate an infinite sum but because for k ≥ n we haveN k = Othis sum is in fact a finite one and e xN = We still need to evaluate the coefficients of thematrix N k for k ∈ 0, 1, 2, • • • , n -1. We have line i and column i + k for i ∈ 1, 2, • • • , n -1 of N k 1 , N k 2 and N k 3 are µ k 1 , µ k 2 andµ k 3 respectively. The other coefficient of N k are null. Then

F

  (x) = 1αe xM 1. (B.14)In this formula we note that we need to sum in fact the coefficients of the first line of e xM . In the basis B we haveαe xM 1 = 3n ∑ k=1 < f (e k )|e 1 >. We have f (e k ) = k ∑ i=1 e -xµ 1 (xµ 1 ) k-i (k-i)! e i for 1 ≤ k ≤ n. Then for 1 ≤ k ≤ n, < f (e k )|e 1 >= e -xµ 1 (xµ 1 ) k-1 (k-1)! . For 2 ≤ k ≤ n, f (e k+n ) = k ∑ i=2 e -xµ 2 (xµ 2 ) k-i (k-i)! e n+i + e -xµ 2 (xµ 2 ) k-1 (k-1)! v 2 . Because v 2 = k , then < v 2 |e 1 >= 1 and for 2 ≤ k ≤ n, < f (e k+n )|e 1 >= e -xµ 2 (xµ 2 ) k-1 (k-1)! . For 2 ≤ k ≤ n, f (e k+2n ) = k ∑ i=2 e -xµ 3 (xµ 3 ) k-i (k-i)! e 2n+i +e -xµ 3 (xµ 3 ) k-1 (k-1)! v 3 .Solution of a cubic equation: Consider the general cubic equation written asx 3 + ax 2 + bx + c =0, (B.22) with x as variable and the real parameters a, b and c. Changing the variable x into z = x -a/3 leads to z 3 + (b -a 2

x 4 +

 4 ax 3 + bx 2 + cx + d = 0, with x as variable and the real parameters a, b, c and d. Changing the variable x into z = x -a/4 leads toz 4 + pz 2 + qz + r = 0, with p = b -3a 2 8 , q = c -ab 2 + a 3 8 and r = d -ac 4 + a 2 b 16 -3a 4 256 .The solving method is based on computing the three real numbers α, β and δ such thatz 4 + pz 2 + qz + r = (z 2 + δ) 2 -(αz + β) 2 . (B.24)Equation (B.24) holds only if(z 2 + δ) 2 -(z 4 + pz 2 + qz + r) = (2δ -p)z 2 -qz + (δ 2 -r), (B.25)is a square. In other words, the discriminant of the right hand side of Equation (B.25) is equals to 0. Thenq 2 -4(δ 2 -r)(2δ -p) = -8δ 3 + 4pδ 2 + 8cδ + q 2 -4pr = 0. (B.26)Equation (B.26), with δ as variable is a cubic equation that can be solved using the method described in the first part of this section. Having in hand δ, we easily deduce α and β. It finally remains to note that the solutions in z of (z 2 + δ) 2 -(αz + β) 2 = 0 are -√

  Figure C.1: Email Throughput (s = 10, µ 0 = µ = 0.2, λ = 1.4)

  , for 0 ≤ u ≤ s -2. Recall that we are trying to prove Inequality (C.2). It then remains to prove that

Table 2

 2

					.1: Fixed point approximation, µ = 0.2
						Wi		PD
		λ	λ0	s	(λ+λ 0 )/n sµ	Simulation Approximation Simulation Approximation
		0.35	0.35	5	70%	0.581	0.581	37.78%	37.78%
		0.475 0.475	5	95%	1.672	1.672	87.78%	87.78%
	n = 1	1.4	1.4	20	70%	0.035	0.035	9.36%	9.36%
		1.9	1.9	20	95%	0.359	0.359	75.54%	75.54%
		3.8	0	20	95%	3.777	3.777	75.54%	75.54%
		0.35	0.7	5	70%	0.436	0.435	29.74%	28.28%
		0.475	0.95	5	95%	1.623	1.623	87.51%	85.23%
	n = 2	1.4	2.8	20	70%	0.003	0.0027	0.74%	0.72%
		1.9	3.8	20	95%	0.320	0.290	61.03%	60.98%
		3.8	0	20	95%	3.777	3.777	75.54%	75.54%
		0.35	1.75	5	70%	0.290	0.288	19.06%	18.76%
		0.475 2.375	5	95%	1.556	1.550	81.84%	81.38%
	n = 5	1.4	7	20	70%	0.001	0.001	0.28%	0.27%
		1.9	9.5	20	95%	0.205	0.204	42.99%	42.97%
		3.8	0	20	95%	3.777	3.777	75.54%	75.54%
		0.35	3.5	5	70%	0.259	0.252	17.21%	16.39%
		0.475	4.75	5	95%	1.516	1.504	79.07%	78.96%
	n = 10	1.4	14	20	70%	0.0001	0.0001	0.23%	0.23%
		1.9	19	20	95%	0.167	0.167	35.23%	35.21%
		3.8	0	20	95%	3.777	3.777	75.54%	75.54%

  1). In particular, we are interested to know, for the different ranges of p, which one of the

	models would be preferred to the other. We choose call center examples with n = 4, i.e., 5 agent
	teams and 5 skills including skill 0. The results are shown in Table 2.2 and Figures 2.5(a) and
	2.5(b).
	Table 2.2: Impact of p (µ i

Table 2

 2

	.2 reveals that chaining performs well for small values

Table 2 .3 and Figure

 2

Table 2 .5: Impact of variability in service times (µ i

 2

Table 2

 2

	.6, and Figures 2.8(a) and 2.8(b).

Table 2 .

 2 8: Probability of abandonment (µ i

Table 2 .

 2 9: Impact of the Call Center Size (µ i

Table 2

 2

.14: Impact of p and U (λ

Table 2 .

 2 15: Impact of p ′ and U (λ i

  20% )

					Chaining			SP	Crossing value
	U	V	t=0% t=5% t=10% t=25% t=50%		(Chaining = SP)
		1/3	50	51.45	52.9	57.25	64.5	55	t=17.24%
		1/2	50	51.55	53.1	57.75	65.5	56	t=19.35%
	1	1	49	50.7	52.4	57.5	66	56	t=20.58%
		2	50	51.55	53.1	57.75	65.5	56	t=19.35%
		3	50	51.45	52.9	57.25	64.5	55	t=17.24%
		1/3	26	26.8	27.6	30	34	34	t=50.00%
		1/2	31	32	33	36	41	37	t=30.00%
	2	1	49	50.6	52.2	57	65	55	t=18.75%
		2	71	72.6	74.2	79	87	76	t=15.63%
		3	81	82.3	83.6	87.5	94	84	t=11.54%
		1/3	20	20.55	21.1	22.75	25.5	26	t=54.55%
		1/2	24	24.65	25.3	27.25	30.5	31	t=53.85%
	3	1	50	51.4	52.8	57	64	53	t=10.71%
		2	79	80.65	82.3	87.25	95.5	82	t=9.09%
		3	93	94.35	95.7	99.75	106.5	95	t=7.41%

Table 2

 2

						Scenarios						Chaining	SP	Crossing value
		λ1	λ2	λ3	λ4	λ0	µ1	µ2	µ3	µ4	µ0	t = 0% t = 10% t = 20%	(Chaining=SP)
	Sc 1	1	2	3	4	5	0.05 0.1 0.2 0.5	1	78	83.5	89	87	16.36%
	Sc 2	2	3	4	5	1	0.05 0.1 0.2 0.5	1	115	122.6	130.2	127	15.79%
	Sc 3	1	2	3	4	5	1	0.5 0.2 0.1 0.05	179	184.9	190.8	184	8.47%
	Sc 4	2	3	4	5	1	1	0.5 0.2 0.1 0.05	111	116.9	122.8	119	13.56%
	Sc 5 0.1 0.2 0.3 0.4 0.5 0.05 0.1 0.2 0.5	1	14	15	16	18	40.00%
	Sc 6 0.2 0.3 0.4 0.5 0.1 0.05 0.1 0.2 0.5	1	18	19.4	20.8	24	42.86%
	Sc 7 0.1 0.2 0.3 0.4 0.5	1	0.5 0.2 0.1 0.05	26	26.7	27.4	30	57.14%
	Sc 8 0.2 0.3 0.4 0.5 0.1	1	0.5 0.2 0.1 0.05	18	18.7	19.4	21	42.86%

.17: Impact of p, p ′ , U and V (W * i = 0.2 for i = 0, • • • 4) ŚĂŝŶŝŶŐ (a) Impact of p and p ′ (µi = µj and λi = λj for

Table 3 .

 3 2: Expressions of T , E(W ) and P D for Models 1,...,4

	Model 1	Model 2

Table 3 .

 3 4: Comparison between approximation and simulationInterval of the call arrival rate λ s = 10, µ 0

Table 4 .

 4 1: Optimal couples of thresholds (s = 10, τ = 30 seconds, α = 80%)

	λ1	λ2	µ	µ0	I1	I2	u * 1	u * 2	(u1, u2) T1	T2	T
	1	1	0.2 0.2 50% 50%	8	8	(8,8)	84.04%	84.04%	84.04%	0.758 0.758 0.758
	1	1.3 0.2 0.2 50% 50%	8	6	(8,7)	84.04%	77.99%	80.62%	0.758 0.401 0.580
	0.5 1.5 0.2 0.2 50% 50%	9	-	(8,4)	96.81%	74.79%	80.30%	1.169 0.055 0.611
	1	1.3 0.2 0.2 67% 33%	8	6	(8,7)	84.04%	77.99%	81.66%	0.758 0.401 0.639
	1	1.3 0.2 0.2 80% 20%	8	6	(8,8)	84.04%	69.15%	80.39%	0.758 0.552 0.711
	0.5 1.5 0.2 0.2 90% 10%	9	-	(9,7)	88.19%	63.94%	82.13%	1.350 0.277 1.243
	1	1.5 0.2 0.2 50% 50%	8	-	(7,5)	90.92%	72.93%	80.13%	0.604 0.111 0.357
	1	1.5 0.2	1	50% 50% 10 -	(10,7)	89.34%	74.94%	80.70%	5.191 0.961 3.076
	1	1.5 0.2	1	80% 20% 10 -	(10,10)	89.34%	67.56%	83.40%	5.191 2.908 4.734
	1.3 1.4 0.2	1	50% 50%	9	8	(9,9)	83.51%	77.09%	80.18%	2.863 2.440 2.652
	1.3 1.4 0.2	1	80% 20%	9	8	(9,10)	83.51%	68.19%	80.26%	2.863 3.621 3.014
	1.3 1.4	1	0.2 50% 50%	9	9	(9,9)	88.63%	87.77%	88.18%	1.616 1.598 1.601
	1.3 1.4	1	0.2 80% 20%	9	9	(9,10)	88.63%	60.45%	82.10%	1.616 1.794 1.742

* P (W1 < τ ) P (W2 < τ ) P (W < τ )

Table 4

 4

				.3: Comparison between ATP and constant step methods
		h	T	SL	r	U		h	T	SL	r	U
		0.1	1.17 80.6% 0.0046 0.71		0.1	1.53 72.15% 0.0782 -6.3
	Sc 1	0.2 0.5	1.12 80.5% 0.0036 0.77 1.04 80.1% 0.0032 0.72	Sc 2	0.2 0.5	1.38 78.7% 0.0201 -0.63 1.23 81.4% 0.0063 0.60
		1	0.98 80.0% 0.0035 0.63		1	1.19 80.7% 0.0062 0.57
		ATP 1.09 80.7% 0.0027 0.82		ATP 1.12 85.6% 0.0008 1.04
		0.1	1,85 80,3% 0,0023 1.62		0.1	3.20 78.9% 0.0314 0.06
	Sc 3	0.2 0.5	1,80 80,3% 0,0017 1.63 1,68 80,3% 0,0013 1.55	Sc 4	0.2 0.5	3.07 79.5% 0.0277 0.30 3.05 79.2% 0.0278 0.27
		1	1,57 80,2% 0,0014 1.43		1	3.14 79.2% 0.0281 0.33
		ATP 1,72 81,0% 0,0003 1.68		ATP 2.95 78.9% 0.0264 0.31
		0.1	1.19 79.9% 0.0067 0.52		0.1	1.05 83.2% 0.0014 0.91
	Sc 5	0.2 0.5	1.13 80.1% 0.0037 0.76 1.08 80.0% 0.0033 0.75	Sc 6	0.2 0.5	1.04 81.7% 0.0021 0.83 1.04 80.8% 0.0025 0.79
		1	1.01 79.9% 0.0033 0.68		1	1.04 80.2% 0.0032 0.72
		ATP 1.12 80.4% 0.0018 0.93		ATP 1.09 82.1% 0.0007 1.02
		0.1	1.38 81.6% 0.0010 1.28		0.1	3,04 78,8% 0,0246 0,59
	Sc 7	0.2 0.5	1.37 81.2% 0.0015 1.21 1.27 80.4% 0.0017 1.10	Sc 8	0.2 0.5	2,84 79,5% 0,0201 0,83 2,72 79,3% 0,0178 0,94
		1	1.24 80.3% 0.0011 1.14		1	2,76 78,9% 0,0188 0,88
		ATP 1.38 81.2% 0.0005 1.33		ATP 2,83 79,7% 0,0172 1,11

Table 4

 4

						.5: Comparison of the methods
		h	T	SL	r	U		h	T	SL	r	U
		M 1 1.01 79.9% 0.0038 0.62		M 1 1.25 80.45% 0.0080 0.44
		M 2 0.99 79.9% 0.0037 0.62		M 2 1.13 81.2% 0.0068 0.45
		M 3 1.45 72.6% 0.0743 -5.98		M 3 1.59 68.0 % 0.0639 -4.79
		M 4a 1.01 80.9% 0.0041 0.59		M 4a 1.24 80.7% 0.0090 0.34
	Sc 1	M 4b 1.06 80.0% 0.0029 0.76 0.1 1.17 80.6% 0.0046 0.71	Sc 2	M 4b 1.20 80.2% 0.0099 0.22 0.1 1.53 72.15% 0.0782 -6.3
		0.2	1.12 80.5% 0.0036 0.77		0.2	1.38 78.7% 0.0201 -0.63
		0.5	1.04 80.1% 0.0032 0.72		0.5	1.23 81.4% 0.0063 0.60
		1	0.98 80.0% 0.0035 0.63		1	1.19 80.7% 0.0062 0.57
		ATP 1.09 80.7% 0.0027 0.82		ATP 1.12 85.6% 0.0008 1.04

Table 5 .

 5 1: E(W 1 ), E(W 2 ), E(W ) and π as a function of c with λ = 1, s = 6, k = 3 and q = 80%

	c	1	2			3	4	5	6
	E(W 1 ) E(W 2 ) E(W ) π	0.65 ∞ ∞ 7.51% 7.76% 8.18% 8.62% 9.20% 10.11% 0.67 0.7 0.75 0.81 0.9 1698 68.33 32.85 20.85 ∞ 142.83 6.63 3.75 2.94 ∞
	probabilities and the relation πλE(W 2 ) =	∞ ∑	∞ ∑	yp x,y . Since the sums of this expression are infinite
				y=0	x=0
	we can find a significant numerical value by space state truncation.
	Figures 5.2, 5.3, 5.4 and Table

Table 5 .

 5 2: Impact of c and K (λ = 1, µ = 0.2, s = 6, k = 3, q = 80%)

	c = 3	K E(W 1 ) E(W 2 ) E(W )	0 2.94 2.94	1 2.69 51.71 57.30 62.84 68.57 73.98 101.57 541.8 2 3 4 5 10 100 2.49 2.33 2.19 2.06 1.68 0.84 3.20 3.51 3.85 4.20 4.57 6.64 44.28 142.83 ∞ 0.7 1698
		π	0.00% 1.04% 1.85% 2.51% 3.04% 3.49% 4.96% 8.03% 8.18%
	c = 4	E(W 1 ) E(W 2 ) E(W )	2.94 2.94	2.66 35.44 37.11 38.90 40.57 42.33 2.39 2.19 2.02 1.88 3.08 3.2 3.36 3.52 3.7	1.36 49.13 4.45	0.76 68.28 6.61	0.75 68.33 6.63
		π	0.00% 1.28% 3.32% 3.18% 3.90% 4.50% 6.48% 8.56% 8.62%
	c = 5	E(W 1 ) E(W 2 ) E(W )	2.94 2.94	2.60 26.36 27.15 27.48 2.36 2.11 2.96 3.03 3.06	1.94 28 3.13	1.78 28.62 3.2	1.24 30.03 3.39	0.82 32.58 3.72	0.81 32.85 3.75
		π	0.00% 1.50% 2.72% 3.74% 4.59% 5.29% 7.46% 9.17% 9.20%

Table 5 .

 5 3: Impact of c and y * (λ = 1, µ = 0.2, s = 6, k = 3, q = 80%) = 16 which leads to E(W 1 ) = 0.7799 and E(W 2 ) = 50.63.

	c = 3	y * E(W 1 ) E(W 2 ) E(W )	0 0.9 20.85 2.94	1 0.88 24.38 27.81 2 0.87 3.22 3.49	3 0.85 32 3.85	4 0.84 36.23 40.72 65.01 515.97 1698 5 10 100 ∞ 0.83 0.8 0.74 0.7 4.19 4.6 6.64 44.5 142.83
		π	10.11% 9.93% 9.75% 9.64% 9.48% 9.44% 9.10% 8.49% 8.18%
	c = 4	E(W 1 ) E(W 2 ) E(W )	0.9 20.85 2.94	0.88 22.78 24.83 27.09 29.05 31.13 41.06 0.86 0.85 0.84 0.83 0.8 3.04 3.19 3.37 3.51 3.65 4.45	0.76 68.38 6.64	0.75 68.33 6.63
		π	10.11% 9.87% 9.71% 9.59% 9.46% 9.33% 9.06% 8.69% 8.62%
	c = 5	E(W 1 ) E(W 2 ) E(W )	0.9 20.85 2.94	0.89 21.88 22.77 23.54 24.42 24.99 28.31 0.87 0.86 0.86 0.85 0.83 2.98 3.03 3.07 3.13 3.16 3.40	0.81 32.57 3.72	0.81 32.85 3.75
		π	10.11% 9.98% 9.84% 9.73% 9.67% 9.57% 9.35% 9.18% 9.20%
	be c = 4 and y

*

Table A .

 A 1: Decreasing greedy in team 0 for single pooling (n = 2, µ i

  Table A.2 provides the simulation results of this algorithm. Table A.2: Increasing greedy for single pooling (n = 2, µ i

Table A .

 A 3: Increasing greedy with no Agents in team 0 for single pooling (n = 2, µ i = µ 0 = 0.2,

	W * i = W * 0 = 0.2, i = 1, 2)
	λ1	λ2	λ0	FF	FD	Single pooling s1 s2 s0	Total s0 + s1 + s2
	1 0.2 0.5 0.5 0.2 1 1 1 1 3 2 1 2 1 3 0.5 0.2 0.1 0.1 0.2 0.5 10 5 15 10 15 5 10 10 10	13 13 20 36 36 8 8 151 151 151	19 19 27 44 44 13 13 171 171 171	9 6 11 10 6 7 22 16 21 16 6 4 4 5 93 68 70 93 81 81	0 0 0 0 0 0 0 0 0 0	15 15 21 38 37 10 9 161 163 162

Table A

 A

.4: Efficiency of algorithm 2 (n = 2, µ i

Table A

 A

	.13: Impact of U (λ 0 = 4, W * 0 = W * i = 0.2, µ 0 = 0.2, i = 1, ..., 4, p = 50%, p ′ = 20%,
	V = 1)
					Chaining			SP	Crossing value
		U t=0% t=5% t=10% t=25% t=50%		(Chaining = SP)
		1	47	47.95	48.9	51.75	56.5	52	t=26.32%
	γ0 = 0.1	2	48	48.8	49.6	52	56	53	t=31.25%
	γi = 0	3	48	48.85	49.7	52.25	56.5	52	t=23.53%
		5	48	48.9	49.8	52.5	57	52	t=22.22%
		1	46	46.9	47.8	50.5	55	51	t=27.78%
	γi = γ0 = 0.1	2	48	48.9	49.8	52.5	57	51	t=23.53%
		3	50	50.85	51.7	54.25	58.5	51	t=5.88%
		5	51	51.75	52.5	54.75	58.5	51	t=0.00%
		1	48	48.85	49.7	52.25	56.5	52	t=23.53%
	γ0 = 0	2	48	48.7	49.4	51.5	55	49	t=7.14%
	γi = 0.1	3	49	49.9	50.8	53.5	58	50	t=5.56%
		5	50	51	52	55	60	51	t=5.00%

Table B .

 B 1: Performance Measures in function of the level of simultaneity (µ = 0.2)

		k	1	2	3	5	10	+∞
	(λ = 0.1, s = 1)	W q W s W q + W s	5 5 10	2.5 7.5 10	1.25 8.75 10	0.31 9.69 10	0.01 9.99 10	0 10 10
	(λ = 0.19, s = 1)	W q W s W q + W s	95 5 100	90.25 85.74 77.38 59.87 9.75 14.26 22.62 40.13 100 0 100 100 100 100 100
	(λ = 0.5, s = 6)	0.068 0.000 0.000 0.000 0.000 0.000 5.000 5.311 5.315 5.317 5.318 5.318 W q + W s 5.068 5.311 5.315 5.317 5.318 5.318 W q W s
	(λ = 1, s = 6)	2.94 5.00 W q + W s 7.94 11.34 13.89 16.74 17.05 17.34 1.95 1.08 0.21 0.02 0.00 W q W s 9.39 12.81 16.53 17.03 17.34
	(λ = 1.15, s = 6)	17.75 16.38 15.03 11.98 6.87 5.00 9.90 15.09 24.51 43.66 64.71 0.00 W q + W s 22.75 26.28 30.12 36.49 50.53 64.71 W q W s

  xµ i e -xµ i xµ i e -xµ i x 2 2 µ 2

	0 . . . 0	e -xµ i . . . 0	i e -xµ i xµ i . . . . . .	. . . e -xµ i x n-2 (n-2)! µ n-2 i . . . e -xµ i x n-3 (n-3)! µ n-3 i . . . . . . e -xµ i e -xµ i xµ i	e -xµ i x n-1 (n-1)! µ n-1 i e -xµ i x n-2 (n-2)! µ n-2 i . . . e -xµ i x 2 i 2 µ 2	          	. (B.13)
							
	0	0	. . .	0	e -xµ i	e -xµ i xµ i	 
							
	0	0	. . .	0	0	e -xµ i

Table C .

 C 1: Comparison under steady-states assumption (θ=15min)

	λ	µ	µ 0	s	γ	Optimal c Optimal T ATP T Difference
	0.005 0.2	0.2	1	0.22	0.80	0.04	0.04	0.00%
	0.02	0.2	0.2	1	0.23	0.13	0.02	0.02	0.00%
	0.05	0.2	0.2	1	0.35	0.06	0.01	0.01	0.00%
	0.02	0.2	1	1	1.68	0.20	0.18	0.18	0.00%
	0.02	1	0.2	1	0.22	0.21	0.04	0.04	0.00%
	0.1	0.2	0.2	5	0.32	4.34	0.76	0.75	1.33%
	0.3	0.2	0.2	5	1.25	3.93	0.55	0.50	10.00%
	0.5	0.2	0.2	5	1.65	2.75	0.23	0.21	9.52%
	0.5	0.2	1	5	6.21	4.29	2.03	1.88	7.98%
	0.5	1	0.2	5	0.39	4.17	0.78	0.73	6.85%
	0.1	0.2	0.2	10	0.52	9.51	1.80	1.72	4.65%
	1	0.2	0.2	10	0.99	8.72	0.80	0.77	3.90%
	1.45	0.2	0.2	10	5.28	1.09	0.01	0.01	0.00%
	8.2	1	1	10	17.89	5.94	0.38	0.36	5.56%
	0.1	0.05 0.05 10	0.06	9.94	0.36	0.33	9.09%
	1.3	0.2	1	10	18.67	9.30	3.46	3.41	1.46%
	1	1	0.2	10	0.66	9.37	1.78	1.77	0.05%
	2	0.2	0.2	28	2.20	27.37	3.29	2.99	10.03%
	4	0.2	0.2	28	1.99	25.49	1.39	1.37	1.46%
	1	0.1	0.1	28	0.24	27.19	1.74	1.71	1.75%
	4	0.27 0.15 28	0.93	26.63	1.89	1.89	0.00%
	4	0.17	1	28 196.31	23.21	2.00	1.79	11.73%
	17.5	0.2	0.2 100	3.83	97.77	2.37	2.16	9.72%
	18	0.2	0.2 100	7.73	93.91	1.65	1.58	4.43%
	18.5	0.2	0.2 100	2.30	85.61	0.58	0.54	7.41%
	95	1	0.2 100 11.13	88.28	0.82	0.73	12.89%
	for the case of different service rates.

  Le contexte est celui de centres d'appels avec des paramètres asymétriques : charge de travail non équilibrée, différents temps de services, prédominance d'une catégorie de clients, taux d'abandons variables et coûts élevés de la multi-compétence. Les architectures les plus connues avec flexibilité limitée comme chaining ne résistent pas à de telles asymétries. Nous proposons une nouvelle architecture, appelée Single Pooling avec seulement deux compétences par agent et nous démontrons son efficience dans diverses situations d'asymétrie.Dans la seconde partie, nous nous intéressons aux problèmes de routage dans les centres d'appels multi-canaux. Dans le premier projet, nous considérons un centre avec des appels arrivant au fil du temps et des emails présents en nombre illimité. Le service des appels se fait en trois étapes dont la seconde est une pause pour l'agent. Nous cherchons à optimiser le routage des emails. L'objectif est de maximiser le débit d'emails traités sous contrainte de temps d'attente pour les appels. De nombreuses recommandations sont proposées au manager. En particulier, nous démontrons que pour obtenir un routage optimal il est nécessaire de fixer à une valeur extrême au moins l'un des deux paramètres définissant le routage des emails.Dans le second projet, nous étudions une politique de seuil de réservation d'agents pour les appels en réception. Nous considérons un cas général de modèle non stationnaire où le processus d'arrivée des appels est Poisson non homogène. Le problème d'optimisation est la maximisation du débit de tâches en émission sous contrainte de qualité de service sur les appels en réception. Nous proposons une méthode efficiente et facile à implémenter de changement adaptatif de seuil. Cette politique est évaluée en comparaison avec les performances optimales trouvées dans le cas particulier de taux d'arrivée constant, et en comparaison avec d'autres méthodes intuitives de changement adaptatif de seuil dans le cas général non stationnaire.Dans le troisième projet, nous considérons un modèle de centre avec option de rappels. Cette option permet de transformer un appel en réception en un appel en attente d'émission. Le problème d'optimisation est la minimisation du temps d'attente des appels en émission sous contrainte de qualité de service pour les appels en réception. Nous proposons une politique de routage à deux seuils, un sur la réservation d'agents pour les appels en réception et un sur le nombre d'appels en attente d'émission. Nous déterminons une courbe optimale entre ces deux seuils.Mots clefs centres d'appels, modèles stochastiques, files d'attente, chaine de Markov, simulation, politique d'ordonnancement, qualité de service, analyse transitoire, routage par compétences, centres d'appels multicanaux

  call centers, stochastic models, queuing systems, Markov chains, simulation, scheduling policies, quality of service, transient analysis, skill-based routing, multi-channel call centersRésuméLes centres d'appels connaissent un grand succès depuis leur introduction dans les entreprises de service. Ils sont le principal point de contact avec les clients, et une composante essentielle de la majorité des entreprises. L'émergence des centres d'appels à grande échelle a suscité de nombreuses problématiques de management. Dans cette thèse, nous considérons des problématiques de management orientées sur les centres multi-canaux et multi-compétences. L'objectif de notre travail est de trouver des résultats qualitatifs et quantitatifs utiles pour le management.Dans la première partie, nous considérons les architectures de centres multi-compétences à flexibilité limitée. Le contexte est celui de centres d'appels avec des paramètres asymétriques : charge de travail non équilibrée, différents temps de services, prédominance d'une catégorie de clients, taux d'abandons variables et coûts élevés de la multi-compétence. Les architectures les plus connues avec flexibilité limitée comme chaining ne résistent pas à de telles asymétries. Nous proposons une nouvelle architecture, appelée Single Pooling avec seulement deux compétences par agent et nous démontrons son efficience dans diverses situations d'asymétrie.Dans la seconde partie, nous nous intéressons aux problèmes de routage dans les centres d'appels multi-canaux. Dans le premier projet, nous considérons un centre avec des appels arrivant au fil du temps et des emails présents en nombre illimité. Le service des appels se fait en trois étapes dont la seconde est une pause pour l'agent. Nous cherchons à optimiser le routage des emails. L'objectif est de maximiser le débit d'emails traités sous contrainte de temps d'attente pour les appels. De nombreuses recommandations sont proposées au manager. En particulier, nous démontrons que pour obtenir un routage optimal il est nécessaire de fixer à une valeur extrême au moins l'un des deux paramètres définissant le routage des emails.Dans le second projet, nous étudions une politique de seuil de réservation d'agents pour les appels en réception. Nous considérons un cas général de modèle non stationnaire où le processus d'arrivée des appels est Poisson non homogène. Le problème d'optimisation est la maximisation du débit de tâches en émission sous contrainte de qualité de service sur les appels en réception. Nous proposons une méthode efficiente et facile à implémenter de changement adaptatif de seuil. Cette politique est évaluée en comparaison avec les performances optimales trouvées dans le cas particulier de taux d'arrivée constant, et en comparaison avec d'autres méthodes intuitives de changement adaptatif de seuil dans le cas général non stationnaire.Dans le troisième projet, nous considérons un modèle de centre avec option de rappels. Cette option permet de transformer un appel en réception en un appel en attente d'émission. Le problème d'optimisation est la minimisation du temps d'attente des appels en émission sous contrainte de qualité de service pour les appels en réception. Nous proposons une politique de routage à deux seuils, un sur la réservation d'agents pour les appels en réception et un sur le nombre d'appels en attente d'émission. Nous déterminons une courbe optimale entre ces deux seuils.Mots clefs centres d'appels, modèles stochastiques, files d'attente, chaine de Markov, simulation, politique d'ordonnancement, qualité de service, analyse transitoire, routage par compétences, centres d'appels multicanaux

(3.11) to easily obtain those of Equation (3.12).

Model 3: The Markov chain associated to Model 3 is given in Figure 3.4. For this model, we have

C, n+1 From the Markov chain of Model 3, we may write c n (λ + µ 3 ) = λc n-1 + µ 0 b ′ n for n ≥ 1, and

for n ≥ 2. We also have from the Markov chain b ′ n (λ + µ 0 ) = λb ′ n-1 + µ 2 b n , for n ≥ 1. Therefore, b n = ρ 2 (1 + ρ 0 )(1 + ρ 3 )u n-1ρ 2 (1 + ρ 0 )ρ 3 u n-2 -ρ 2 ρ 0 (1 + ρ 3 )u n-2 + ρ 2 ρ 0 ρ 3 u n-3 , for n ≥ 3. From a state (B, n) for n ≥ 1, we may write b n (λ + µ 2 ) = λb n-1 + µ 1 a n , which leads to Chapter 5

Call Centers with a Callback Option

We consider a call center model with a callback option, which allows to transform an inbound call into an outbound one. The optimization problem consists on minimizing the expected waiting time of the outbound calls while respecting a service level constraint on the inbound ones. We propose a routing policy with two thresholds, one on the reservation of agents for inbound calls, and another on the number of waiting outbound calls. The purpose of this study is to determine a curve relating the two thresholds. The paper version of this chapter is the ongoing paper [START_REF] Legros | Call centers with a call back option[END_REF].

Appendix A

Appendix of Chapter 2 A.1 Optimization Heuristic

A coherent comparison between chaining and single pooling requires first the optimization of their total cost. In this section, we develop a greedy heuristic for the optimization step of the two models, and prove their efficiency. The heuristic is a simulation based optimization method. Recall that for each model, we optimize the total staffing cost under the constraints W i ≤ W * i , for i = 0, 1, ..., n.

In the numerical examples below, we consider Markovian assumptions for inter-arrival and service times. The analysis can be applied in a similar way to any other assumption.

The question addressed here is how can we compute the optimal number of agents in each team?

In some particular cases the answer is simple. For example when the arrival and service rates are identical for all skills and when all skills have the same costs, we would create teams with the same number of agents. A more difficult situation is in the case of asymmetric arrival or service rates.

For the simulation based optimization considered here, some information about the simulation process are as follows. We use C++ program codes. For a given simulation with a given set of parameters, we consider a single replication that we run for a sufficiently long time. The lengths of the confidence intervals for the different performance measures derived by simulation are in the order of 10 -4 . To obtain such confidence intervals, we simply gradually increase the replication a high number of small steps).

The method is as follows. In each team i (i = 0, 1, ..., n), we start with the worst (overestimated) staffing level s i computed from an FD model. In order to take into account the chaining configuration, i.e., the fact that customers type i -1 can be routed to team i and customers type i can be routed to team i+1, we adjust the initial staffing levels from s i to s ′ i for team i (i = 0, 1, ..., n).

We use the method suggested by [START_REF] Wallace | A Staffing Algorithm for Call Centers with Skill-Based Routing[END_REF]. The corrected staffing level s ′ i is given by

for i = 0, 1, ..., n, where R i,j = s i s j s-s i . This number is that of agents of team i who could go to team j, i, j = 0, 1, ..., n and i ̸ = j. Using Equation (A.1), s ′ i may not be an integer. We then round it to the nearest integer above.

A.2 Metric Comparison

In Table A.5 and Figures A.1(a)-A.1(d), we compare between the expected waiting times of SP and chaining for a given total staffing level. We optimize the staffing of the various teams in the two models for the case t = 0, i.e., no incremental cost for regular skills. The results show the same qualitative conclusions as those in Chapter 2. In this section we focus on the possibility to work simultaneously on different tasks (usually back office tasks). This way of working is not common in call centers as we know that a human resource is limited in the simultaneous treatment. However, we wonder if the possibility of simultaneity goes with an increasing of the performance measures. We consider a model with no increasing in productivity when an agent treats tasks simultaneously. A simple classical modeling of this situation is to suppose that an agent can treat at most k tasks simultaneously. The tasks arrive one by one according to a Poisson arrival process of rate λ. We also suppose that the service time distribution follows an exponential distribution of rate µ/n for 1 ≤ n ≤ k when n tasks are done simultaneously by an agent. When an agent treats more than one task at the same time we suppose that the service rate changes automatically in function of the number of tasks in service when one task finishes or starts the service. We denote by s the number of agents in the call center. When a task arrives, if all agents are busy with k tasks, the new task waits in the queue. If not this new task is routed to the less busy agent. If two or more agents are the least busy then the new task is routed to one of these agents with the same probability. When an agent finishes a task she starts a new one if there is one in the queue, if not she finishes her remaining tasks in service if any.

) .

We may also write

) .

for n ≥ 1. Finally, note that all the other quantities can be computed in the same way.

B.3 Method for Deriving P (W < t) in a Three-Stages Hypoexponential Distribution Service

In this section we propose a method to derive P (W < t) for an arriving customer in a system with a single agent and distributed according to a three-stages distribution with rates µ 1 , µ 2 and µ 3 . Recall that the queue is infinite, we do note consider any abandonment or retrial and customers are served according to a FCFS rule. The method is based on the Dunford decomposition for the sub-generator matrix. This method conducts to close form formulas which are easier to compute than the general formula of a distribution (when there are distinct sums of exponential distributions with rates λ 1 ,λ 2 ,• • • ,λ a and a number of terms in each sum equals to r 1 , r 2 , • • • , r a respectively. The cumulative distribution function for t ≥ 0 is given by

). Yet this method can hardly be generalized to other phase-type distribution which explain why we do not present it as a main contribution.

The distribution is a phase-type distribution. If we denote by M the sub-generator matrix of this distribution. The cumulative distribution function is given by F (x) = 1αe xM 1 where 1 is a column vector of ones of the size k (k stages), e A is the matrix exponential of A and α = (1, 0, . . . , 0). The difficulty is to evaluate e xM . We propose a method based on the Dunford decomposition. We consider an arriving customer when n (n ≥ 1) customers are already in the system, one in service in stage 1 and n -1 is the queue. Thus the distribution is composed with n exponential distributions of rate µ 1 , n exponential distributions of rate µ 2 and n exponential distributions of rate µ 3 . Note that when the customer in service is in stage 2 or 3 we find one stage less of rate µ 1 and/or µ 2 . The sub-generator matrix of this distribution is of dimension 3n×3n. We denote by m i,j its coefficient on line i and column j for i, j ∈ {1, 2

When µ 1 = µ 2 = µ 3 , we have the well known Erlang distribution. Thus we consider here

We still need to evaluate < f (e n+1 )|e 1 > and

This equality is equivalent to

) .

This implies that

Consequently when n customers are in the system and one customer is in the first stage of the service we have,

This finishes the demonstration of the method.

B.4 Reminder for the Cardan-Ferrari Method

In Sections 3.4.1 and 3.4.2, we use the Cardan-Ferrari method to solve cubic and quadric equations.

Below, we describe the principles of these methods. We also refer the reader to the textbook [START_REF] Gourdon | Les maths en tête : Algèbre[END_REF] for more details.

B.5 Proof of Proposition 2

We start by proving the third statement. It consists on solving the inequality λ 2 < λ 3 as a function of w * . This inequality is equivalent to

We then need to know the sign of w * 2 + w *

(

1 µ i µ j = 0 with w * as variable has two real solutions denoted by w * 1 and w * 2 . They are given by

Then, the numerator in Inequality (B.27) is strictly negative if and only if w * is strictly between w * 1 and w * 2 . It is easy to see that w * 1 < 0. We next prove that the second solution w * 2 is strictly higher than 1 µ 0 . We have

Multiplying the right hand side of Equation (B.28) by the following positive quantity 2

Since w * 2 > 1 µ 0 and defining w * as w * = w * 2 , Inequality (B.27) holds if and only if 1 µ 0 < w * < w * .

This finishes the proof of the third statement.

Let us now prove the first statement. It suffices to prove that Model 2 is not optimal for all λ > 0 if and only if λ 2 ≤ 0.

⇐ If λ 2 ≤ 0, Model 2 is obviously not optimal for all λ > 0.

⇒ Assume Model 2 is not optimal for all λ > 0. Then there are two possible reasons:

We prove next that the second inequality is impossible. The inequality

From the third statement of Proposition 3, we have:

The inequality R ≤ λ 4 is equivalent to

Having R ≤ λ 4 and λ 2 ≤ λ 3 implies

This previous inequality is impossible because

. The proof is as follows.

We have

Multiplying the last right hand side by the following positive quantity 2

This proves that the inequality R ≤ λ 4 ≤ λ 2 ≤ λ 3 is impossible and finishes the proof of the first statement.

The proof of the second statement of Proposition 3 is straightforward. The details are then omitted. It suffices to see that we choose Model 3 if λ 2 < λ < λ 3 . If not, Model 3 is chosen if R ≤ λ < λ 3 . This finishes the proof of the second statement and completes the proof of the proposition. ✷

B.6 Proof of Lemma 1

We want to solve the following inequality in ρ 0 :

B.7 Expression of the Probability q 2

This section is related to the proof of Theorem 1. Equation (B.30) gives the expression of the probability q 2 .

Appendix C

Appendix of Chapter 4 C.1 Results for the performance Measures

From the expressions of the performance measures found in Section 4.3.1, we can derive some monotonicity results. We conjecture that the email throughput is increasing and neither convex nor concave in u and that the call service level is decreasing and concave in u (for 0 ≤ u ≤ s). In the case of equal service rates we can prove Proposition 5.

Proposition 5 The following holds:

1. The email throughput T is strictly increasing and neither convex nor concave in u, for 0 ≤ u ≤ s. However the end of the email throughput, for s -2 ≤ u ≤ s and s ≥ 2, is concave in u.

2. The call service level P (W < τ ) is strictly decreasing and concave in u, for 0 ≤ u ≤ s and a < 1 or u + 1 ≤ a < s.

Proof. Let us prove the first statement. From Equation (4.2), we have

Abstract

Call centers have been introduced with great success by many service-oriented companies. They become the main point of contact with the customer, and an integral part of the majority of corporations. The large-scale emergence of call centers has created a fertile source of management issues. In this PhD thesis, we focus on various operations management issues of multi-skill and multichannel call centers. The objective of our work is to derive, both qualitative and quantitative, results for practical management.

In the first part, we focus on architectures with limited flexibility for multi-skill call centers. The context is that of call centers with asymmetric parameters: unbalanced workload, different service requirements, a predominant customer type, unbalanced abandonments and high costs of crosstraining. The most knowing architectures with limited flexibility such as chaining fail against such asymmetry. We propose a new architecture referred to as single pooling with only two skills per agent and we demonstrate its efficiency under various situations of asymmetry.

In the second part, we focus on routing problems in multi-channel call centers. In the first study, we consider a blended call center with calls arriving over time and an infinitely backlogged queue of emails. The call service is characterized by three successive stages where the second one is a break. We focus on the optimization of the email routing to agents. The objective is to maximize the throughput of emails subject to a constraint on the call waiting time. Various guidelines to call center managers are provided. In particular, we prove for the optimal routing that all the time at least one of the two email routing parameters has an extreme value.

In the second study, we examine a threshold policy on the reservation of agents for the inbound calls. We study a general non-stationary model where the call arrival follows a non-homogeneous Poisson process. The optimization problem consists on maximizing the throughput of outbound tasks under a constraint on the waiting time of inbound calls. We propose an efficient adaptive threshold policy easy to implement. This scheduling policy is evaluated through a comparison with the optimal performance measures found in the case of a constant stationary arrival rate, and also a comparison with other intuitive adaptive threshold policies in the general non-stationary case.

In the third study, we consider a call center model with a call back option, which allows to transform an inbound call into an outbound one. The optimization problem consists on minimizing the expected waiting time of the outbound calls while respecting a service level constraint on the inbound ones. We propose a routing policy with two thresholds, one on the reservation of the agents for inbound calls, and another on the number of waiting outbound calls. A curve relating the two thresholds is determined.