Using Differential Evolution Method to Solve Crew Rostering Problem

Airline crew rostering is the assignment problem of crew members to planned rotations/pairings for certain month. Airline companies have the monthly task of constructing personalized monthly schedules (roster) for crew members. This problem became more complex and difficult while the aspirations/criterias to assess the quality of roster grew and the constraints increased excessively. This paper proposed the differential evolution (DE) method to solve the airline rostering problem. Different from the common DE, this paper presented random swap as mutation operator. The DE algorithm is proven to be able to find the near optimal solution accurately for the optimization problem. Through numerical experiments with some real datasets, DE showed more competitive results than two other methods, column generation and MOSI (the one used by the Airline). DE produced good results for small and medium datasets, but it still showed reasonable results for large dataset. For large crew rostering problem, we proposed decomposition procedure to solve it in more efficient manner using DE.


Introduction
Development of crews rostering plan which be able to produce the high utility of crews become the priority in human resources department in airline industry.It is estimated that the use of optimization software for airline could save more than US $20 million per year [1].Saving 1% in crew utilization can save cost largerly.Though airline crews scheduling became attention in many operation research literature such as [1][2][3][4][5][6][7] but airline crews scheduling remains to become the main attention for many researchers due to its level of complexity and difficulty to solve.Therefore, methods and approaches which are used to solve it are continuously developed to get better result both in optimality side and speed of computational time.Generally, solving airline crew scheduling is done by decomposition approach [8][9][10]), it devides problem into crew pairing and crew rostering.Crew pairing is done to get initial feasible solution, that is sequence of flight which begin and end at the same home base.Crew rostering assigns pairings which were arranged for the certain month to set of crews based on individual calender.
Decomposition approach is very effective to solve the difficult and complex problem but this method loss the global treatment since crew pairing and crew rostering done separately.Some other researchers developed the integrated approach to overcome obstacle, such as Souai and Thegem [5], where crew pairing and rostering were done simultantly to get a better optimality level.Many optimization methods have been developed to solve crew scheduling to increase roster quality and to improve computational time such as simulated annealing [11], genetic algorithm [5], tree search algorithm [12], hybrid genetic algorithm [13] and GASA hybrid algorithm [14].
This research focused on developing differential evolution algorithm applied on intelligent airline crew rostering system.This paper is organized as follows.The second section reviews differential evolution (DE).Section 3 describes the problem statement.In Section 4, we describe our methodology.Section 5 explaines the experimental setting and the results.Section 6 concludes the results.

Differential Evolution
Differential evolution is an evolutionary population-based algorithm proposed by Storn and Price [15,16].Since its initiation in 1995, DE has shown its performance as a very effective global optimizer.DE originated with Genetic Annealing (GA) Algorithm.Since GA was very slow and effective control parameters were hard to determine, the modification of the GA algorithm were made.DE uses a floating-point instead of bit-string encoding and arithmatic operations instead of logical ones.DE differs significantly from the evolutionary algorithms in the sense that distance and direction information from the current population is used to guide the search process.DE uses the differences between two randomly chosen vectors (individuals) as the base to form a third vector (individual), referred to as the target vector.Trial solutions are generated by adding weighted difference vectors to the target vector.This process is referred to as the mutation operator where the target vector is mutated.The next step is recombination or crossover which is applied to produce an offspring.This new individual is accepted only if it improves on the fitness of the parent individual.The basic DE algorithm is described in more detail below [16].

Initialization
In this step, a set of initial solutions are genereted randomly.A random number generator assigns each variable of each vector value from within specified range, lower bound, b L , and upper bound, b U .For example the initial value (g = 0) of the j th variable of i th vector is: where, rand j (0,1), returns a uniformly distributed random number from within the range [0,1].

Mutation
DE mutates and recombines the population to produce a population of N-trial vectors.In particular, differential mutation adds a scaled, randomly sampled, vector difrence of third vector.To combine three different ranmly chosen vectors to create the mutant vector, v i,g , the following equation is used: where the scale factor, F (0, 1+) is a positive real number that control the rate at which the population evolve.The vector index, r 0 , r 1 , and r 2 , can be chosen randomly and meet r 0  r 1  r 2  i.

