Staffing at Axxicom Airport Caddy Optimization with queueing theoretic models

(1)

MSc Mathematics

Master Thesis

Staffing at Axxicom Airport Caddy

Optimization with queueing theoretic models

Author: Supervisors:

Isabel Bevort

prof. dr. R. N´

u˜

nez Queija

dr. M. Remerova

Examination date:

(2)

Abstract

Axxicom Airport Caddy handles the transportation of Passengers with Reduced Mo-bility at Schiphol Airport, which is a service to passengers provided by the airport. A performance requirement of 98% is prescribed by Schiphol Airport regarding the on-time delivery of departing passengers at the gate of their flight. The goal of this thesis is to optimize capacity planning with the constraint of achieving this performance level. This is done by the development of a mathematical model that captures the relationship be-tween the inflow of passengers, the allocated capacity and the experienced performance. Several variations of queueing theoretic models capturing various characteristics of the actual logistic processes were investigated. The best-fitting model is based on an M/D/c queueing system that incorporates the practically highly relevant property that several jobs can be handled as one, thus tremendously improving the efficiency of the trans-port service. The sojourn time distribution was key in the mathematical performance analysis. The result of this project is a general framework that can be used to assess the correct level of capacity, accompanied with specific recommendations regarding the required minimal capacity per hour per weekday for three time periods depending on the season (high, low, mid).

Title: Staffing at Axxicom Airport Caddy – Optimization with queueing theoretic mod-els

Author: Isabel Bevort, isabel.bevort@student.uva.nl, 5881420 Supervisors: prof. dr. R. N´u˜nez Queija, dr. M. Remerova Internship Supervisors at Axxicom A.C.: D. Held, G. Dekkers Second Examiner: dr.ir. E. Winands

Examination date: August 30th, 2016

Korteweg-de Vries Institute for Mathematics University of Amsterdam

Science Park 105-107, 1098 XG Amsterdam http://kdvi.uva.nl

Axxicom Airport Caddy

Havenmeesterweg 43, 1118 CB Schiphol http://trigion.nl/airport-caddy

(3)

Acknowledgements

First and foremost I owe many thanks to my supervisors Sindo N´u˜nez Queija and Maria Remerova. They have patiently guided me in the right direction during our long but pleasant brainstorming sessions and supported me on occasions where I was feeling stressed. When programming errors led to me say things like “my probabilities are imaginary”, they encouraged me to try again and keep smiling. I have learned a lot from them as their guidance and supervision led me to the completion of this graduation project.

Furthermore I would like to thank Axxicom Airport Caddy for welcoming me into their offices. During the first months of the project, the office was located at the seventh floor, overlooking the airport and runways, which was a fantastic view and made me even more enthusiastic about the project. Next, the supervision of Gijs Dekkers and Desi Held has moved my research in the right direction. I owe thanks to Evert de Jong for providing me with the tools to extract the relevant data from AAC’s vast database and thinking with me about the correctness of certain assumptions. I would also like to thank Hubald de Heij, the other graduation intern at AAC and current employee, for the companionship and teaching me about airplane models. It was very pleasant to do this graduation project in an environment where my enthusiasm about the aviation industry and Schiphol Airport is shared.

Lastly I would like to thank the people in my personal environment. My loved ones have supported me, encouraged me, and thought with me. They have always trusted in my abilities, reassured me when it was necessary, and celebrated with me when I made improvements.

(4)

1. Introduction

As one of the major European hubs, the safe transportation of passengers and guests is a key priority of Schiphol Airport. In this respect, Passengers with Reduced Mobility (PRMs) deserve special attention. Travel assistance to PRMs at Schiphol Airport is provided by Axxicom Airport Caddy (AAC). This company has the challenging task of assisting up to 1500 passengers per day, of whom the nature of their reduced mobility varies from being wheelchair bound to being blind or deaf. Performance requirements at the level of 98% are imposed by Schiphol Airport regarding, for instance, the on-time delivery of passengers at their destination (e.g., departure gate). This graduation project was done as an internship for Axxicom Airport Caddy and its purpose is to develop a model with which we can optimize capacity for the departing passengers, considering the performance requirements and the inflow of passengers.

This chapter is devoted to the description of the problem that was investigated for this research project. We will start with a management summary aimed at the direction of Axxicom Airport Caddy. Note that the management summary may contain terms that are unfamiliar to readers outside Axxicom Airport Caddy, but which will be explained later. Further, we will give an overview of this thesis, followed by a precise problem description, an explanation regarding AAC’s key performance indicators (KPIs) and an overview of the available data that is provided to us by AAC.

1.1. Management summary

When the project started in February 2016, Axxicom Aiport Caddy was encountering an overcapacity while at the same time not meeting their performance targets imposed by Schiphol Airport. In this thesis we will report on the steps taken to optimize both capacity and performance in the work area Departures, including the mathematical work as part of this graduation project. Furthermore, some data analysis will be done for the two other work areas, Screened and Unscreened.

The start of the mathematical research was to analyze data from the work area Depar-tures in order to find the bottlenecks of AAC’s operations. It was important to recognize which variables would be relevant for the capacity requirements, making it possible to develop a mathematical model to mimic the working of Axxicom Airport Caddy’s busi-ness. From the data analysis we uncovered that the distribution of the capacity over the day was far from optimal. The busiest moments in this work area happen between 08:00 AM and 14:00 PM, distributed over two peaks with a small decrease in between at around 11:00 AM. However, the second peak was not always staffed accordingly, re-sulting in a disparity between inflow of customers and available capacity to handle the

(7)

customers appropriately.

Starting with a simple model, and expanding it gradually with the help of my supervisors both at Axxicom and the UvA and through literature analysis, we have developed a model that mimics the Departures work area well. There are some reservations that will be discussed later on, however, we feel confident that the description is quite accurate. This model has captured the dependence between the capacity, the inflow of passengers and the KPI-performance. In fact, the results of the model include a precise expression of the success probability of KPI-1, which is the indicator for AAC’s performance regarding the departing customers. Using this relation, we have optimized the capacity such that the performance of KPI-1 is at least 98%, which is the requirement from Schiphol Airport. We do this for three time periods with different inflows of customers: high season, low season, and in between. The result from this optimization is a concrete advice on the number of passengers assistants that should be scheduled at every hour of the day, for each weekday in these three periods. These results are enclosed in the document created for this purpose, see [4] in the reference list and appendix B.

Recommendations

The main results from this analysis are the following recommendations.

1. Firstly, we would advise AAC to reserve a part of their capacity for the pre-notified customers, which are the customers that meet the requirements to be included in the KPI-1 calculations. These conditions can be summarized as: the passengers need to be at the Axxicom desk on time and the processing of their task should consist of the transport from desk to gate; see section 3.3 for an elaborate description. Note that various but infrequent special cases are excluded by these requirements. All other tasks will be called ad hoc tasks. These tasks arrive at a far smaller rate at the desk than the pre-notified tasks. However, since they often represent customers that arrive quite late and need to catch their flight, they tend to get priority, leading to a decrease in available capacity for the tasks that influence this KPI.

2. Secondly, we would advise a different distribution of capacity over the day. On one hand, there should be more capacity during the second peak in the busiest time frame. On the other hand, keeping more capacity available after significant peaks will ensure that the backlog occurring from these peaks will not hinder the performance of KPI-1 in the period following them.