Crossover
DE employs the uniformly crossover.Crossover builds trial vectors out of parameter value that have been copied from two different vectors.In particular, DE crosses each vector with a mutant vector to creat u i,g . , The crosover probability, Cr ∊ [0,1], is user-defined value that controls the fraction of parameter value that are copied from the mutant.

Selection
If the trial vector, u i,g , has an equal or lower objective function value (better fitness value) than that of its target vector, x i,g , it replaces the target vector in the next generation; otherwise, the target retains its place in the population at least one more generation.
, , , , , Once the new population created, the process of mutation, recombination, and selection are repeated until the optimum solution is achieved or prespesified termina-tion criterion is satisfied, e.g., the number of generations reaches preset maximum, g max .
DE has been applied in many field successfully.In 1995, DE has been used by Ken to solve 5-dimension Chebyshev model.By the time, Ken modified genetic annealing algorithm with differential mutation operator.Different from genetic annealing, DE has not found some difficulty to find the coefficient even 33-dimension Chebyshev.
Tasgetiren [15] used a discrete differential evolution (DDE) algorithm to solve the single machine total earliss and tardiness penalties with a common due date.A new binary swap mutation operator called Bswap is prented.In addition, the DDE algorithm is hybridized with a local search algorithm to further improve the performance of the DDE algorithm.The performance of the proposed DDE algorithm is tested on 280 benchmark instances ranging from 10 to 1000 jobs from the OR Library.The computational experiments showed that the proposed DDE algorithm has generated better results than those in the literature in terms of both solution quality and computational time.
A genetic differential evolution (GDE) was derived from the differential evolution (DE) and incorporated with the genetic reproduction mechanisms, namely crossover and mutation used to solve traveling salesman problems (TSP).The Greedy Subtour Crossover (GSX) was employed to generate an offspring to denote the difference of the parents.A modified ordered crossover (MOX) was employed to perform mutation to generate trial vector with a user defined parameter, the parameter were used to control the rates of the target vector components and the mutated vector components in the trial vector.Moreover, a 2-opt local search was implemented to enhance local search performance.GDE was implemented to the well-known TSP with 52, 100 and 200 cities with variable parameters.Based on analysis and discussion on the results, typical values of the parameters were given, with which GDE provided effective and robust performance [18].
Omran and Salman [19] improved Differential evolution by combining with chaotic search, opposition-based learning, and quantum mechanics, called CODEQ, to solve constrained optimization problems.The performance of the proposed approach when applied to five constrained benchmark problems is investigated and compared with other approaches proposed in the literature.The experiments conducted show that CODEQ provides excellent results with the added advantage of no parameter tuning.

Problem Statement
There are two main processes in the airline crew planning.These two process are pairing and rostering.Pairing, is the step where the flying activities are created.The flight timetable is used as input to form sequences of flights, known as pairings.The timetable horizon usually covers a period of 4-6 weeks.The main objective of this process is to utilize the minimum number of crew to cover the complete timetable.Rostering is the step where the created pairings are assigned to actual crew (pilot and stewardess), with regard to the qualifications and previously assigned activities, referred to as pre-assignments, of the crew.The objective is to find feasible assignments that minimize costs and match the notion of quality of life for the crew imposed by the airline [9].In this context, problem of airline crew scheduling generally varies among different countries.Especially in rostering problem since each country may impose different set of regulations, rules and policies.The speed of roster construction is a critical matter in airline crew schedule- ing.In [9], it is stated that for monthly planning, solution should be obtained within 15-20 minutes.While, for shorter horizon planning, such as daily planning, where there are some changes to roster, solutions should be obtained within 1-5 minutes.Therefore, solving optimization problem of airline crew scheduling in time manner become very important.
In this paper two main problems are addressed.First, how to mathematically formulate the problem of airline crew scheduling.This formulation includes the construction of objective function and spesific set of constraints which influence the roster quality.We modified the model which is developed by Lučić and Teodorović [11] based on the real condition of rostering in the airline under study.We modified the single crew and single aircraft scheduling model becomes multiple crews and multiple aircrafts model.We also added open time criteria to the objective function.Second, how to solve the model using differential evolution approach with considering the simplicity of the model, quality of the resulting roster, and computational time.