3. Thirdly, having a higher variability in service times than is currently the case would entail a limitation on the obtainable performance that is not dependent on the capacity. In other words, deploying more PAs may, in that case, not increase the KPI.

The third recommendation may be rather surprising. Axxicom Airport Caddy has com-mented as a reaction to this observation that they should therefore strive to keep the

(8)

All relevant results and interpretation are given in chapter 5. The remainder of the thesis will be structured as follows. The problem description is given in this chapter. The choice for a certain mathematical model is justified in chapter 2, where a short introduction on the mathematical theory is given. Chapters 3 and 4 contain the mathematical analysis concerning the model development. Chapter 6 contains the mathematical conclusion and the popular summary (meant for first year mathematics students) completes this thesis.

1.2. Overview

In this thesis we will investigate AAC’s performance concerning the departing passengers. The main mathematical model that we use is an extension of the standard M/D/c model, which will be explained in section 3.4. The expressions used are derived by Franx, see [10]. He derives a closed form expression for the waiting time distribution, enabling us to get a description of the success probability of the performance measure for the departing passengers. It is one of the six Key Performance Indicators (KPIs) that will be explained in section 1.4. The success probability for departing passengers is given by the probability that a departure deadline is not violated, denoted as P(S ≤ T ), where S represents the sojourn time and T represents the amount of time that is available to AAC to complete the job before the deadline. The state probabilities of the M/D/c are needed for this probability, and they are computed using the Fast Fourier Transform method [20, 21]. The arrival rate parameter is estimated using the Maximum Likelihood Estimator for the exponential distribution [23], incorporating the possibility that jobs may be combined. The data samples are obtained after analysis of data extracted from Axxicom Airport Caddy’s vast database.

The goal of this graduation project was to find an optimal distribution of capacity over the day, such that the performance requirement imposed by Schiphol Airport is met. By capacity we mean the number of deployed Passengers Assistants (PAs). The optimal distribution was obtained by extending the stationary M/D/c model that is developed in chapter 3 to a non-stationary model with the use of the effective load defined by Massey [13]. In the non-stationary model, see chapter 4, the success probability can be computed similarly to the success probability in the stationary model, with the addition that the queue’s ability to drain the workload is included.

The approach taken for this graduation project was as follows. We have started with data analysis in order to uncover the bottlenecks for AAC’s low performance, while simultaneously investigating suitable models through literature research. The distribu-tion of capacity throughout the day turned out to be of importance. The next step, in chapter 3, was to develop a stationary model that could mimic the situation at Axxicom Airport Caddy well. This large chapter describes the development of the stationary model and covers several models that were investigated for this purpose. The term “sta-tionary” refers to the assumption that we focus on a period of the day during which the characteristics of passengers flows do not fluctuate too much (we will incorporate such fluctuations later). The first model consists of a standard M/M/c queueing system

(9)

with exponential interarrival times and service time and c identical parallel servers. This model was improved upon by incorporating into the parameter estimation the fact that jobs may be combined. Next, we have derived the success probability of the performance measure under investigation. We found that the assumption of exponential service times did not describe the operational reality accurately enough. In fact, the results showed that it represented a lower bound on the performance that can actually be attained. Therefore, we moved on in section 3.4 to an M/D/c queueing model with deterministic service times, in which the combination of jobs was included as well. This model gave a good representation of AAC’s situation. Next, we chose to extend the M/D/c model to a non-stationary model in chapter 4 using the effective load as defined by Massey, see [13]. This model was then used for the optimization procedure, yielding a recommended minimal amount of capacity per hour for each weekday in three time periods, the re-sults of which are gathered in the internal company manuscript created for this purpose, see [4].

There have been several publications on PRM management at airports, see [5, 18, 2, 3]. To the best of our knowledge, no mathematical articles have been published on the development of a mathematical model for PRM management in the aviation industry. Before deciding on the queueing theoretical model, we have investigated several modeling possibilities, including networks of queues [15, ch. 5] and call center applications [11].

1.3. Precise problem description

In this section we will elaborate on the problem description. At the start of this project, Axxicom Airport Caddy experienced overcapacity while at the same time not meeting the performance requirements. It was not clear what caused this disparity. Therefore, it would be interesting to investigate the relevant variables in determining the level of required capacity, with which a model could be developed to mimic Axxicom’s situa-tion. Finding the relationship between the inflow of passengers, the capacity and the performance was crucial in being able to optimize both capacity and performance. This was done for the performance concerning the departing passengers only. For the arriving passengers, the performance is measured differently, which can possibly lead to interest-ing research for a next graduation thesis. We did include data analysis to provide AAC with interesting insights concerning the inflow of passengers over the day, the average duration of tasks, and the combination factor of passengers into one task.

Let us introduce Axxicom Airport Caddy in more detail. It has already been mentioned that this graduation project was done as an internship for this company, which is located at Schiphol Airport and handles the transportation of Passengers with Reduced Mobil-ity. The PRMs consist of various types of passengers, e.g. passengers in a wheelchair, Unaccompanied Minors (UMs) or blind passengers. The transportation of these passen-gers is done from gate to the Arrivals hall, from the Departures hall to the gate, or from gate to gate when it concerns a transfer passenger.

(10)

350.000 PRMs a year. Their performance is measured using six Key Performance Indi-cators (KPIs) which averaged around 92 % when this project started. Schiphol Airport has imposed a performance requirement of 98%, and AAC has accepted this graduation project as part of their efforts to improve the performance.

As mentioned, the travel assistance consists of many different types of transport inside the airport. The jobs that are generated by these PRMs are therefore divided over so-called work areas. The work areas are defined in such a way that there is hardly any overlap. The passenger assistants are assigned to one work area during the entire duration of their shift. It may occur that PAs are interchanged between work areas at the beginning of their shift, but in theory, once a shift has started, a PA will not do any jobs generated in another work area. Naturally, re-assigment of PAs at an operational level can still be an effective measure to counter unexpected loads; but such measures must by nature be taken on an ad hoc basis.

When this research project started there were four work areas:

• “Departures”: for departing passengers from either the reception desks (called “ASBA” for “ASsistentieBAlie”) or the buffer location (waiting area, called “Ser-vice Centrum West”);

• “Terminal”: for passengers that are either strictly arriving or transferring, coming from a “clean” origin;

• “Platform”: for both arriving and departing passengers, for the gates that are not connected (i.e., gates that have no jet bridge, which means the airplane is accessed through stairs);

• “Unscreened”: for strictly arriving passengers or transferring passengers, coming from a “dirty” origin; the arriving passengers are brought directly to the exit hall, but the transferring passengers need to go through a security filter before they can be brought to the gate of their next flight or to a buffer location.

The terms “clean” and “dirty” have originated from the installation of the OneXS secu-rity system, centralizing the secusecu-rity checks at Schiphol Airport. As mentioned, trans-ferring passengers from a “dirty” origin need to go through an additional security filter before being allowed onto their connecting flight.

For completeness, it should be mentioned that a new work area, called “India”, has been added in the mean time due to the great inflow of PRMs from the flights coming from India, and the addition of Jet Airways flights1 to the flight schedule.

This project will be restricted to the analysis of the work area Departures because we will focus on AAC’s performance concerning departing passengers. While this thesis includes partial data analysis for the Screened and Unscreend work areas, we leave the Platform work area out of consideration, mainly due to the fact that AAC will not take any risks on the platform and supplies a standard overcapacity there for this reason.

1

Jet Airways operations on Schiphol have started in March 2016, with daily flights to Delhi, Mumbai and Toronto

(11)

The Departures work area handles the passengers on local departing flights and the de-parture part of transfer passengers that are located at the Service Centrum West (SCW), the buffer location acting as the waiting area for passengers with a long connection time. Let us give a simplified version of the current staffing planning at AAC. For a given time period, the expected number of PRMs per day is based on the number of PRMs in that time period in the previous year, yielding expected staffing requirements as follows. Each PRM will need 55 minutes of service from one Passenger Assistant, so for instance, on a day on which 1000 PRMs are expected, they plan to have 1000×55 available service minutes. These 55 minutes include breaks for the PAs. However, the data analysis and experience have shown that the distribution of the capacity over the day is not optimal. In this graduation project, we have identified the relevant variables for the performance and staffing and investigated how to optimize staffing with respect to the available passenger flow data. The main goal in this optimization will be to obtain a performance of 98% for KPI-1, which we describe in more detail in the next section.

1.4. Key performance indicators

AAC measures its performance using six Key Performance Indicators (KPIs), which are success percentages. For completeness, we will introduce all six here, but only KPI-1 will be optimized in this thesis.

Firstly, for departing passengers, the performance is measured with KPIs 1 and 2. KPI-1 represents the percentage of PRMs that are delivered on time at the gate. A boolean is recorded for each passenger, with a 1 if the passenger is delivered on time, and a 0 if not. The notion of on-time delivery is dependent on the type of flight, either European (EUR) or Intercontinental (ICA), see below. As mentioned in the management summary, not all departing passengers are included in the KPI-1 calculations. See section 3.3 for the conditions that need to be satisfied.

1. KPI-1 records the performance for departing passengers that are pre-notified and whose booking times are less than 60 minutes before the STD2 for EUR flights and less than 90 minutes before the STD for ICA flights.

2. KPI-2 records the performance for departing passengers whose booking times do not meet these requirements.

Secondly, for transfer passengers, the performance is measured with KPI-3. It represents the percentage of PRMs that have caught the transferring flight.

Finally, for arriving passengers, the performance is measured with KPIs 4 and 5. They represent the percentage of times that a PA is present on time at the gate of the ar-riving flight. Here, on-time presence is dependent on the moment that the assistance is requested. For arriving flights this is measured using the creation time, which is the moment that the airline notifies AAC that they will be carrying passengers in need of assistance.

(12)

1. KPI 4 records the performance for arriving passengers that are pre-notified more than 15 minutes before the AIBT3.

2. KPI 5 records the performance for arriving passengers whose assistance require-ments were notified less than 15 minutes before the AIBT.

Note that transferring passengers are both departing and arriving as they have both an arriving and a departing flight, therefore three KPIs will be recorded: for the arriving flight either KPI-4 or KPI-5, for the departing flight either KPI-1 or KPI-2, and KPI-3. The sixth and last KPI gives the percentage of PRMs that have left the airport within two hours.

1.5. Available data

Axxicom Airport Caddy works with a software program called Groundstar, which was specifically created for them by INFORM. It is a scheduling program that dynamically assigns incoming tasks to the personnel. There are four types of shifts:

1. The PA (“PassagiersAssistent” or passenger assistant) is responsible for transport-ing the PRM from the start location to the end location. This is, for instance, a transport from reception desk to gate, a gate to the Arrivals hall, or a transport from an arriving flight to a connecting flight.

2. The TA (“TilAssistent” or lifting assistent) is responsible for the lifting part of passengers that are wheelchair-bound. This concerns the lifting in and out of the airplane seat.

3. The PLAM (“‘PlatformMederwerker” or platform employee) is responsible for the transport of PRMs whose flights depart from a gate that is not connected to the terminal by a so-called jet bridge.

4. The CHF (“Chauffeur” or chauffeur) drives the buses that transport the PRMs on the platform, which covers all gates that are not connected to the terminal by jet bridges.

The program Groundstar combines data from three parties:

1. From the digital recording device (PDA) that the PAs carry. The PAs are required to record multiple time stamps for each task that they complete. See below for more explanation.

2. From Schiphol Airport: the airport provides flight information data to AAC, which is recorded in Groundstar and used, among other things, to determine deadlines for the passengers.

3

AIBT: Actual In-Blocks Time, or ONB: On-Blocks. This is the moment that a flight has arrived and the airplane is being connected to the gate

(13)

3. From the airlines: the airlines are requested to send information about the passen-gers on board of their flights that require assistance. This information is sent for arriving flights.

Each PA carries a PDA during their shift with which they record multiple moments for each task. A task is defined as the transport of a PRM from one location to another. A unique passenger ID is assigned to each PRM, which is the same for each of the tasks that is generated by the passenger. One PRM can generate multiple tasks, for example a transferring passenger that is dependent on a wheelchair with a long connection time:

1. From the seat in the airplane to the gate,

2. From the gate to the buffer location (i.e., waiting area), 3. From the buffer location to the gate of the departing flight, 4. From the gate into the airplane seat.

Note that the lifting part of passengers is not done by a PA, but by a TA. This particular PRM has thus generated 4 different tasks, two of which are processed by a PA, and the other two by a TA. We have mentioned that only certain tasks are included in the KPI-1 calculations. In this example, the only task that would count for KPI-1 is the third task: transporting the passenger from the buffer location to the gate of their departing flight. Multiple similar tasks may be combined into one combined task. To each combination of jobs, a unique Combined Task Number is given as well, even when the task is not combined with another task. So each PRM has a unique passenger ID and may share a combined task number.

Next, each task consists of several steps, where each step is recorded. See figure 1.1 for the recorded time stamps of an example task that was handled on February 10thbetween 08:15 and 09:16.

Figure 1.1.: Time stamps recorded during the processing of one task The records that will be interesting to us are:

(14)

• Started: the PA has arrived at the start location and finds the PRM. The service of this job has started.

• Arr Delivery: the PRM has been delivered at the end location. This recorded moment is the deadline for KPI1.

• Finished: the Pa has finished the transportation of the PRM and is ready for the next job.

The time stamps Confirmed and Finished are interesting because the time between them is the total time that the PA is busy with this job and unable to process any other jobs. The time stamp Started is mostly interesting for the data analysis as it gives AAC an indication of how long the actual transport takes, and how long it takes the PA to arrive at the start location of the job. Note that Schiphol Airport is a large airport and it can take up to 30 minutes to walk from one side to the other side. Lastly, the Arr Delivery is interesting in the analysis as this is the moment recorded for the KPI-1 success or failure. Note that the time stamps Arr Delivery and Finished do not always coincide. This is due to the fact that the “unloading” of the PRM off the caddy or wheelchair may occur in between these two moments.

A tool4 was designed for the purpose of this thesis by Axxicom Airport Caddy, allowing us to extract the data for a given time frame, work area and type of PRM, for selected task statuses. We have chosen to export the data in Excel-format, with one row per changing task status. The excel-files contain a number of columns, with the most relevant being

• the Task Status Change Time, which allowed for the computation of task comple-tion duracomple-tions, for both each job and the combined jobs,

• the Combined Task Number, which allowed for the computation of the number of occurrences per Combined Task Number in a separate column,

• the Task Start time in combination with the Passenger Unique ID, which allowed for the derivation of the arrival intensities per hour (using the pivot table tool from Excel).