Index
k is index for kind of crews k = 1, ,K.For example, k = 1 is for Pilot F-100; k = 2 is for Pilot CN-235; and k = 8 is for stewardess and etc.
i is index for numbers of crew members (1, ,m k ).For example m 5 = 17 is the number of crew members for Boeing 737-200, and m 8 = 55 is the number of stewardess of Boeing 737-200.j is index rotation/pairing which assigned to crew members (1, ,n k ).For example n 5 = 82, is the number of pairings for Being 737-200.
l is number of days in a month (1, ,31).

Parameters
d jk is length of rotation-j which assigned to crew k. (in hours).
1,if member from crew can be assigned to day 0, otherwise

Objectives Function
The objective function of this airline crew rostering is minimizing three terms of criterias.

Cost of roster
Cost of roster paid by the airline company to crew is variable cost.By assumption that salary per hour is same to all crew, cost of roster can be represented by actual flying hours.

Deviation of flying days between crew members
Let k t be the average flying days per month for crew member k, then The deviation of total flying days per month can be formulated as where p is positive integer.In this paper we use p = 1.

Open time
Open time is days when a crew member does not have flying duty.If there are 31 days in a scheduling month, then open time for crew member k can be formulated as

Constraints
There are some constraints which must be satisfied when constructing a roster.The following are the constraints used:

Flight time constraint
Maximum flying hours for pilot and co-pilot is 110 hours per month, and for cabin crew is 120 hours per month.So d max,k = 110 for k = 1, ,7 and

Duty period constraint
Maximum duty period allowed to crew member k is 21 days.

Numbers of take-off
Numbers of maximum take off allowed to pilot is 90, then v max,k = 90.But, cabin crew have no take off constraint.

Numbers of crew reqirement
Every rotation needs minimum numbers of crew.

Rotation without free day
When crew members complete this rotation not allowed to have free day.
No overlap constraint Two rotations in series may not be overlap each other.It means that precedence rotation must be finished when following rotation will start.
Airline rostering problem is optimization problem with many constraints.It falls to constrained optimization problem.To use DE to solve this original problem we have to transform the problem into an unconstrained optimization problem.We use external penalty function method [20,21] to do this transformation.Basically, this method incurs big penalty while the solution violated any constraints.The resulting constrained optimization problem is as follows where: where β 1 , β 2 , and β 3 are weight coeficients of the objective function, r 1 ,…,r 8 are penalties which are given since model violate any constraints where r 1 ,…,r 8 → ∞ will assure that algorithm will satisfy constraint early before considering objective function.Equation ( 16) becomes the fitness function of DE method.If we assume cost of roster more important than cost of deviation of flying day and more important than open time, then β 1 , β 2 , and β 3 are selected carefully where β 1  β 2  β 3 and assure followed inequality is true : Term (17) will assure hirarchical ordering of solving iteratively by DE.

Solving the Model Using DE
In this paper we applied this model on a real case taken from MNA (Merpati Nusantara Airlines), Ltd., a private airline company based in Indonesia [22,23].To solve the above problem we pursue the following steps of DE.

Initialization
Initial solution x ijk is defined by generating n_pop binary random numbers (0 or 1) of dimension n k × m k .Let X np,0 be the initial solution population np, then ,0 Rand np [0;1] is binary random 0 or 1 of population np.