Regarding capacity, another tool was created by Axxicom Airport Caddy, allowing us to extract data about the shifts. This information is delivered in Excel-format with one row per shift, therefore it is known which PA’s had what shifts in which time periods. Through analysis on the Excel-file, this information has been used to determine the number of shifts per hour that was used, for each day of the week.

We will use data from August 2015 and later. The earlier data concerning task duration are not relevant in the current situation due to the introduction of the OneXS security system, which has centralized security at Schiphol Airport, and has led to the separation of jobs into different work areas.

(15)

1.6. Structure

Now that we have described the problem, the remainder of thesis is structured as follows: 1. Chapter 2 links the real situation at Axxicom Airport Caddy to the mathematical

modeling choice.

2. Chapter 3 is a large chapter describing the development of the stationary model which is used for the final model in the next chapter.

3. In chapter 4 we extend the analysis to a non-stationary model, which is then used for the capacity optimization.

4. Chapter 5 lists the main results and an interpretation of the optimization proce-dure.

(16)

2. Mathematical modeling

In this chapter we translate the problem description at Axxicom Airport Caddy into a mathematical framework. We will begin by giving a short introduction on queueing theory, after which we will look closer at the Departures work area and finish with a justification of the modeling choice.

2.1. Queueing theory

Let us give a short introduction to queueing theory. A basic queueing model is shown in figure 2.1. It represents the flow of tasks through a servicing system that is characterized by uncertainty in the arriving of tasks and uncertainty in the processing times of tasks. The actual flow further depends on the number of parallel tasks that can be performed and the rules according to which the order of service is determined. The basic queueing model is characterized by the following, see[1]:

• The arrival process of customers

It is assumed that interarrival times are independent and have a common distribu-tion. Oftentimes, customers arrive randomly in practical situations, in which case the Poisson distribution is used.

• The behavior of customers

Customers are either patient and willing to wait, or impatient and will leave after a certain amount of time passed since their arrival.

• The service times

In most models it is assumed that the service times are independent and have a common distribution, and that they are independent of the interarrival times. • The service discipline

There are many different service disciplines, but the most used are first-come-first-served, last-come-first-first-come-first-served, processor sharing, priorities.

• The service capacity

There can be either one single server or multiple parallel servers providing service to the customers.

• The waiting room

The number of customers in the system can have an upper bound due to a limited waiting area space.

(17)

Figure 2.1.: Basic queueing model

The typical notation for this type of models was introduced by Kendall and has the form a/b/c. The first letter characterizes the arrival process with the interarrival time dis-tribution, the second letter characterizes the service time distribution. The third letter specifies the number of servers. For example, the letter G is used for a general distribu-tion, the letter M is used for the exponential distribution (M stands for memoryless), and the letter D is used for deterministic times. The notation M/D/c then represents a model with exponential interarrival times, deterministic service times and c servers. This notation can be extended for other types of queueing models. For example, if the waiting room is limited, a fourth letter is added, so M/D/c/N is as above, but with room for only N customers, including those in service. The service discipline is usually mentioned with the model, as is the customers’ behavior.

The occupation rate ρ represents the fraction of time that the server is working. For stability, it is required to have ρ < 1. A system is stable when the queue cannot grow to infinity. The arrival rate is the average number of tasks arriving per unit of time. In a single-server system G/G/1 with arrival rate λ and mean service time E(B), the occupation rate or load is given by ρ = λE(B). In a multi-server system G/G/c, we require that λE(B) < c, and the occupation rate per server is ρ = λE(B)/c.

Let us define two further notions that are of special interest. Let the random variable L(t) denote the number of customers in the system at time t, and let Sn denote the sojourn time of the n-th customer in the system. When the system is assumed to be stable, i.e. the occupation rate per server is less than one, these random variables have a limiting distribution as t → ∞ and n → ∞, which are independent of the initial condition of the system, see [1, p. 25-26]. Let the random variables L and S have the limiting distributions of L(t) and Sn, respectively. Then

pk= P(L = k) = lim

t→∞P(L(t) = k), FS(x) = P(S ≤ x) = limn→∞P(Sn≤ x), where pk is the probability that there are k customers in the system, and FS(x) gives the probability that the sojourn time of an arbitrary customer entering the system is nog greater than x units of time.

(18)

entering the system per unit time is called Little’s Law. It states that E(L) = λ E(S),

under the assumption that the system is stable.

Before justifying our choice to model the work area as a queueing system, we will give several relevant performance measures in the analysis of queueing models:

• the distribution of the waiting time and the sojourn time of a customer, where the sojourn time is defined as the waiting time plus the service time,

• the distribution of the number of customers in the system,

• the distribution of the amount of work in the system, which is the sum of the service times of the waiting customers and the residual service time(s) of the customer(s) in service,

• the distribution of the busy period of the server, where the busy period is defined as a period of time during which the server is continually working.

2.2. From practice to mathematics: the work area as a

queueing system

This research project is focused on investigating the performance of the work area Depar-tures. The process of the Departures work area can be well described using a multi-server queueing theory model. Let us first describe the process of this work area.

As a passenger arrives at the desk to request their travel assistance, the time in which Axxicom can complete their assistance starts. AAC is assigned a deadline to complete this task, which is defined as the Scheduled Time of Departure (STD) minus 30 or 45 minutes, for European or Intercontinental flights, respectively. The next part of a passenger’s time in Axxicom’s system constitutes their waiting time, which means that they wait in the waiting area near the desk for their travel assistance to begin. There are multiple passenger assistants working at the same time. At the moment that one of them becomes available to start the transportation of a specific passenger, they begin the travel assistance for the passenger. Note that the time that a PA walks to the start location of the task – between time stamps Confirmed and Started – is included in the total assistance time. During the transport of the passenger, the passenger assistant records several moments in their PDA and completes the task when the passenger has been delivered at the end location (a gate) and is transferred into the care of the airline crew. If the passenger was delivered before the passenger’s deadline, a success is registered for the KPI-1 of this passenger. If not, it is recorded as a failure.

(19)

With this description it is natural to use a queueing model to mimic the work area. Indeed:

• there is an arrival process of customers (passengers),

• the customers are patient and will not leave no matter how long their waiting time is,

• the duration of the transport constitutes the service time of a customer, • the service discipline will be assumed to be first-come-first-served, • there are multiple PAs (servers) providing service to the customers, and

• the waiting room will be assumed to be infinitely large, which is not unrealistic given the magnitude of Schiphol Airport.

The following chapter will describe the process and result of the development of a sta-tionary model that mimics the Departures area on an hourly basis. In chapter 4 we will extend the model to a non-stationary version in order to get an optimal distribution of capacity over the day.

(20)

3. Stationary Markovian Analysis

This chapter is devoted to the description of the development of a suitable model through stationary Markovian analysis. The models in this chapter are all assumed to be sta-tionary, i..e, time-homogeneous Markov chains; the suitable model will be extended to a non-stationary model in chapter 4. In this chapter, the mean number of customers is used as a performance measure in the starting phases of model development, mainly to assess the models. The more advanced and therefore more suitable models will use the sojourn time distribution as performance measure.

In section 3.1, we begin with the most simple model: an M/M/c queueing model to which no extra assumptions have been added. Section 3.2 will have the added assumption that jobs may be combined included in the parameter estimation. The sojourn time distribution is added in section 3.3, with an interesting conclusion. This conclusion leads us to investigate a model with not exponential, but deterministic service times in section 3.4. This turns out to be the most suited model. This M/D/c queueing model includes the combination of jobs through the parameter estimation as well. This is the model that will be extended in chapter 4 to a non-stationary model, with which both capacity and performance will be optimized.

The models in this chapter will be tested using parameters that are estimated with the data from the hour between 08:00 AM and 09:00 AM on the Tuesdays of February 2016. We have chosen Tuesdays because Axxicom Airport Caddy has reported this to be the busiest day of the week for them.

3.1. Single job processing (model 0)

The first model we will investigate is the standard M/M/c model, a multi-server model with exponential interarrival times and exponential service times. The arrival rate is denoted by λ and the service rate by µ. In this first model we will assume that each passenger is handled by one PA, i.e., that the processing is done 1-to-1. We have al-ready noticed that this assumption is not accurate for the actual logistics of AAC, but this model will be instrumental in the development of more accurate ones, eventually culminating in the model of chapter 4.

Let us examine the adequateness of this model for the work area. In order to do so, we will determine its performance depending on the number of servers. As mentioned, for this first model we will use the mean number of customers as a performance measure. We suspect that having the expected total number of PRMs in the work area on average equal to the total number of PAs on shift in the work area will lead to a system that is performing well. Indeed, if that is the case, there will on average be no queue forming

(21)

while at the same time, no PAs will have to be idle. In mathematical terms, we will express the expected number of customers in the system, E(L), as a function of the number of servers in the system, Ec(L). We suspect that having the number of servers equal or close to a fixed point of the function E(L), i.e., Ec(L) ≈ c, will lead to a system that is performing well. We will compare the minimal number of servers for such a well-performing system to the number of PAs that were actually used at that moment. The number of customers in an M/M/c queue is equal to the sum of the number of customers in the queue and the number of customers in service. The expectation of the latter is given by the product of the occupation rate per server, ρ, and the number of servers, c. The former has an expected value that is given by the following expression [1, p. 44]:

E(Lq) = ΠW · ρ 1 − ρ,

where ΠW represents the probability that a job has to wait, given by Πw= pc 1 − ρ = (cρ) c c! (1 − ρ) c−1 X n=0 (cρ)n n! + (cρ)c c! !−1 .

Combining the above, we get the following expression for the mean number of customers in the system:

Ec(L) = ΠW · ρ

1 − ρ+ c · ρ. (3.1)

Note that all terms in this expression can be expressed in terms of the number of servers, c, and on the occupation rate per server, ρ. From the definition of this load, namely: ρ = λ

cµ, it follows that Ec(L) can be expressed in terms of c, λ and µ. The arrival and service rates will be estimated with the provided data. So, for a given arrival rate and service rate, the mean number of servers is indeed a function of only c, and we can identify the c that makes Ec(L) closest to c., i.e., the c for which Ec(L) and c close. 3.1.1. Parameter estimation

As mentioned, we assume that the interarrival times and the service times are expo-nentially distributed. The following two propositions are well-known results regarding parameter estimation using the Maximum Likelihood Estimator, see e.g. [23].

Proposition 1. Let X be an exponentially distributed random variable with parameter λ, and let x1, x2, . . . , xn be realizations from X, with sample mean ¯x =

1 n n P i=0 xi. The

(22)

Proposition 2. Let X be a random variable with a Poisson(λ) distribution. Given a realization of n measured values ki= 0, 1, 2, . . . for i = 1, . . . , n, the maximum likelihood estimator of the parameter λ is ˆλ = 1

n n P

i=1 ki

By assumption the system has service times that are exponentially distributed. The parameter can therefore be estimated with the MLE for the exponential distribution:

ˆ

µM LE := 1 ¯ x,

where ¯x is the sample mean of n independent samples from a exponentially distributed random variable. The data provides us with realizations of the service times.

By assumption the interarrivals times are exponentially distributed. This means that the arrival process is a Poisson process, i.e., the number of arrivals in a time unit has a Poisson distribution. Therefore, we can estimate the parameter λ with the Maximum Likelihood Estimator for the Poisson distribution:

ˆ λM LE := 1 n n X i=0 ki,

where a sample of n realizations ki = 0, 1, 2, . . . for i = 1, . . . , n of a Poisson distributed random variable is given. The data provides us with numbers of arrivals per hour. For the Tuesdays of February 20161, this yields the following estimators:

Hours λˆ µˆ Hours λˆ µˆ 0 2.00 0.96 12 35.25 1.27 1 1.50 1.04 13 27.00 1.24 2 0 0 14 15.00 1.22 3 0 0 15 14.50 1.40 4 1.25 1.28 16 11.25 1.48 5 10.25 1.17 17 9.25 1.36 6 9.25 1.45 18 13.75 1.21 7 18.50 1.31 19 17.50 1.37 8 29.75 1.37 20 9.75 1.18 9 23.75 1.25 21 2.00 1.35 10 31.75 1.49 22 1.50 1.77 11 32.00 1.31 23 3.25 1.18

Table 3.1.: Estimators ˆλ and ˆµ for each hour of the Tuesdays in February 2016

1

Recall that we have chosen Tuesdays as this is the busiest day of the week for AAC. We have tested the models on data from the month of February because it was available to us earliest.

(23)

It was already mentioned that for the analysis is this chapter, we will use data from the hour between 08:00 AM and 09:00 AM, for which the estimators are highlighted in Table 3.1: the arrival rate is ˆλ1 = 29, 75 and the service rate is ˆµ1= 1, 37.

3.1.2. Performance analysis

We have seen that Ec(L) depends on λ, µ and c, so we have defined a function expnocust (for expected number of customers) of c in Mathematica using equation (3.1), represent-ing the expected number of customers.

Note that since we are assuming stability, the load has to be less than 1, which poses an extra condition on the number of servers, indeed:

ρ = λ

c µ < 1 =⇒ c > λ µ.

We will use λ = ˆλ1 and µ = ˆµ1 from the previous section, hence we will only consider

c > ˆ λ1 ˆ µ1 = 29.75 1.37 = 21.7153.

Plotting this function on the domain yields the following graph:

22 24 26 28 30 c 10 20 30 40 50 Number of customers L expnocust(c) c Figure 3.1.: Graph of Ec(L)

Note that continuous values for c are (and will be in the following sections) allowed for ease of computation.This is possible since it is possible in Mathematica to evaluate a sum up to a non-integer summation bound: it will be summed up until the greatest integer less than or equal to that value.

(24)

The plot also includes the line L = c, since our goal is to find c such that of Ec(L) is closest to c. As we can see, Ec(L) has only one fixed point for c > 21.7153, which we can approximate using Mathematica’s FindRoot function2:

FindRoot[expnocust[c]==c, {c,25}].

This yields c → 24.7561. We can conclude that on the February Tuesdays mornings from 08:00 to 09:00 AM, if each PA handled only one PRM at a time, the expected number of customers in the system would be equal to the number of servers if this number of servers was equal to approximately 25.

3.1.3. Evaluation

Axxicom Airport Caddy has provided us with data concerning their capacity. In Febru-ary, there were on average 15.75 PAs working in the work area Departures on Tuesday morning in the hour between 08:00 and 09:00 AM, and the system still seemed stable. This is a rather large discrepancy with what model 0 suggests, with stability occurring for c > 21.7153. From this we can conclude that model 0 is not describing the real situation accurately enough, which is what we expected. In the next section we will improve this model by incorporating the possibility of combining jobs.

3.2. Processing with combined jobs (model 1)