Mutation
Different from those which usually used as a mutation operator in DE as indicated in Equation ( 2), this paper introduce a random swap as a mutation operator.Let r 0 be random number between 0 and 1 with n k × m k dimension for each population np, and v np,r0,g be element of solution V at column r 0 at generation g.If W np,g-1 is the best population for generation g-1 and w np,r0,g-1 is the element at column r 0 of W np,g-1 , then mutant of generation g is defined as , 0, 1 0 , 0, , where c m is mutation probability (0 < c m <1) which represents mutation power imposed to the best population of previous generation.Term W np,r0,g-1 mod( 2) means what value left after dividing W np,r0,g-1 by 2. The selection of c m must be done carefully because too small c m can cause old solution is difficult to exit from local optimum.While, too big c m causes noise solution such that fast convergence toward global optima can not be achieved.Selecting c m accurately becomes the successful key of implementing this algorithm.
As an illustration how this mutation operator works, notice the following example.Suppose we have 3 pilots and 4 pairings, use c m = 0.2.
Suppose we have W 0 (current solution): Entry (i,j) = 1 indicates a pilot i be placed in pairing j, entry(i,j) = 0, otherwise.After generating randomly, suppose we got: 0 0,10 0, 40 0, 08 0, 70 0,30 0, 20 0,90 0,15 0,85 0, 05 0, 25 0,55 To have a new solution, apply Equation (19).By comparing r 0 and c m for each corresponding entry in W 0 and r 0 , we have matrix W' 0 , where each entry (typed in bold) should be changed since r 0 < c m .0 1 0 0 0 The changes should be made as follows The other entries of W 0 remains the same as r o ≥ c m ,.Therefore, we obtain the new solution 0 0 1 0 We follow this process each time the algorithm is on the mutation step.

Crossover
Crossover changes over parent solution X np,g by mutant solution V np,g to construct a new solution U np,g .Crossover is done by defining threshold probability (0 < c r < 1) for mutant to change current solution.Then we generate n_pop random numbers (0,1).If the random number < c r , then the mutant replaces the current solution and otherwise.

Selection
This process is done by comparing the fitnesses of the parent solution and the new solution which is produced from crossover process.The new population will replace the old population only if the new population has better fitness than that of the old population.The fitness function refers to Equation (16).Then, the solution of the next generation X np,g+1 can be obtained from this formula: , where f(U np,g ) and f(X np,g ) are the fitness of U np,g and X np,g respectively.The constructed mathematical model consists of some aspirations/criterias as the objective function and some constraints.The objective function includes minimum flying times, deviation of flying days, and open time.Some constraints which are considered when constructing a roster include overlap, crew requirement of pairing, free days before seven days of flying days, maximum flying times, and maximum numbers of take off.

Experiments and Analysis
We used Matlab to implement DE algorithm.The experiments were done using the datasets shown in Table 1.The datasets consist of pairings, numbers of crews, type of fleed and rules.Datasets are divided into two categories: small and large datasets.Small dataset consists of assignment of F-100, CN-235, DHC-6, and Cassa 212 pilot.While, large dataset consists of assignment of Boeing 737-200 pilot and stewardess.
The experiments aim to assign crew members in which pairings.The results of the experiments on small datasets were compared with column generation [22,23] and MOSI (method used by the company).For large datasets we included exact decomposition as a comparative method [23].
Probability of mutation (c m ) and probability of crosser (c r ) are the key succes factors for DE implementation.Therefore, these two parameters should be determined carefully.The precise parameter values will lead to the global optimal solution.In this paper, these parameters settings were done through trial process.The best c m , is determined by using one dataset (Cassa 212) with npop = 50, gmax = 50, c r = 0.7 and the experiments were done with 5 replications.Table 2 shows the objective function values for different c m .The best value is typed in bold.We found that c m = 0.05 is the best one.
Using c m = 0.05, the same procedure was done to find the best c r .The results are shown in Table 3.
Next, using combination of c m and c r we tried to determine the best values of these two parameters.Through experiments, we obtained the best combination between c m and c r are 0.1 and 0.5 as shown in bold in Table 4.
Criteria in the objective function are weighted based on their significances.In this paper we put roster cost minimization as the most significance followed by deviation of minimum flying days and open time with weights of 10 4 , 10 2 and 1 respectively.Other parameter needs to set up is pinalty for the constraints.Penalties for overlap and rotation without free day, number of pilot requirements and day off constraints are 10 15 , 10 13 and 10 11 respectively based on their signifances.While, the penalties for flying day, flying times, and numbers of take off are set to 10 6 as these constraints are assumed to be the same important.
Experimental results for small datasets are indicated in Table 5.For the flying time criterion, DE spent more Table 6 shows the total deviation of average flying days for three methods.We see that DE is superior for all assignments except for Cassa 212.This proves that DE produced good roster quality.From Table 6, we see that DE produced better deviation of flying days than other methods for three assignments and MOSI is superior for Cassa 212.For the open time criteria, as shown in Table 7, DE produced superior results compared to the other methods for Cassa 212 and the same results for other assignments.
Dataset of pilots and stewardess assignment on Boeing 737-200 aircraft has a larger size.This dataset consists of 17 pilots with 82 pairings and 55 stewardess with 114 pairings.It required high number of iterations and computational time to obtain near optimal solution if random initial solutions were used.We proposed a sequence of partial optimization and total optimization techniques.In partial optimization, the dataset is splitted into several smaller sets and the corresponding optimization problems were solved by DE method.At the end of this stage, all of solutions were combined all together to form initial solution for the total optimization problem.Setting the initial solution through this process decreased the computational time significantly.
For the pilot assignment of Boeing 737-200, number of pairings is much higher than the available pilots. Therefore the resulting flight schedules violate several constraints.The results for this assignment are shown in Table 8.The table presents the flying days and the actual difference produced by each pilot.In this case, we compared 6 methods namely differential evolution (DE), column generation, MOSI, exact decomposition with 21 flying days (DEC21), exact decomposition with 20 flying days (DEC20) and exact decomposition with 19 flying days (DEC19).The results are presented in Table 8.

Conclusions
We have investigated the use of DE in solving airline The best mutation and crossover probability are 0.1 and 0.5 respectively.Different from using DE in general, in this paper we introduced random swap as mutation operator.For small datasets, DE was proven more competitive then column generation and MOSI, based on constraints fulfillment and roster quality.In the assignment of Boeing 737 pilot, DE produced smoother flying days and the least pilots which violated flying days constraint compared to other methods.For the stewardess assignment of Boeing 737-200, DE violated the least constraints compared to column generation and MOSI.DE produced superior deviation of flying days criterion but it is still inferiror for iwo other criteria.

Tabel 9 .
Violation of flying days constraint for assignment of Boeing 737-200 stewardess.scheduling problem.Generally, the rostering problem in MNA has characteristic indifferent from other airline companies in terms of roster constraint and roster quality.Selecting mutation and crossover probability accurately became the successful key to implement DE.

Table 8 . Comparison flying days and differences.
We see from the table that all of the methods violated flying days constraint.DE assigned smoother flying days than other methods based on the standard deviation values.DE reached the smallest standard deviation.We also recorded that DE produced the smallest number of pilots violating flying days constraint, that is 3 pilots.While, column generation, MOSI, DEC21, DEC20 and DEC19 respectively produced 6, 7, 7, 8 and 8 pilots violating the constraint.Different from those of pilot assignments, in stewardess assignment of Boeing 737-200, we compared DE, column generation and MOSI.Table9shows that DE only violated crew requirement constraint while column generation and MOSI violated both crew requirement and flying days constraints.We can see in Table9that column generation and MOSI assigned 6 and 7 stewardesses exceeding the maximum flying days, 21 days.While, in the pilot assignment of Boeing 737-200, as shown in

Table 10 ,
DE is superior for minimum deviation of flying days.In the assignment of Boeing 737-200 stewardess, DE is superior by its minimum deviation of flying days compared to other methods.But, DE is inferior for flying time and open time criteria.