In this section, we will improve the previous model by assuming that jobs may be combined (but do not need to be). This happens in reality when for example two PRMs that cannot walk long distances travel together or have neighboring departure gates: they need to or can be transported at the same time. These passengers can be handled by the same PA who will transport them on a so-called caddy, see figure 3.2.

Depending on the passenger arrival intensity, more or less jobs will be combined. We will incorporate this flexibility by adapting the load definition.

3.2.1. Adapted load definition

In order to incorporate the fact that jobs may be combined, let us see each combination of jobs as a single job. By this we mean that we will interpret the jobs such that there are only incoming single jobs, some of which will represent a combination of jobs, and some of which will represent a job that was not combined. For ease of discussion, let us call the jobs that represent a combination of jobs a combined job and the jobs that represent a job that was not combined a 1-to-1 job.

For the 1-to-1 jobs, there is no need to adapt the data concerning their service time in any way, since the duration of their service is just the end time minus the start time (which is at the time stamp “confirmed”). For the combined jobs we will adapt the data.

2

This function takes as input a function or function and {x, x0}, then searches for a numerical root to

(25)

Figure 3.2.: When multiple PRMs that are not capable of walking long distances need to be transported, a PA drives them around in a caddy.

Doing so will allow us to estimate the arrival and service rate of the incoming jobs, both combined and single, each seen as one job.

Service rate

From the data it is known what each individual job’s service duration (regardless of combination) is from the PA’s confirmation until the completion of service. We will define the duration of the service of a combined job as the duration between the first confirmation and the last completion. In other words, it will be the total time that the PA is busy with this particular combined job.

Arrival rate

There are two possibilities for adapting the arrival rate:

1. The first option is to count the incoming jobs, which are either 1-to-1 or combined, each of these being counted exactly once, and use the results from this count to estimate the arrival rate of single jobs.

(26)

2. The second possibility is to use the percentages of combinations, and use ˆλ1 from the previous section in an equation to account for these combinations: we put ˆ

λcombined = a1 · λ + · · · + an· _n1 · λ where a1 + · · · + an = 1. The combination terms are then ai·1_i · λ for i ≥ 2, where ai represents the fraction of jobs that is handled with a 1 to i processing, and a1 represents the fraction of jobs that is not combined, i.e., that are handled 1-to-1.

The advantage of defining the estimators with the second possibility is that we can decide how much jobs we want to be combined and by how much. It would make it possible to analyze for various choices what the effect of the load and the performance would be. However, it should be noted that it is unclear what the second possibility would entail for the assumption of Poisson arrivals. We will therefore only investigate the first adaptation of the arrival rate in this analysis. The consequences of the second possibility will be discussed in section 6.2.

3.2.2. Parameter estimation

The data regarding the number of incoming single jobs (again, these are either combined or 1-to-1), will serve as observations for estimating the Poisson parameter for the arrival process. For the service rate, we have computed the duration of each single job. For 1-to-1 jobs, nothing changes. For the combined jobs, we will use the aforementioned total time one PA is busy with the combined job, i.e., the duration between the first confirmation until the last end. For example, let us say that a PA has processed a combination of three jobs, for which we know the following durations:

Confirmation End

11:00 11:20

11:00 11:40

10:50 11:15

The service time of this combined job is then the duration between 10:50 and 11:40, i.e., 50 minutes, i.e., 0.83 hours. We have done this for all combined jobs in the data set, yielding the following estimators:

(27)

Hours λˆ µˆ Hours λˆ µˆ 0 1.25 1.04 12 15.75 1.27 1 0.50 0.96 13 14.25 1.26 2 0 0 14 9.75 1.35 3 0 0 15 6.75 1.37 4 1.00 1.42 16 8.00 1.47 5 6.25 1.15 17 6.00 1.52 6 7.00 1.44 18 7.75 1.09 7 11.50 1.31 19 9.00 1.39 8 15.75 1.30 20 6.00 1.23 9 13.00 1.28 21 1.25 1.23 10 11.00 1.27 22 1.50 1.77 11 14.50 1.24 23 2.25 1.33

Table 3.2.: Estimators ˆλ and ˆµ for each hour of the Tuesdays of February 2016

Again, we start by analyzing the hour between 08:00 and 09:00 AM for the Tuesday mornings from Feburary 2016, for which the estimators are highlighted in table 3.2: the arrival rate is ˆλ2 = 15.75 and the service rate is ˆλ2 = 1.30. As we can see, in this case the arrival rate of jobs is much smaller than the arrival rate in the previous section, ˆ

λ1 = 29, 75, when combination of jobs was not included yet. 3.2.3. Performance analysis

For this section we will do the same computations as we did in subsection 3.1.2. We will now use λ = ˆλ2 and µ = ˆµ2. The extra condition imposed by the stability assumption (ρ > 1) leads in this case to c > λˆ2

ˆ

µ2 > 12.1155. With these parameters, plotting

(28)

15 20 25 30 c 10 20 30 40 Number of customers expnocust(c) c Figure 3.3.: Graph of Ec(L)

The fixed point is then found in a similar way as in section 3.1.2 by using Mathematica’s FindRoot function

FindRoot[expnocust[c]==c, {c,15}],

which yields c → 14.4015. The number of customers is therefore closest to the number of servers when there are approximately 14 servers.

Hence, including the possibility of job combination in the parameter estimation yields a number of servers that is much closer to the actual number of servers that was used, which was 15.75 PAs.

3.2.4. Evaluation

With model 1 we have used the performance measure of the expected number of cus-tomers with the incorporation of combined jobs. The result suggests that it could be a good model since the difference between the actual number of servers and the capacity that the model suggests is not too great. Let us examine the sojourn time distribution in the next section with model 2.

3.3. Sojourn time distribution (model 2)

In this section we will elaborate on the previous model by adding and computing an additional performance measure: the sojourn time distribution. With this performance measure we will be able to compute the probability that a customer stays longer in the system than the time available to process the job before the deadline expires. Note that since we are elaborating on the model from the previous section, we will use the word “job” both for combined and 1-to-1 jobs, and count each of these as a single incoming job.

(29)

In this section, we will estimate the parameters using a different data set. We will compute the new performance measures for a subset of all jobs, namely for all jobs that satisfy the requirements to be included in the KPI-1 measurement, see section 1.4 and below in section 3.3.1.

It makes sense to consider only these jobs, because we want to know the effect of the capacity on the performance of KPI-1, which is measured only for the jobs that satisfy all of the conditions. Not excluding the data that is not considered for the performance calculation could entail that the estimators are negatively influenced by the jobs that do not satisfy the conditions. We will elaborate on this in subsection 3.3.1 It should be kept in mind that a certain capacity will need to be added on top of this for the processing of jobs that do not satisfy the KPI-requirements. This will be discussed in section 3.4.3. The performance measure that we will investigate in this section can be determined from a well-known result for M/M/c queues, see [1, p.46], and is oftentimes used to assess the performance of queueing systems. It makes sense to use the sojourn time distribution P(S > t) because evaluating it for t = T , the available processing time dictated by the deadline (see subsection 3.3.2), gives us the probability that a job exceeds its deadline, i.e., it represents the KPI-1 failure probability.

We will begin with the parameter estimation, taking into account the possibility to combine jobs, before applying the estimated parameters to determine the sojourn time probability and the probability of deadline violation.

3.3.1. Parameter estimation and comparison

Let us give the conditions for jobs to be included in the KPI-1 calculations. There are six of such conditions:

1. The creation time of the jobs needs to be at least 1 hour before STD3for European flights and at least 1.5 hours before STD for intercontinental flights, where the creation time is reasonably self-explanatory: it is the moment that the assistance was first requested.

2. The booking time of the jobs needs to be at least 1 hour before STD for European flights and at least 1.5 hours before STD for intercontinental flights, where the booking time is the moment that the passenger has reported him-/herself at the desk4.

3. The passenger type needs to be “PRM”.

4. The final task status needs to be “FINISHED”, meaning that the job needs to be finished at the time of KPI calculation. This condition is important for the daily KPI calculations.

(30)

5. The pick-up location needs to be different than the drop-off location, meaning that the job should not consist of gate boarding assistance (i.e., assisting the passenger onto their seat).

6. The end location needs to be a gate, meaning that the job should not consist of assistance on the platform.

When all of these conditions are met, the job is included in the KPI-calculation. From the 680 jobs in February 2016, 496 satisfied all of these conditions. With these jobs, we estimated the parameters similarly as before with the maximum likelihood estimators to obtain:

ˆ

λgood= 13.00 (3.2)

ˆ

µgood= 1.2730 (3.3)

In order to be able to do comparisons that make sense, we will repeat the numerical experiment in section 3.2.3 with the estimated parameters in this case. Plotting the expnocust function with these parameters yields the following graph:

15 20 25 30 c 10 20 30 40 Number of customers expnocust(c) c Figure 3.4.: Graph of Ec(L)

The fixed point is then at c → 12.3156. Therefore, the number of customers is then closest to the number of servers when there are approximately 12 servers.

We can summarize the results for the February Tuesday mornings from 08:00 to 09:00 AM in the following table:

(31)

Model Dataset Value of c such that Ec(L) is close to c

Model 0 All jobs (no jobs are combined) 25

Model 1 All jobs (combined jobs) 14

Model 2 Only the good jobs 12

Table 3.3.: Results for the Tuesdays from February in the hour between 08:00 and 09:00 AM.

3.3.2. Probability of missing a deadline

In this section we will derive the probability of missing a deadline, i.e., the KPI-1 failure probability. The distribution of the sojourn time in an M/M/c model is known and is given by the following expression, see [1, p. 46]:

P(S > t) = ΠW 1 − c(1 − ρ) e cρ(1−ρ)t₊ 1 − ΠW 1 − c(1 − ρ) e−µt. We need to compute this probability for t = T , where T is defined as follows:

Definition. The available processing time T represents the amount of time that the server has available to complete the job before the deadline expires, i.e., it is the duration between booking time and the deadline.

The notion of available processing time is inspired by Lehoczky’s lead-time in [12], which is defined as the difference between the customer’s absolute deadline and the current time. In his article, the lead-time is dynamic and decreases with time. Such dynamical data is not available to us. However, we do know the moment that customers (physically) enter the system, which is recorded as the booking time. This time stamp is recorded for each passenger as they arrive at the desk to request their assistance, whether it was pre-notified or not. We have therefore chosen for this deterministic version of the lead-time, which we call the available processing time. Note that transfer passengers do not have a booking time as they do not go to a desk first. For these passengers, Axxicom has agreed that the following moments are sensible replacements for their “booking time”:

• If the STA5 _{of their arriving flight is after the ONB}6 _{time, then we use the STA} as booking time;

• If the STA of their arriving flight is before the ONB time, i.e., the flight was delayed, then we use ONB + 15 minutes as booking time.

We are still only interested in the hour between 08:00 and 09:00 AM on the Tuesdays in February. During this period on those days, the average available processing time was 2 hours, 5 minutes and 48 seconds: T = 2.0967 hours.

(32)

Notice that all terms in the expression for P(S > t) depend on λ, µ and c. Therefore, we can, again, evaluate this expression as a function of c, where we use the estimated values ˆλgood and ˆµgood. Using the stability condition ρ < 1, this restricts the domain to c > 10.2121. 10 12 14 16 18 20 22 24 c 0.2 0.4 0.6 0.8 1.0 probsojourntime P(S>T)

Figure 3.5.: Probability P(S > T ) depending on c

We see that this probability decreases for increasing c, but seems to stagnate at a certain level. In fact, this level corresponds to the limit of the probability for c → ∞. This makes sense because when there are infinitely many servers, there is always a server available for each arriving customer. Hence, there is no waiting time and the sojourn time is then equal to the service time. For the exponential service times with parameter ˆµgood, this gives the following limit:

lim

c→∞P(S > T )(c) = e

−ˆµgoodT _{= 0.0693146.} _(3.4)

This result suggests that it will be impossible for Axxicom to have a performance of 98% on average (the maximum is at 93% approximately). Indeed, the above probability describes the fraction of passengers that miss the deadline. The goal is to obtain a performance of 98%, however we found that it can only get below 0.02 when we adapt the available processing time or the service rate.

This suggests that it will not be possible for Axxicom to obtain a performance of 98%, unless the PRMs enter the system at about 3.3 hours before the deadline (i.e., their booking time needs to be 3.3 hours before the deadline, which is then about 4 hours before STD).

In the next section we will evaluate and comment on this result. 3.3.3. Evaluation

Model 2 seems to give a good representation of the real situation at Axxicom. However, we have made assumptions on the interarrival times and service times that need to be

(33)

verified. The assumption of Poisson arrivals is quite reasonable as the passengers appear to arrive randomly. We do not have the data to check this, so we are limited to verifying the assumption on the service times. In particular, the assumption of exponential service times needs to be checked. When we plot the service times of the Tuesdays in February (for the entire day so that we have enough data points), we see that the exponential distribution is not a very good fit. Numerically, it seems that the model is a worse case than the data shows us. We have tried to fit the Erlang-k distribution to the service time distribution, however, this was not successful and we will discuss it further in section 6.2. The service times are plotted in the following bar chart, see Figure 3.6.

Figure 3.6.: Bar chart of the task durations from the 496 KPI-1 jobs of the Tuesdays of February 2016

As we can see in figure 3.6, most of the service times of the 496 good jobs are between 30 and 80 minutes, increasing slowly but steadily. There are some exceptions but these rep-resent a very small portion of the total. Therefore we have decided to use the assumption of deterministic service times, which is the topic of the next section.

3.4. Deterministic service times (model 3)

As we have seen in section 3.3.3, the assumption of exponential service times leads to some results that do not correctly reflect reality. In this section we will therefore model the Departures work area as an M/D/c queueing system, i.e., we assume that the service times are deterministic. This assumption makes sense as we have uncovered throug data analysis that the service times do not vary too much over the course of the day. There are some exceptions in the very early and very late hours, but these do not represent a big portion of the jobs. Even though it is not a perfect description of the reality because the service times are not exactly deterministic, this assumption will prove to mimic reality the most accurately of all the models we have studied and will therefore be extended to a non-stationary model in chapter 4. It may be worthwhile to consider Erlang-k distributed service times, but that did not fit into this project. We will comment a bit more on this later.

(34)

and compute the probability that a job’s sojourn time is bigger or smaller than the average available processing time. This probability represents the failure probability, respectively success probability, and describes the KPI-1 performance well.

We will first describe the method to obtain this probability, which is more involved than for the M/M/c model.

3.4.1. Method

For the M/D/c model, the occupation rate per server is defined as ρ = λD

c ,

where λ is the arrival rate, D is the deterministic service time and c is the number of servers. For stability we again require ρ < 1, i.e. λD < c.

There are several relevant results for the M/D/c queueing model. Crommelin described an exact algorithm analysis in 1932, see [7]. We will use the explicit expression for the waiting time distribution given by Franx in his 2001 paper on the M/D/c queue, see [10]. Since the service times are deterministic by assumption, we can derive the sojourn time distribution, which is given by the following expression:

P(S ≤ t) = P(W + B ≤ t) = P(W + D ≤ t) = P(W ≤ t − D). The waiting time distribution is given by

P(W ≤ x) = e−λ(kD−x) kc−1 X j=0 Qkc−j−1 λj(kD − x)j j , for (k − 1)D ≤ x ≤ kD. This is an explicit expression, except for the quantities Qm, with m = 0, 1, 2, . . . . Those are given by Qm =Pmi=0qi, where qiis equal to the probability that there are i customers in the queue. Let pi denote the state probability, i.e., the probability that there are i customers in the system. Then q0=Pc_j=0pj and qi= pi+cfor i > 0. Franx shows in his paper, see [10, p. 3] that these state probabilities are given by the following expression:

pj = c X k=0 pk (λD)j j! e −λD₊ j+c X k=c+1 pk (λD)j−k+c (j − k + c)!, for j = 0, 1, 2, . . . (3.5) This yields an infinite system of linear equations that can be solved in several ways: the geometric tail approach and the Fast Fourier Transform method are both quite efficient methods, see [22]. We will solve this system using the discrete Fast Fourier Transform (FFT) method as described by Tijms in both [20, p. 29] and [21, p. 457].

Fast Fourier Transform method

This subsection is devoted to describing the Fast Fourier Transform method that will yield the state probabilities, see [20, p. 29]. Let P (z) = P∞

(35)

function of the probabilities pj. When an explicit expression for P (z) is available, the discrete FFT method is a powerful method to recover numerically the values of the unknown probabilities pj, see [21, p. 455]. Let us describe how the FFT method computes the probabilities from the generating function.

1. First, we choose an integer N such that the tail probability P∞

j=Npj is negligibly small and so that

P (z) ≈ N −1

X

j=0

pjzj, |z| ≤ 1. (3.6)

It is required for the FFT method that N = 2m for some integer m. 2. Next we define the complex numbers

fk= P (e2πik/N), k = 0, 1, . . . , N − 1,

with i2 = −1. With the assumption that an explicit expression for P (z) is available, the numerical computation of the fk’s is not demanding.

3. Using (3.6), it follows that

fk≈ N −1

X

j=0

pje2πikj/N, k = 0, 1, . . . , N − 1. (3.7)

4. The FFT recovers the probabilities pj from the fk’s by using the inverse formula

pj = 1 N N −1 X k=0 fke−2πikj/N, j = 0, 1, . . . , N − 1.

5. The inverse formula follows from the fact that the matrix WN whose (n, k)-th element equals e2πikj/N can be shown to have (1/N ) ¯WN as inverse matrix, where the elements of WN are the complex conjugates of the elements of WN. The linear equations in (3.7) can be conveniently written in matrix notation as

f = p · WN, and hence p = f · 1 NWN = 1 Nf · WN.

6. The key explanation for the speed of this method is the observation that the discrete Fourier transform of length N can be written as the sum of two discrete

(36)

Fourier transforms, each of length N/2: N −1 X k=0 fke−2πikj/N = N/2−1 X k=0 f2ke−2πi(2k)j/N + N/2−1 X k=0 f2k+1e−2πi(2k+1)j/N = N/2−1 X k=0 f2ke−2πikj/(N/2)+ wj N/2−1 X k=0 f2k+1e−2πikj/(N/2),

where w = e−2πi/N. This can be implemented recursively, which explains the requirement of N being equal to 2m for some integer m.

The recovery of the probabilities as done in step 4. has actually been implemented in Mathematica using the very efficient functions Fourier and InverseFourier. This implementation has been done in such a way that it is no longer required for N to be a power of 2. The remaining steps were implemented by us, see appendix A.1.

In order to use the FFT method, the generating function must be in the form of an explicit expression. Given the form of the state probabilities in (3.5), this requires some extra work.

Generating function

In this subsection we will derive an explicit expression of the generating function of the state probabilities. By definition, the generating function is given by

P (z) = ∞ X j=0 pjzj = ∞ X j=0 e−λD(λD) j j! c X k=0 pk+ j+c X k=c+1 pke−λD (λD)j−k+c (j − k + c)! ! · zj = c X k=0 pke−λD ∞ X j=0 (λD)j j! z j ₊ ∞ X k=c+1 ∞ X j=k−c pke−λD (λD)j−k+c (j − k + c)!z j = c X k=0 pke−λD(1−z)+ ∞ X k=c+1 ∞ X `=0 pke−λD (λD)` `! z `−c+k = c X k=0 pke−λD(1−z)+ ∞ X k=c+1 pkzk e−λD(1−z) zc = c X k=0 pke−λD(1−z)+ P (z) −Pc k=0pkzk zc_eλD(1−z) .

(37)

Hence P (z) 1 − 1 zc_eλD(1−z) = Pc k=0pkzc zc_eλD(1−z) − Pc k=0pkzk zc_eλD(1−z) P (z) ·z c_eλD(1−z)_{− 1} zc_eλD(1−z) = Pc k=0 pk(zc− zk) zc_eλD(1−z) P (z) = Pc k=0 pk(zc− zk) zc_eλD(1−z)_{− 1} P (z) = Pc−1 k=0 pk(zk− zc) 1 − zc_eλD(1−z) . (3.8)

This expression still involves the unknowns p0, p1, . . . , pc−1, so we cannot use it for the FFT method just yet. Tijms describes how to obtain an explicit expression for P (z) from (3.8), see [21, p. 380]. The first step is to compute the zeros of the denominator 1 − zceλD(1−z) in the region |z| ≤ 1 in the complex plane.

Proposition 3. The 1−zceλD(1−z)= 0 has exactly c distinct zeros z0, z1, . . . , zc−1inside or on the unit circle, with z0= 1.

Proof. The proof uses Rouch´e’s theorem which we will state here, see [9, p. 105]: Theorem 1. Given two functions f and g analytic in A with C the boundary of A. If |g(z)| < |f (z)| for all z ∈ C, then f and f + g have the same number of roots inside C. Let f (z) = −zceλD(1−z) and g(z) = 1. Let D = {z | |z| ≤ 1}, and so C = {z | |z| = 1}. Since f (1) = −1, it does not hold that |f (z)| > |g(z)| for all z ∈ C. Define

D= {z | |z| ≤ 1 + },

where the value of will be determined later on. We will show that |f (z)| > |g(z)| = 1 for all z ∈ D. First note that |zc| > 1 for all z ∈ D. Suppose z = 1 + for some . We will determine which values of we want to consider. Then

|zceλD(1−z)| ≥ zce−λD = (1 + )ce−λD When = 0, it holds that (1 + )ce−λD = 1. Also, we have that

d d (1 + )ce−λD= c(1 + )c−1e−λD− λD(1 + )ce−λD = (1 + )c−1e−λD(c − λD(1 + )) (3.9) Now, if = 0, we get d (1 + )ce−λD = 1 · 1 · (c − λD) > 0.

Staffing at Axxicom Airport Caddy Optimization with queueing theoretic models

MSc Mathematics

Master Thesis