Multiple objective solution approaches for aircraft rerouting under the disruption of multi-aircraft

Expert Systems With Applications 83 (2017) 283–299

Contents lists available at ScienceDirect

Expert Systems With Applications journal homepage: www.elsevier.com/locate/eswa

Multiple objective solution approaches for aircraft rerouting under the disruption of multi-aircraft Yuzhen Hu∗, Hong Liao, Song Zhang, Yan Song School of Economics and Management, Harbin Engineering University, 150001 Heilongjiang, China

a r t i c l e

i n f o

Article history: Received 23 September 2016 Revised 4 April 2017 Accepted 5 April 2017 Available online 18 April 2017 Keywords: Aircraft rerouting Multiple objective Polynomial-time algorithm ε -constraints method Neighborhood search algorithm

a b s t r a c t This paper considers a multi-objective aircraft recovery problem for airline disruption. An integer programming formulation is first established based on connection network with three conflicting objectives, where the first objective minimizes the total deviation from original flight schedules, the second objective minimizes the maximal flight delay time, and the third objective minimizes the number of aircraft actually attended in swapping. Optimal analysis is provided for a small scale aircraft recovery problem with the last two objectives and then one polynomial time algorithm for aircraft recovery problem after the disruption to multi-aircraft in a fleet at an airport. One heuristic combined ε -constraints method and neighborhood search algorithm is designed for large scale disruption recovery problem. The results obtained from computational experiment illustrate effectiveness and efficiency of the heuristic. The outcome of this research could provide theoretical and practical supports for airlines to reduce flight delays. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction The airline industry and academia have got great theoretical and practical advances at the planning level, such as flight, fleet assignment, aircraft scheduling and crew scheduling. The application of optimization-based approach for airline scheduling has been proved to be a great achievement. Nevertheless, actual daily operation may fail due to quite a few disturbances such as severe weather or mechanical issues. These failures are defined as disruption. If these disruptions are not mitigated appropriately and promptly, they will affect crew connections and passenger itineraries and result in significant damage to airline image and a substantial loss of profit. The process of monitoring and rescheduling relative resources to recover from the failure of the daily operation is denoted as disruption management. Given a disruption management process to the existing schedule, the airline is said to be in a recovery operation. Although disruption management in the airline industry has been a topic of concern for an increasing number of researchers; for example, see Ball, Barnhart, Nemhauser and Odoni (2007), Clausen, Larsen, Larsen, and Rezanova (2010), Etschmaier and Mathaisel (1985), Filar, Manyem, and White (2001), Kohl, Larsen, Larsen, Ross, and Tiourine (2007), and Castro and

∗

Corresponding author. E-mail addresses: [email protected] (Y. Hu), [email protected] (H. Liao), [email protected] (S. Zhang), [email protected] (Y. Song). http://dx.doi.org/10.1016/j.eswa.2017.04.031 0957-4174/© 2017 Elsevier Ltd. All rights reserved.

Oliveira (2011), there has been relatively little achievement on the disruption recovery operations. In spite of disruption recovery problems of the operational phase similar to flight scheduling problem of the planning phase, the former’s problems are aggravated by four things as follows. First it is more complicated operational factors that operations control center (OCC) needs to consider in real time environment. For example, when disruption occurs, operations control center (OCC) coordinators will typically recover disrupted flights taking flight schedules, aircraft routings and airport curfew limitation, especially the disruption reasons and disturbance scale into account. If the disruption level is not severe, then the dispatcher will reschedule flights only by swapping aircraft. If severe disruption, such as airport closure or severe weather, occurs, OCC has to cancel some flights to promise completion of total airline production. The second problem is that the staffs of OCC often wish to make decisions in as close to a real-time environment as possible. For example, flight waves are often used to increase convergence opportunities between adjacent flights especially in hub airports. Faced with small scale disruption for some flight departure, for instance, OCC will try their best recover the flights according to real-time situation before next peak hours of flight departures. If there are still some delays, then flight rescheduling willcontinue via flight waves. The third problem is that the daily operation in practice is derived from comprehensive consideration of multiple objectives instead of single objective like total cost of airline operations that is commonly used for planning environment. For example, like other

284

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

profit-making organizations, airline always seek maximum benefits in the planning level. In the daily operation, however, airline cannot ignore one phenomenon that passengers’ satisfactions with airline symbolize future potential benefits of the airline. If single objective of total delay cost is only considered to recover from one disruption, extreme delays in certain flights may occur. It can easily lead to passengers’ dissatisfaction of some flights on the airline with extremely long delays. It even can lead to some advert event such as passengers damaging airport facilities and refusing to board the plane. Therefore, extremely long delay time of flights is often avoided by OCC as much as possible. Additionally, OCC also wishes to recover flightss from disruption by swapping fewer aircraft, because it means that the flight schedules originally implemented by the aircraft may be indirectly affected if one more available aircraft is attended in the recovery process. As a result, unlike the flight planning phase, airlines do not generally incline to deal with total flight information by traditional mathematical programming method in the presence of a disruption to airline daily operations. To this day, there have been considerable achievements made in the use of mathematical programming to airline planning and the operation of flight schedules, while little advancements have been made practically in disruption management of airline operations. One possible explanation is that the model formulation with single objective has similar characters with the network flow model for airline planning phase and can be easily solved by corresponding optimization techniques. Most researchers are ones who have studied airline flight scheduling problems, and they have been accustomed to using mathematical programming methods similar with the approach on flight scheduling. However, they ignored the complexities of real time environment, stage of flight waves during daily operations, the diversity of airline objectives and the control of aircraft number in the recovery operations. Such a solution scheme may not be implementable by an OCC based on the airline’s operation situation, for it prevents choice diversity of OCC coordinators depending on airline’s practice and implementation of flexible recovery options in real time environment. The principle goal of this paper is to define, formulate, and solve an aircraft recovery problem with multiple objectives (MOARP) that is amenable to the constraints imposed by an OCC coordinator. We formally describe the process of flight rescheduling depending on the disruption scheme and present an integer programming model formulation with three objectives, which includes minimizing deviation from the original schedules, minimizing extreme flight delays, and miniziming the number of aircraft duty swap. In the context of solving the problem we also introduce some theoretical analysis for other two objectives except for the objective of total deviation from original flight schedules. Solution approaches are then developed to obtain Pareto front near optimal solution combined with ε -constraints methods for multi-objectives computation and neighborhood search algorithms for improvement of initial recovery solutions. Finally, some simulated disturbance experiments are made to validate the proposed approach. The remainder of the paper is organized as follows. After a review of relevant work done within irregular airline operations in the following content of Section 1, the problem and model are formally defined in Section 2. Section 3 discusses the theoretical result of small disruption recovery problem based on flight waves. Heuristic is outlined for aircraft recovery problem with multiple objectives in Section 4. Computational results are shown in Section 5 that validates our approach. Here we observe the flexibility of multi-objective approach for choice of recovery solution in practical application. Section 6 summarizes our work and suggests related future work.

1.1. Literature review 1.1.1. Network models Clausen et al. (2010) report that although substantial achievements have been made in developing solution methods that support the stand-alone recovery of aircraft and crew since the mid 1980 s (Teodorovic and Guberinic, 1984 & Gershkoff, 1987), the majority of the mathematical models and solution methods for solving the airline recovery problems are similar to the methods applied for planning purposes. Tools for planning as well as for recovery are, in most research cases, based on a network flow model that describes how flights can be sequenced in a rotation. There are three types of networks: Time space network, time-band network, and connection network. Time space network and Time band network are applied widely between the end of 20th century and the beginning of 21th century. Yan and Yang (1996) developed a time space network model formulation to deal with flight rescheduling due to the breakdown of aircraft with the objective of minimizing flight total cost of delays and cancellation, and solved it by combining network simplex algorithm and sub-gradient method with Lagrangian relaxation. Their work was extended by Yan and Lin (1997) to address airport closures, by Yan and Tu (1997) to address multiple fleet substitutions, by Thengvall, Bard, and Yu (20 0 0) to reduce passenger dissatisfaction using protection arcs and by Thengvall, Yu, and Bard (2001) to integrate all recovery options in a model. Zhang, Yu, Desai, and Lau (2016) represented a two stage iterative heuristic algorithm to solve the integrated aircraft and crew schedule recovery problem, and the model in each stage is constructed by the multicommodity network flow model for aircraft rescheduling and crew rescheduling respectively. Bard, Yu, and Argüello (2001) represented a time-band network model to solve aircraft recovery problem. Bratu and Barnhart (2006) described two models, abbreviated DPM and PDM, based on a time-band network. The test data set is limited into 4 aircraft. Three different scenarios with different levels of disruptions were analyzed only with DPM since PDM was not suitable for operational use. Eggenberg, Bierlaire, and Salani (2010) solved the similar problem using column generation algorithm, where each aircraft route is computed by time band network model. The model was then extended by Hu, Xu, Bard, Chi, and Gao (2015, 2016) to solve integrated recovery problem of aircraft and passengers. Another important achievement of airline recovery research is the construction of set partitioning model formulation based on connection network. Rosenberger, Johnson, and Nemhauser (2003) added additional time slot and capacity constraints into the connection network model and then apply column generation algorithm and CPLEX 6.0 to solve the aircraft recovery problem. Andersson and Varbrand (2004) solved the similar problem by Lagrangian relaxation-based heuristic and Dantzig-Wolfe decomposition technique. All above literatures only consist on single objective programming to minimize total recovery operation cost instead of on multi-objective programming to avoid extreme cases and to reduce the aircraft swapping scale. The methodologies we develop in this work address precisely these limitations. 1.1.2. Solution approaches In addition to similarity of network model formulation, researchers apply solution methods similar to the state-of-the-art within solution approaches for the tactical planning problems. Clausen et al. (2010) concluded that so far no planning tools have been able to cope with the complexity of re-planning all airline operations at the same time during disruptions. Most literatures try to reduce the problem space by limiting the recovery time period and the number of crew members or aircraft included into recovery.

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

Jafari and Zegordi (2010) constructed a model for integrated recovery of aircraft and passengers based on aircraft rotations with the objective of minimizing the total cost of recovering all flights, aircraft, and passengers. But the data set in this paper only contained 13 aircraft divided into 2 fleets. Bisaillon, Cordeau, Laporte, and Pasin (2011) presented a neighborhood search heuristic (LNS) to recover aircraft and passengers simultaneously with the objective of minimizing operating costs and impacts on passengers. The core components of the heuristic include a construction phase, a repair phrase and an improvement phase. In each iteration, aircraft routes optimization is obtained first and then passenger itineraries are reassigned based on the aircraft routes. Sinclair, Cordeau, and Laporte (2014) improved LNS heuristic via adding some additional steps in each phase. Then Sinclair, Cordeau, and Laporte (2016) presented a column generation post-optimization heuristic based on LNS heuristic in Sinclair et al. (2014) to obtain improved recovery solutions within reasonable computing time. The approach was verified by some instances for different limitations of recovery period, number of aircraft or disruption situations, such as flight disruption, aircraft disruption and airport disruptions. Maher (2015b) developed column-and-row generation framework based on general column generation approach and Bender’s Decomposition to improve the solution runtime and quality. Maher (2015a) applied the column-and-row generation method for integrated recovery problem with passenger reallocations. It is limited in the computation result solutions by 6 hours, 8 hours and a single day of recovery period, and 48 aircraft in a single fleet type. Although each one of the above solution approaches can give one feasible or near optimal solution through adding the constraints of recovery period length and the constraints of the number of available resources, it is still not sufficient for the OCC staff to choose one recovery solution that he/she wants according to real time environment. However, our solution methods can provide a set of recovery solutions, which meets the selection requirements of the OCC staff, by moving the corresponding constraints to the objective function and computing Pareto front of the multiobjective programming.

1.1.3. Flight waves Although OCC coordinators stated that they often reschedule flights by flight waves in actual situations, the principle was used in few academic researches. Recovery stage is defined in Abdelgahny, Ekollu, Narisimhan, and Abdelgahny (2004), where the flights in each set are resources independent, i.e. two flights in each stage cannot be covered by the same aircraft or crew members, while one aircraft implements the flight in earlier stage can also serve the flight in later stage. The definition of recovery stage is similar with that of flight waves, especially in hub airports. Brute and Barhart (2006) reported that in their simulation the daily operation was divided into different simulated system states. Recovery model was populated with data reflecting the current simulated system state, and then solved to determine new departure times and cancellations for the remaining flight legs in current day. Although solutions were generated for the rest daily operations, it is the decisions in the current system stage that were only implemented. In the subsequent stage, the simulated airline system state was then updated, recovery model was resolved, solutions were generated for the rest of the day, and associated decisions are implemented for this stage. The process repeated until the implements of last-stage decisions had been finished. However, the literatures only described the idea of recovery stage, lack of application of the stage in the detailed analysis and proof of optimization results.

285

1.1.4. Multi-objective As stated above, most studies for traditional airline recovery problems only solve single objective problems. Even faced with multi-objective model, they often propagate a linear-sum on objectives to change them into the single one. However, It is stated that the single objective solution cannot overall reflect the true cases with practical multi-objective needs by Fieldsend and Singh (2005). Multi-objective solution approaches have been successfully applied in robust airline scheduling problems (Burke, Causmaecker, Maere, Mulder, Paelinck, & Berghe, 2010; Chou, Liu, & Lee, 2008; Lee, Lee, & Tan, 2007). But, there are very few literatures to study the multi-objective optimization problems of the aircraft schedule recovery. Liu, Chen, and Chou (2010) constructed a multiobjective combinational optimization model formulation for daily short-haul recovery problems, and developed a hybrid evolutionary algorithm composed of an adaptive evaluated vector and the method of inequality-based multi-objective genetic algorithm. A simulated disturbance experiment of daily domestic airline plans in Taiwan with only 7 aircraft and 39 flights was tested. However, the objectives are lack of formal description of flight cancellation, which is a common phenomenon in most airlines. The short haul recovery approach is not suitable for general scale disruption in bigger airlines. 1.1.5. Robust aircraft planning In addition to the reactive ways of handling disruption as stated above, another proactive method to deal with airline disturbances is robust airline planning, which makes flight, crew schedules and aircraft routings less sensitive to disruptions. It can either absorb disruption completely or make an easy recovery from disruption using different types of recovery approaches. In recent years, many publications on robust scheduling planning appear in the literature and one of the classical topics is robust aircraft routings. Lan, Clarke, and Barnhart (2006) first presented two independent models to reduce passenger missed connection and achieve robust airline schedules. Then Dunbar, Froyland, and Wu (2012, 2014) extended the research by integrating optimizing aircraft routings and crew pairing into one model formulation and analyzing how these dependencies affect the propagation of delays through the flight network and minimizing the expected total propagated delay. While Yan and Kung (2016) established a new model to minimize the maximal possible total propagated delay and evaluated it based on the benchmark approach (Dunbar, Froyland, & Wu, 2014). It should be mentioned that robust airline planning is a difficult task, and it cannot always keep the plan feasible particularly due to the unpredictability of disturbances. Sometimes, in case of major disruption, other good recovery measures are also required in order to restore the planning operations. 1.2. Contributions The main contributions of our research can fall broadly into three categories as follows. (1) The integer programming model constructed in this study reflects many practical objectives, such as minimizing total cost of flight delays and cancellations, minimizing the maximal flight delay time and minimizing the number of swapping aircraft, which are usually involved by aircraft recovery problems. (2) It is the first time that a recovery solution approach for resolving multi-objective airline rescheduling problems is developed combined ε -constraints method and neighborhood search algorithm, where the ε -constraints method is in charge of seeking the Pareto front for ARPMO and the

286

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

neighborhood search algorithm is responsible for improving the local feasible solutions of ARP in each iteration of the ε constraints method. This approach is handier for the choice of OCC operator in real time situations. (3) This paper firstly analyzes the optimality of MOARP with the objectives of minimizing the maximal flight delay time and minimizing the number of swapping aircraft according to a common flight rescheduling options of flight waves, and designs one polynomial algorithm for it. The optimal analysis give a major concussion for the expression that airline rescheduling problem is NP- hard, although traditional airline rescheduling problem with minimizing the total recovery cost has been proved NP- hard.

2. Problem statements 2.1. Disruption scenarios Given a set of flights F originally covered by a set of aircraft P, if a subset of aircraft P’ ⊆ P is disrupted during the time horizon (Tbegd , Tendd ), where Tbegd denotes the beginning time of the disruption and Tendd denotes the end time of the disruption, then the flights schedules implemented by P’ will be failure. OCC operators need to take a series of recovery measures to get back to the initial schedule as soon as possible based on information of a planned (or initial) schedule and disruption state, i.e. the location and available time of the planes at the moment the disruption occurs. According to the number of aircraft in the set P’ and the difference value between Tendd and Tbegd , disruption scale can be generally divided into several types. For example, based on regulation of one big airline in China, if more than 10 aircrafts in one airport are disrupted in successive 2 hours, then it will be defined as large scale flight delays. In this paper, we do not necessarily give an exact division of different scales of the airline disruption for our research, instead we give a principle as follows. If flight delays were only permitted for disruption recovery, airport curfew time can still be satisfied, the disruption is defined as small scale disruption. Otherwise, it is called as large scale disruptions.

2.2. Recovery period Recovery from disruption cannot be done at one go. It needs time horizon to recover the initial flight schedules. We refer to it as the recovery period denoted by T = (Tbegr ,Tendr ), where Tbegr represents the beginning time of the recovery horizon and Tendr represents the end time of the recovery horizon. The time horizon is longer than the disruption horizon, i.e. Tbegr ≤ Tbegd and Tendr > Tendd . The length of the recovery period is not deterministic. The recovery horizon lasts until all changes caused by the disruption have been carried out. Recovery horizon starts at the same time or earlier than the departure time of the first updated flight and ends with or later than the new arrival time of the last updated flight. Generally speaking, the end time of the recovery period should be controlled in the current day to avoid the regular airline operation of the next day. Therefore, it is defined that Tendr is earlier than the airport curfew time in this study no matter the disruption scale. For small scale disruption, the recovery horizon may ends in advance, so the recovery techniques should be altered according to the flight delays situations. For large scale disruption, some measures should also be taken to promise the finish of the recovery horizon by the end of the current day.

2.3. Recovery options When a disruption occurs, some recovery options are considered generally: swapping aircraft (exchanging aircrafts that have few assignments with grounded aircrafts of the same type or different type), delays (providing new departure times for the flights according to the available times of the delayed aircrafts or the available times of newly assigned aircrafts), cancellations (cancelling flights and reassigning their passengers), and calling reserved aircraft (using spare aircrafts) and so on. During the recovery period, the previously mentioned recovery options are utilized according to real-time environment. As flights can be finished before curfew time even though there is no swapping between aircraft, aircraft swapping and flight delays are sufficient for recovery from a small scale disruption. What the OCC operators concerned about is only how to reduce the delay time of the flights by aircraft swapping in airports successfully. After all other recovery options, such as flight cancellations and calling reserved aircraft, are cost-wasting relatively. Faced with large scale disruption, however, the measures only including the two recovery options may lead to violation of airport curfew time limitation. Therefore, the above four recovery options are all applied for flight rescheduling from a large scale disruption to avoid the extension of the disruption into the following day. In addition, aircraft swapping only within single fleet will be considered for small scale disruption, and crew availability will not affect the generation of alternate routings. However, large scale disruption might need to aircraft swapping between different fleet types for recovery as soon as possible. In practice, aircraft swapping between the same aircraft fleet type has higher priority over different fleet types, as each aircraft type has different characteristics that limit cross-assignments. In particular, different types have different seating capacities. Substitution is allowed between some types in accordance with airline policies and aircraft configurations. Intuitively, aircraft with larger capacities can substitute for smaller ones but not vice versa. 2.4. Flight waves Flight waves are generated by alternation of peak-valleys of flight departures and peak-valley of flight arrivals. In a time horizon a large number of flights arrive and a few number of flights leave from the airport, it is defined as arrival waves. While in next time horizon a few number of flights arrive and a large number of flights leave from the airport, it is then defined as departure waves. If some incoming flights are delayed due to some disruptions at this airport, then OCC operators can reroute the aircraft at the airport for significantly reducing departure flight delays, applying the flight connection between arrival waves and departure waves. The formation of flight waves can increase the possibility of swapping aircraft in a relative short time horizon and promotes the recovery from the disruption especially in some big airlines. Daily operation can be divided into several stages by flight waves. The recovery decision is implemented successively. In the current stage, information of disruptions scheme of flights and available resources is collected and recovery solution for the rest daily operation is generated. Then in next stage, all information is collected again and new recovery approach is applied for the rest of the whole day. The update of recovery stages is continued until no flight delay occurs in the following flight operation of the rest of the day before the end of the recovery period. 2.5. Connection networks model In this study, the aircraft recovery problem will be modeled as to allocate limited resources (aircraft) P to activities (flights) F with

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

multiple objectives. This model is similar with the set partitioning problem with additional side constraints, while the difference mainly lies in the multiple objectives of the model. A route is defined as a series of flights connected in time and in space that can be flown by a single aircraft. An aircraft route refers to the matching of one route with one aircraft which can fulfill curfew time limitation of the airport where the last flight in the aircraft route arrives. Note that an aircraft route may include the original scheduled flights, as well as possible delayed flights with specified delay times. As each aircraft has its own available time, the delay time of the flight in the same route but different aircraft routes may be totally different. Therefore, for each flight, we can estimate its delay cost based on the aircraft route it belongs to, and can also evaluate its cancellation cost by average ticket price of the flight.

2.5.1. Objectives of the problem The problem investigated in this paper reflects the practical operation of an airline when irregularities force its operating schedule to be altered. When some aircraft are delayed or grounded due to unforeseen events, the flights covered by those flights will be disrupted. To recover the schedule as soon as possible, available aircraft can be rerouted, delays can be imposed, and surplus aircraft can be called into service. Several objectives have been proposed to guide the recovery, depending on the data available at the time of the disruption and the priorities of the airline. In our paper, the objectives include the following three aspects. (i) Minimization of the total deviation from the original flight schedules: The deviation is composed of flight delays and flight cancellations, which are both unified by total cost based on the deviation. It is a conventional objective in most literatures and can generally dominate the realization of airline disruption management. (ii) Minimization of the maximal flight delay time during recovery. It is important for airline to provide a reliable service to passengers in order to attract high-value passengers who are sensitive to on-time reliability airline, to raise customer loyalty. However, recovery solution based on the single objective of minimizing total deviation may lead to extreme long time delay of some flights. The long delay times might not only cause extra cost to airlines by transferring passengers to other airlines or requiring the provision of meals and other services, but also reduce customer retention rates. To fulfill the promised service level, minimizing the maximum flight delay time objective can neutralize the extreme phenomenon derived from the first objective. (iii) Minimization of the number of swapping aircraft applied during recovery. It ensures that the original flight schedule is assigned to the same aircraft as much as possible. This objective reduces the additional cost of the non-profit duty swapping for flights and the inconvenience of crew change. Considering the above three objectives of the MOARP problem, the objectives are normally combined into a single scalar function in conventional approaches. However, a set of appropriate weights to combine the objectives is difficult to find. Furthermore, because each objective has a different definition of optimality, dimensional disunity of the three objectives impedes the progress of consolidation in the three objectives by simple weighted sum methods. This study formulates the MOARP problem as a multi-objective network model in this section and then resolved it in next section applying the well-known ε -constraint technique to obtain a set of non-dominated, alternative airline recovery solutions, known as the Pareto-optimal set.

287

2.5.2. Constraints of the problem As the MOARP problem that this paper addresses closely reflects the problem that must be solved in practice, besides considering the above three objectives, each aircraft route in the connection network must be feasible with respect to the following flow balance and operational constraints as well: (i) flights for each feasible aircraft route can only be adjusted once or cancelled; (ii) each aircraft can cover at most one route to form one feasible aircraft route; (iii) each flight in a new schedule should not depart earlier than its originally scheduled departure time prior to the disruption; (iv) subsequent departures for each aircraft should be scheduled only after the previous flight arrives and only after a minimum turnaround time, that is, aircraft flow should be continuous; and (v) the airport arrival curfew time should be observed. (vi) no aircraft scheduled for maintenance service during the recovery period will have its original route altered. (vii) at the end of the recovery period, enough number of appropriate aircraft types should be positioned at each airport in the network to ensure that the published schedule can be executed. If this is not the case, the disruption will extend into the following day – a situation that we wish to avoid.

2.5.3. Mathematical formulation In this subsection, a multiple objective integer programming formulation, based on connection network model, is constructed to describe the MOARP problem more accurately. In this formulation, Aircraft will be assigned to the routes and the match of aircraft and routes must fulfill constraints (iii), (iv), (v) and (vi). Constraints (i), (ii) and (vii) are enforced explicitly in the following formulation. The following notations are used in the formulation. Index I Aircraft index e Fleet index F Flight index r Route index s Airport station index Sets P E F R S P(e) R(f) R(s)

Set of disrupted aircraft and available aircraft for swapping in the recovery period (Tbegr ,Tendr ) Set of fleet types for the set of aircraft P Set of flights originally covered by aircraft set P Set of routes which can cover by single aircraft in P Set of airport stations where all flights in F fly through Set of aircraft that belong to fleet e; e∈E Set of routes covering flight f; f ∈F Set of routes end in the airport station s; s ∈S

Parameters drfi Delay cost of flight f covered by route r and aircraft i; r ∈R, f ∈F, i ∈P yri parameter that represents whether route r is originally covered by aircraft i: = 0 means route r is originally covered by aircraft i, = 1 otherwise; r ∈R, i ∈P cf average cost of cancelling flight f; f ∈F trfp Delay time of flight f covered by route r and aircraft i; r ∈R(f), f ∈F, p ∈P mse Number of aircraft belonging to fleet e that should terminate at airport station s; s ∈S, e ∈E Decision variables xri = 1 if route r to which aircraft i that is assigned is considered in the solution, = 0 otherwise; r∈R, i ∈P The first objective of the multiple objective aircraft recovery problem is illustrated as follows.

288

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

min

 

Z1 =

f ∈F i∈P r∈R( f )

dr f i xri +

 f ∈F

 cf 1 −

  i∈P r∈R( f )

 xri

(1)

Where the first term shows the flight delay cost and the second term gives the flight cancellation cost. If multipliers of xri are merged, then we can get the following formulation.

min



Z1 =

f ∈F

cf −

 

f ∈F i∈P r∈R( f )

(c f − dr f i )xri

Definition of swapping stage

one flight circle of available aircraft corresponds to one swapping stage at most.

Lemma 1

The optimal solution of SP1 ca be got in swapping stage 1.

Lemma 2

IN5 and IN6 is not useful for obtaining the optimal solution.

Lemma 3 and Lemma 4

(2)

The second and the third objectives in the formulation are given in (3) and (4), where Z2 refers to the maximal flight delay time in a feasible recovery solution and Z3 the number of aircraft that is assigned to the routes not representing the original flight schedules covered by the aircraft.

min

min

Z2 = max f

Z3 =

 i∈P r∈R

 

tr f i xri

(3)

i∈P r∈R( f )

yri xri

(4) BSMCF is designed for SP1 and SP2

As stated above, the multiple objective aircraft recovery problem is subject to the following detailed constraints.

 

∀f ∈ F

xri ≤ 1

i∈P r∈R( f )



∀i ∈ P

xri ≤ 1

(6)

r∈R

 

xri ≤ mse

∀s ∈ S,

e∈E

(7)

i∈P (e ) r∈R(s )

xri = {0,

1}

∀i ∈ P,

Fig. 1. Roadmap of the proofs for small scale problem.

(5)

r∈R

(8)

Constraint (5) promises that each flight is either cancelled or reassigned to one aircraft route. Constraint (6) ensures that each aircraft can at most be assigned to one aircraft route. Constraint (7) stipulates the balance of aircraft at each airport station by the end of the recovery period so that flight schedules for the next period/day can then begin as scheduled. Constraint (8) defines the decision variables as being integer or binary. MOARP can be modeled accurately with the formulation listed in Eqs (2)–(8), however, it is computationally very difficult to solve. The difficultly results from the fact that, for each aircraft i ∈P, all possible routes that it can be covered by i must be generated. Note that the number of feasible routes may be very large—it grows exponentially with the number of aircraft for swapping and time length between Tbegr and Tendr planning time horizon. In an attempt to generate feasible recovery flight schedules for aircraft reroutings guided by the three objectives, we develop a heuristic based on neighborhood search algorithm and ε -constraints method to address this problem. Before giving the practical heuristic, we first give some optimal analysis for the MOARP problem under small scale disruption guided by the last two objectives: minimization of the maximal flight delay time and minimization of the number of aircraft participating in swapping for airline recovery. The two objectives will be solved hierarchically in next section. 3. Recovery from small scale disruption of aircraft within one airport As we all know that traditional aircraft recovery problem with the objective of minimizing total deviation from the original schedules has been proved to be NP-hard (Hu et al., 2015), therefore, it is impossible to obtain the optimal solutions by polynomial algorithms. However, some optimal analysis for other two objectives will be given in this section based on flight waves. The corresponding roadmap is given in Fig 1.

First, Hierarchical goal programming model (SP1 and SP2 ) are established, and swapping stage and swapping intervals are defined under small scale disruption and recovery option of aircraft swapping in Section 3.1. In Section 3.2, the set of swapping intervals of available aircraft for swapping are divided into several subsets in the first swapping stage and then Lemma 1 try to prove that one flight circle of available aircraft corresponds to one swapping stage at most, i.e. there is at most one flight circle of one aircraft in a single swapping stage. Lemma 2 shows that the optimal solution of SP1 can be got in swapping stage 1. IN5 and IN6 are not useful for obtaining the optimal solution. Therefore, the interval subsets IN1 , IN2 , IN3 and IN4 are enough for obtaining the optimal solution of small scale problem. In Section 3.3, finally, the proposed algorithm BSMCF is designed for resolving the optimal solution to the problem. 3.1. Definition of small scale disruption and recovery situation Let ni denotes the number of flight originally covered by aircraft i∈P, then ri ={fi 1 ,fi 2 , …, fini } represents the route originally covered by aircraft i∈P, where {fi 1 ,fi 2 , …, fini } is defined a flight sequence originally covered by aircraft i∈P. In the flight sequence {fiu ,fi ( u +1) , …, fiv }⊆{fi 1 ,fi 2 , …, fini }, if the departure airport of fiu (denoted as dep(fiu ))is the same as the arrival airport of fiv (denoted as arr(fiv )), then the flight sequence {fiu ,fi ( u +1) , …, fiv } is named as a flight circle (denoted as c). Set tc is the duration of the flight circle c, i.e. the difference between the arrival time of the flight fiv and the departure time of the flight fiu . Set dij is the new departure time for flight fij and sij is the original scheduled departure time, and then the delay time is denoted as eij (eij = dij −sij , ∀i∈P, j = 1,2,…, ni ). What this section is hoped to solve is the model formulations of (3–8) using hierarchical sequence method, where model formulations of (3–8) is divided into two models (SP1 and SP2 ) with single objective respectively. SP1 is composed of formulations (3) and (5–8) and it can also be represented as formulation (9). Set Z2 ∗ represent the optimal value of SP1 . SP2 can be represented as one equation stated in (10).

SP1

minZ2 (X )

s.t.

X ∈D

  SP2

min Z3 (X )

s.t.

i∈P r∈R( f )

(9) tr f i xri ≤ Z2∗ X ∈D

∀f ∈ F

(10)

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

Fig. 2. Definition of swapping stage.

Where X = {xri } present the set of decision variables defined in the mathematical model formulation (2–8), and D is the feasible region defined by Constraints (5-8). For the disruption resulting from temporary delays of some aircraft (aircraft set P’⊆P), flights in the set Fi , i∈ P’ will not operate as scheduled. As flights can be finished before curfew time even though there is no swapping between aircraft under small scale disruption, then flight cancellation is not necessary. What the OCC operators concerned about is only how to reduce the delay time of the flights by aircraft swapping in airports successfully. It is then important to reschedule flights to minimize the maximal flight delay time using as fewer as possible swapping aircraft by two recovery options. (1) rerouting aircrafts; or (2) delaying the flight according to the estimated availability of the disrupted or the newly assigned aircraft. Generally speaking, a small scale disruption often occurs in a single airport instead of in several airports simultaneously, therefore, swapping aircraft must happen in this airport first and the airports where the delayed flights arrive successively via flight waves. For all aircraft in set P’ disrupted in the airport ap, a set of available aircraft (denoted as P’’ ⊆P/ P’) in the airport ap will be sought for swapping with the disrupted aircraft firstly. This is defined as swapping stage 1. After swapping stage 1, if ∃ei 2 > 0, for any i∈ P’ ∪ P’’ (It is noted that the flight sequence {fik | i∈ P’ ∪ P’’, k = 1,2,…, ni } may be not actually covered by aircraft i after swapping), some aircraft in the departure airports of fi 2 (dep(fi 2 )), will then be sought for swapping with the aircraft which covers{fik |k = 2,3,…, ni }. The swapping stages repeat until there is no delayed flight in the following flight sequence, and the last stage is denoted as swapping stage N. An example of swapping stages is shown in Fig. 2, where each node presents a feasible solution, and the number in the node represents the swapping stage. Each aircraft i∈ P’’ will stay in the disrupted airport during the time period [ai , sij ], where ai means the available time of aircraft indexed i for next flight and sij means the scheduled departure time of the next flight fij and they satisfy ai < sij . During this time period, the aircraft i may have a chance to swap with the disrupted aircraft to reduce the delay time of the following flights. The period of time is defined as swapping interval, and denoted as ini . One example of swapping intervals is presented in Fig. 3. There are three sample aircraft routes. The first route, originally assigned to aircraft indexed 1, originates in PEK with flight f11 and terminates in PEK with flight f12 . The second route, originally covered by aircraft indexed 2, begins from PEK with flight f21 and ends in TAO with flight f22 . The third route, scheduled to aircraft indexed 3 originally, begins with flight f31 from PEK to WEH, and ends with flight f33 from PEK to YNJ. Aircraft indexed 1 and 2 are assumed to be disrupted in airport PEK until time a1 and time a2 respectively. The maximal delay time of flight f11 and flight f21 is the highest if no aircraft swapping is used to handle this disruption. Thus only

289

swapping aircraft in the airport PEK can reduce the delay time of flight f11 or flight f21 effectively. Any swapping in other airports is futile. The departure airport of both flights f31 and f33 , denoted as dep31 and dep33 , is PEK, that is dep31 = dep33 = dep11 = dep21 = PEK. Therefore, there are two swapping intervals corresponding to aircraft indexed 1 and 2 in airport PEK, one time interval is the time period [a3 , s31 ] and the other one is [a3 ’, s33 ]. For simplicity, the disruption intervals of aircraft 1 and 2 are defined as [s11 , a1 ] and [s21 , a2 ] respectively. Disrupted aircraft 1(disrupted interval [s11 , a1 ]), disrupted aircraft 2 (disrupted interval [s21 , a2 ]) and aircraft 3 (swapping interval [a3 , s31 ] or swapping interval [a3 ’, s33 ]) will swap with each other in the swapping stage 1. After swapping stage 1, aircraft 1, 2 and 3 may not cover their original flight strings. If we assume the swapping result in the swapping stage 1 is that aircraft 1 covers flight string {f11 , f12 }, aircraft 2 covers flight string {f31 , f32 , f33 } and aircraft 3 covers flight string {f21 , f22 }. Flight f11 and flight f31 are still delayed and the delay times are a1 −s11 and a2 −s31 . Through delay propagation, if the scheduled departure time of flight f12 is earlier than the available time of aircraft 1 in the airport SHE, then swapping stage 2 will implement in SHE, where the following flight string covered by aircraft 1 is {f12 } now. Similarly, through delay propagation, if the scheduled departure time of flight f32 is earlier than the available time of aircraft 2 in the airport WEH, then swapping stage 2 will implement in WEH, where the following flight string covered by aircraft 2 is {f32 , f33 }now. From the above definitions, we know that swapping stage m does follow the swapping stage m − 1. Different swapping strategies in swapping stage m − 1 will lead to totally different swapping intervals in swapping stage m. To make sure the solution feasibility of the MOARP, we also make the following assumptions as well. The duration of the disruption period is less than the minimum value of the time length of all flight circles.

ma − ms < minc tc

(11)

ma = max{ai |i ∈ P ’}

(12)

ms = min{si 1 |i ∈ P ’}

(13)

Where ai is the available time of the disrupted aircraft indexed i, i∈P’, and tc refers to the duration of the flight circle c. 3.2. Characteristic analysis This section is concerned about the characteristic analysis of the problem formulated by SP1 . In swapping stage 1, the set of swapping intervals is denoted as IN, then IN is composed of six subsets according to the definition of swapping intervals and the relationship between swapping intervals and disruption period.

IN = IN1 ∪ IN2 ∪ IN3 ∪ IN4 ∪ IN5 ∪ IN6

(14)

IN1 = {ini |ai ≤ ms, ma ≤ si 1 , i ∈ P ’}

(15)

IN2 = {ini |ai ≤ ms < si 1 < ma, i ∈ P ’}

(16)

IN3 = {ini |ms < ai ≤ si 1 < ma, i ∈ P ’}

(17)

IN4 = {ini |ms < ai < ma ≤ si 1 , i ∈ P ’}

(18)

IN5 = {ini |si 1 ≤ ms, i ∈ P ’}

(19)

290

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

Fig. 3. Sample of swapping intervals.

ࢺWe have the conclusion

IN6 = {ini |ai ≥ ma, i ∈ P ’}

(20)

Z2∗ = min{ min m

Lemma 1. ∀i∈P’’, let Mi ={in | in ∈ IN1 ∪ IN2 ∪IN3 ∪IN4 }, then |Mi | ≤ 1. Proof. If ∃i∈P’’, |Mi | > 1, assuming |Mi | = 2, that is Mi ={in, in’}. Then according to the definition of swapping intervals, there must be one flight circle c, between in and in’, and it satisfies tc < ma− ms. However, it contradicts the assumption illustrated in formulation (11–13).  Let SOLm denotes the set of feasible solutions after swapping stage m, m = 1,2,…,N. Let valuesol denote the value of objective formulation (9) of the feasible solution where sol∈SOLm , then we have the following lemmas. Lemma 2. Z2∗ = minsol∈SOL1 valuesol . Proof. According to the definition of swapping stages, we easily know that

Z2∗ = min{ min m

sol∈SOLm

valuesol }

(21)

After swapping stage m, for any one sol ∈SOLm , m = 1,2,…, if ei ( m +1) > 0, then the set of swapping intervals will be sought in depi ( m +1) to carry out the swapping stage m + 1. It should be noted that the eim stay stable during swapping stage m + 1. It is assumed that the feasible solution sol’ ∈SOLm+ 1 is obtained following sol ∈SOLm, then for ∀sol ∈ SOLm+1 ,

∃eim > 0

(1) If valuesol = eim , as eim remains unchanged, then

valuesol = max{eim , maxei(m+1) } i

i

(23)

From the formulations (22) and (23) we get

val uesol = max{val uesol , maxei(m+1) } i

∴ valuesol ≥ valuesol ࢺ∀sol ∈ SOLm+1 ,valuesol ≥ minsol∈SOLm valuesol ∴ minsol∈SOLm+1 val uesol ≥ minsol∈SOLm val uesol ≥ ... ≥ minsol∈SOL1 valuesol ∴ minsol∈SOL1 val uesol = minm {minsol∈SOLm val uesol }

valuesol } = min valuesol sol∈SOL1

(24)

(25) 

∗

From Lemma 2, we know that Z2 can be obtained after swapping stage 1. From Lemma 1, we know that in swapping stage 1, the aircraft for swapping and the swapping interval follows a one-to-one mapping. The following Lemmas are all described in swapping stage 1. As aircraft swapping of swapping stage 1 in different disrupted airports is independent of each other, Z2 ∗ can be obtained by choosing the minimal value of all optimal solution of each disrupted airport. It is our major task in this paper to get the optimal value in each disrupted airport independently. In the following part of Section 3, all statements are described based on only one disrupted airport, and they are similar in other disrupted airports. The swapping interval set IN5 is excluded from the optimal solution by the proof of Lemma 3 and the swapping interval set IN6 is excluded from the optimal solution by the proof of Lemma 4. As aircraft swapping In the following discussions, we use (inl , fi 1 ) to denote that flight sequence {fij | j = 1,2,…,ni } are covered by aircraft indexed l, l,i ∈P’ ’, OPT to denote the set of optimal solutions of the model formulations SP1 , INsol represent the swapping intervals which participate in swapping actually. Lemma 3. ࢜sol ∈OPT satisfying IN5 ∩ INsol = ø and valuesol = Z2 ∗ . Proof. We assume ∃sol∈OPT satisfy

IND = IN5 ∩ INsol = {in1 , in2 , . . . , ink }, s11 ≥ s21 ≥ . . . ≥ sk 1 (22)

(2) If valuesol = ei’m , then we know that ei’m >eim . As ei’m and eim remains unchanged, so

valuesol = max{ei m , maxei(m+1) }

sol∈SOLm

(26)

and

valuesol = Z2 ∗ The above assumptions fail to hold, as we can always find one feasible solution, the value of which is not higher than valuesol . It is assumed that part of the matching in the feasible solution sol are (ink ,fu 1 ) and (inv ,fk 1 ). The delay time of flight fu 1 and fk 1 are

eu 1 = du 1 − su 1 = max{ak , su 1 } − su 1 = max{ak − su 1 , 0}

(27)

ek 1 = dk 1 − sk 1 = max{av , sk 1 } − sk 1 = max{av − sk 1 , 0}

(28)

If inu ∈IND, then su 1 ≥ sk 1 ≥ ak If inu ∈IND, as inu ∈INsol , then inu ∈ IN5 and su 1 > s01 ≥ sk 1 ≥ ak From (27) we get eu 1 = 0

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

A feasible solution sol’ is created, where ink ∈IND and (inv ,fu 1 ). Other matching is the same as that in sol. As ink does not take part in any swapping in the feasible solution sol’, the delay time of the flight fk 1 , denoted as e’ k 1 = 0. The delay time of fu 1 in sol’, denoted as e’ u 1 , is

e’u 1 = du 1 − su 1 = max{av , su 1 } − su 1 = max{av − su 1 , 0}

(29)

ࢻsu 1 ≥ sk 1 ࢺav −su 1 ≤ av −sk 1 ࢺfrom the formulations (28) and (29) we get e’u 1 ≤ ek 1 , where equality can then be used if av ≥ sk 1 = su 1 . ࢺvaluesol’ ≤valuesol A contradiction is shown compared with valuesol = Z2 ∗ , so the initial assumption must be false. Overall, the above assumptions is proved to fail.  Lemma 4. ࢜sol ∈OPT satisfying IN6 ∩ INsol = ø and valuesol = valueopt . Proof. We assume that ∃sol ∈ OPT satisfy

IND = IN6 ∩ INsol = {in1 , in2 , . . . , ink }, a1 ≤ a2 ≤ . . . ≤ ak

(30)

and

val uesol = val ueopt The above assumptions can be proved to fail, as we can always find one feasible solution the value of which is not over than valuesol . It is assumed that part of the matching in the feasible solution sol are (ink ,fu 1 ) and (inv ,fk 1 ). The delay time of flight fu and fk are

eu 1 = du 1 − su 1 = max{ak , su 1 } − su 1 = max{ak − su 1 , 0}

(31)

ek 1 = dk 1 − sk 1 = max{av , sk 1 } − sk 1 = max{av − sn 1 , 0}

(32)

If inv ∈IND, then av ≤ak ≤sk 1 If inv ∈IND, as inv ∈INsol , then inv ∈ IN6 and av < a0 ≤ ak ≤ sk 1 From the formulation (32) we get ek 1 = 0. One feasible solution sol’ is created, where ink ∈INsol and (inv , fu 1 ). Other matching is the same as that in sol. As ink do not take part in any swapping in the feasible solution sol’, the delay time of the flight fk 1 in sol’, denoted as e’k 1 = 0 The delay time of the flight fu 1 in the feasible solution sol’, denoted as e’u 1 , is

e’u 1 = du 1 − su 1 = max{av , su 1 } − su 1 = max{av − su 1 , 0}

(33)

ࢻav ≤ak ࢺav −su1 ≤ ak −su1 ࢺFrom the formulations (31) and (33) we get e’u1 ≤eu1 ࢺ A contradiction is shown compared with the assumption valuesol = Z2 ∗ , so the initial assumption must be false. Overall, the above assumptions is proved to be false.  From Lemmas 3 and 4, we know that if the swapping interval set IN1 does not exist, the optimal solution can be obtained by swapping within the swapping intervals belongs to the interval set IN2 , IN3 and IN4 . In next section, we will propose one polynomial algorithm to find the optimal solution. 3.3. The proposed algorithm In this section, one algorithm for SP1 and SP2 , called Binary search - minimum cost flow algorithm (BSMCF), is designed based on the special structure analysis of SP1 and SP2 , followed by the complexity analysis and optimum proof. As the implementation of swapping stage m (m = 2,3,…, N) has no impact on Z2 ∗ , it is not given the detailed swapping method for swapping stage m (m = 2,3,…, N). The recovery solutions derived from swapping

291

stage 1 (the following flight sequence of fj 1 is covered by the same aircraft that is reassigned to fj 1 ) are the optimal solutions of SP1 . In terms of the optimality of SP2 , swapping stage m (m ≥ 2), or the aircraft of IN5 and IN6 participating in swapping will apparently increase the number of swapping aircraft in the recovery solutions. Therefore, aircraft rerouting in swapping stage 1 with the swapping intervals IN1 , IN2 , IN3 and IN4 can also get the optimal solution of SP2 . In swapping stage 1, reassignment of aircraft to the flights sequence can be defined as the 0–1 arc flow from one aircraft to one flight which is the first flight of the corresponding flight sequence. Therefore, given one original node o and one termination node t for the flow, a new network can be defined as CN = (o,t,V,A,ω,y,b), where V = I∪J and A = A1 ∪A2 ∪A3 . I and J represent the aircraft set and the flight sequence set, which are originally covered by the aircraft set, of the corresponding swapping intervals. A1 = {(o,i)|i∈I} represent the beginning of reassigning the aircraft to the flight sequences, and A3 = {(j,t)|j∈J} represent the completion of the aircraft to the flight sequences. A2 = {(i,j), i∈I, j∈J}, where each arc (i,j)∈A2 represent a possible reassignment of the aircraft denoted as i∈I to the flight sequence denoted j∈J. There are three labels in each arc (i,j)∈A2 . The first label (denoted as ωij ) refers to the delay time of the first flight in the flight sequence j which is re-covered by the aircraft i, i.e. ωij =max{ai −sj 1 ,0}. The second label (denoted as yij ) represents whether it is the aircraft i, which is reassigned to the flight sequence j, that originally cover the flight sequence j. yij = 0 if the flight sequence j is originally covered by the aircraft i, = 1 otherwise. The third label (denoted as bij ) is the flow capacity of the arc (i,j)∈A2 and bij = 1, ∀(i,j)∈A2 . Similarly, two labels are defined for (o,i) ∈A1 (denoted as yoi and boi )and (j,t)∈A3 (denoted as yjt and bjt ). yoi = 0, boi = 1, ∀i∈I and yjt = 0, bjt = 1, ∀j∈J. In the algorithm, we use binary search method to promote the solution as close to the optimal solution as possible. In each iteration of the binary search algorithm, the arc set A2 is updated by comparison of new threshold values (denoted as ωmin and ωmax) and ωij , (i,j)∈A2 . A new sub-network is then obtained in each iteration. A minimum cost flow problem can be defined with the cost yij of each arc in the sub-network. We then apply the minimum cost path algorithm to the minimum cost maximum network flow problem for the sub-network of each iteration. The pseudo-code of algorithm BSMCF is as follows. Input. IN1 , IN2 , IN3 , IN4 from IN and the corresponding flight sequences originally covered by the aircraft in INi , i = 1, 2, 3, 4. Step 1. Set n =|P’| + |IN1 | + |IN2 | + |IN3 | + |IN4 |, which represents the number of nodes in I or J. Step 2. Construct CN=(o,t,V,A,w,y,b). Step 3. Set ωmax =max{ai −sj 1 | i = j, (i, j)∈A2 }, ωmin = 0, sol= {(i, j)| i = j, (i, j)∈A2 }. Step 4. If ࢜(i, j)∈ A2 satisfying ωmin ≤ωij < ωmax, then let  value1BSMCF = ωmax, value2BSMCF = yi j , then proceed to (i, j )∈sol

Step 12. Step 5. Set a sub-network CN’ = (o,t,V, A1 ∪A2 ’∪A3 , ω,y,b), where A2 ’= A2 /{ωij >(ωmin + ωmax)/2, (i,j)∈A2 }. Step 6. Set val fl = 0 as the initial flow of CN’ from o to t. Step 7. If val fl=n, then proceed to Step 10, otherwise, proceed to Step 11. Step 8. Construct an incremental network CN’(fl) and find path, denoted by M, with minimum cost of yij , from o to t; Step 9. Let C(M) denote the minimum arc capacity of path M, and let θ =min{C(M), n −val fl}; add θ to fl along M in network CN’, and obtain a new flow fl; return to Step 7. Step 10. update sol = {(i,j)| (i,j) ∈ fl ∩ A2 ’} and set ωmax =max( i , j )∈ sol ω ij , to Step 4. Step 11. set ωmin = (ωmax + ωmin) / 2, to Step 4. Step 12. Stop the algorithm.

292

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

Output: sol, value1BSMCF and value2BSMCF . The complexity of the algorithm is O(n4 log(n)). One can consider Step 3, Step 4, Step 5, Step 10 and Step 11 as a binary search algorithm. There are n2 connections between I and J. It requires O(log(n)) to find the optimal value by implementing the binary search algorithm. Step 7, Step 8 and Step 9 describe the minimum cost flow algorithm. The addition of flow is implemented as most n times in a minimum cost path. The minimum cost path can be obtained through Shortest path algorithm, the complexity is not over than O(n3 ). Therefore, the algorithm needs at most O(n4 log(n)). Theorem 2. BSMCF can obtain the optimal solution of SP1 and SP2 , Z2 ∗ =value1BSMCF , Z3 ∗ =value2BSMCF . Proof. First, Step 3, Step 4, Step 5, Step 10 and Step 11 give a binary search algorithm to solve SP1 and SP2 . In each iteration of the binary search algorithm, the problem can be defined as a minimum cost maximum network flow problem with the flow cost yij of each arc (i,j). value1BSMCF that represent the biggest delay time of the flights in the flight sequence (denoted as j) covered by the aircraft (denoted as i) corresponds to the objective of SP1 . value2BSMCF that means the number of aircraft not covering its original flight sequence in the solution of each binary search iteration corresponds to the objective of SP2 . The arc set {(i,j) ∈ sol ∩ A2 }of the solution derived from the minimum cost path algorithm corresponds to the reassignment of the aircraft denoted as i to the flight sequence denoted as j in current iteration of the binary search algorithm. Step 6, Step 7, Step 8 and Step 9 solve this minimum cost maximum network flow problem to get value1BSMCF and value2BSMCF applying the minimum cost path algorithm, which has been proved in Tomizawa (1972). Therefore, the optimal reassignment of aircraft to flight sequence can be obtained by Step 6, Step 7, Step 8, and Step 9. Second, the following proof by contradiction confirms that binary search method can get Z2 ∗ and Z3 ∗ . If we assumed that ∃sol’ ∈SOL, value1sol’ = Z2 ∗ < value1BSMCF ., or value2sol’ = Z3 ∗ < value2BSMCF . From Step 4, ࢻ࢜(i, j) ∈A2 , ωmin ≤ωij < ωmax, ࢺ max( i , j ) ∈ sol‘ ω ij < ω min Set CN” =(o, t, V, A1 ∪A2 ”∪A3 ,w,y,b), where A2 ” = A2 / {(i, j) | ωij > (ωmax + ωmin)/2, (i, j) ∈ A2 }. However, from Step 4 and Step 12, we know that there does not exist a feasible flow fl from o to t in CN”, val fl = n. Therefore a contradiction is shown.  It is mainly concerned about the recovery method from the disruption occurred in a single airport in this section. Faced with the disruptions happening in more than one airport simultaneously, if the recovery process in airports where the disruptions occur are independent with each other, the above recovery method presented in this section can also be applied to get the optimal solution for SP1 and SP2 . However, with the propagation of flight delays, it is more likely that the flights propagated by two disruptions might fly through the same airports before or after swapping aircraft. The recovery processes for different disruptions in more than one airport may have impact on each other. It will lead to a chaotic situation in aircraft swapping. Therefore, no literature can prove that the optimal recovery solution can be derived from recovery processes which may affect to each other under the disruptions in several airports. In terms of the complicated situation for disruptions, especially large scale disruptions, which occur in several airports, metaheuristics will be given for obtaining the near optimal solution set in next section.

4. Solution approach for MOARP In this section, a solution approach (denoted as CNS) for MOARP is described based on a ε -constraints method and heuristics. The ε-constraint method solves a set of constraint single-objective problems obtained by choosing one objective as the only objective to optimize and incorporating inequality constraints for the remaining objectives. The set of problems is obtained by assigning different values to the components of the ε –vector. For each problem in this set, a near optimal solution will be obtained applying some heuristics based on neighborhood search algorithm. First an initial feasible solution is given and then a heuristic based on neighborhood search algorithm is designed to get the near optimal solution. MOARP problem is a three-objective combinatorial optimization problem for which we proposed the mathematical formulation (2)(8), consist in minimizing, respectively, total deviation from the original schedules, the maximal flight delay time and the number of aircraft which do attend in any swapping operations. In the algorithm CNS, we choose objective Z1 as objective function to optimize, as it is the most complicated of the three objectives. Although the objective Z2 that is valued by time may also take noninteger value theoretically, in practical, the delay time can be divided into several equal intervals. If time length of each interval is small enough, such a process is realistic in practical airline operations. Hence, the constrained single-objective problem SP(ε 2 , ε 3 ) considered in the proposed ε -constraint based solution approach is defined as follows.



SP(ε2 , ε3 ) min Z3 (X )

Z2 ( X ) ≤ ε2 s.t. Z3 (X ) ≤ ε3 X ∈D

(34)

Where X = {xri }present the set of decision variables defined in the mathematical model formulation (2-8), and D is the feasible region defined by Constraints (5-8). The proposed method consists in defining a series of ε constraint problems based on a gradual variation of parameters ε 2 and ε 3 . In the algorithm CNS, The optimal solution of each SP(ε 2 , ε3 ) is denoted as opt(Z1 , ε 2 , ε 3 ). Although Abounacer, Rekik, and Renaud (2014) has proved that the optimal solutions of some of these ε -constraint problems are Pareto optimal. Here we give approximate Pareto front of MOARP by generating approximate optimal solution of the identified ε -constraint problems. The following algorithm hereafter show the way of fixing parameters ε 2 and ε 3 to yield the approximate Pareto front. In the algorithm CNS, the parameter ε 3 , which is used in the inequality constraint Z3 (X) ≤ ε 3 , is assigned to different values incremented by one starting the minimum value (Z3 _min) and ending the maximum value (Z3 _max). For each value of ε 3 , the parameter ε2 is then assigned to m different values decremented by ε 2 from the maximum value (Z2 _max) to the minimum value (Z2 _min). The resulting problem SP(ε 2 , ε 3 ) is solved for a given pair (ε 2 , ε 3 ) by heuristics to get a near optimal solution X∗ . Solution X∗ is Pareto optimal if there no exists Pareto optimal solution X’ already identified in the previous iterations such that Z1 (X∗ ) ≥ Z1 (X’ ), Z2 (X∗ ) ≥ Z2 (X’ ) and Z3 (X∗ ) ≥ Z3 (X’ )and. The pseudo-code of algorithm CNS is as follows. Input: original schedules, P, Z2 _min, Z2 _max, Z3 _min, Z3 _max Step 1. Set S =ø, ε 3 = Z3 _min, ε 2 = (Z2 _max − Z2 _min)/m; Step 2. If ε 3 > Z3 _max, then stop the algorithm; otherwise, ε 2 = Z2 _max. Step 3. If ε 2 < Z2 _min, then ε 3 = ε 3 + 1, to Step 2; otherwise, to Step 4.

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

293

Fig. 4. sample for insert the temporary cancelled flight circle into an available aircraft.

Step 4. Set X∗ =opt(Z1 , ε 2 , ε 3 ), if ∃X’∈S, Z1 (X∗ ) ≥ Z1 (X’), Z2 (X∗ ) ≥ Z2 (X’) and Z3 (X∗ ) ≥ Z3 (X’), then to Step 5, otherwise, set S =S∪{X∗ }, to Step 5; Step 5. If ∃X’∈S, Z1 (X∗ ) < Z1 (X’), Z2 (X∗ ) < Z2 (X’) and Z3 (X∗ ) < Z3 (X’), then set S =S/{X’} and ε 2 = ε 2 −ε 2 , to Step 3, otherwise, set ε 2 = ε 2 −ε 2 , to Step 3. Output: S. Based on analysis of the algorithm CNS, we need to address three issues: (1) how to determine certain values of Z2 _min, Z2 _max, Z3 _min, Z3 _max; (2) how to generate an initial feasible solution for the resulting problem SP(ε 2 , ε 3 ) is solved for a given pair (ε 2 , ε 3 ); (3) how to design a local search heuristic to promote improvement of the initial feasible solution. We describe our approaches in the following subsections. 4.1. Evaluation of some parameters Values of parameters Z2 _min, Z2 _max, Z3 _min and Z3 _max refer to search scale of the solution approach and the number of Pareto optimal solutions. Regulations will be given to determine the values of theses parameters by some rules in real time environment. For the range of ε 2 , value of Z2 _min equals to 0, and value of Z2 _max is often evaluated by the flight delay limitation of passenger tolerance. For example, Hu et al., 2015), 2016) give the flight delay time limitation in the mathematical model formulation, and some values of 4 hours, 8 hours, and 12 hours are tested according to the practical situation. In this study, the value of Z2 _max is determined applying similar methods. As for range of the parameter ε 3 , value of Z3 _min is stipulated by the number of aircraft disrupted since the routes of the disrupted aircraft have to be reassigned totally. It is apparent that Z3 _max is not over than the total number of aircraft. 4.2. Initial solution generation The direct results of disruptions are flight delay and/or cancellation if no other available aircraft are attended in swapping. What this section hope to obtain is a basic feasible solution by flight delays and cancellation. A flight can be temporarily cancelled if its delay is greater than ε 2 or if this delay causes a violation of its

origin/destination airport curfew limitation. If such a flight belongs to a flight circle, this flight circle is temporarily cancelled. A flight circle can also be temporarily cancelled if it precedes a flight, not included in a flight circle, whose delay leads to a new departure or arrival time within the curfew period of its origin/destination airport. The resulting initial solution ensures both inequality constraints Z2 (X) ≤ ε 2 and Z3 (X) ≤ ε 3 no matter what is the value of ε 2 and ε 3 . 4.3. Heuristics for solution improvement The heuristic of improving the initial solution to obtain Z1 ∗ is designed based on neighborhood search algorithm. In this section, there are two issues need to described. One is the neighborhood generation and the other is heuristic steps. 4.3.1. Neighborhood generation Before describing the heuristic, some methods of generating neighborhood solution are first presented in details. Method (1). Inserting temporary cancelled flight circles into other available aircraft routes This method can promote the objective value of Z1 decrease apparently. Each circle will tried to be inserted into available aircraft routes (disrupted aircraft routes and other scheduled aircraft routes) holding the constraints Z2 (X) ≤ ε 2 , Z3 (X) ≤ ε 3 and airport curfew time limitation. Generally speaking, there three ways to insert flight circles to available aircraft routes: before the first flight of one route, between flights in one route and after the last flight in one route. Fig. 4 gives a sample for method(1), where the temporary cancelled flight circle {f11 , f12 } can be inserted into total five positions: before f21 , between f22 and f23 , before f31 , between f32 and f33 , and after f34 . However, at most three insert activities (see Fig. 4(a)) can lead to feasible neighborhood solutions: route {f11 , f12 , f21 , f22 , f23 } covered by aircraft 2 (See Fig. 4(b)), route {f11 , f12 , f31 , f32 , f33 , f34 } covered by aircraft 3, and route {f31 , f32 , f11 , f12 , f33 , f34 } covered by aircraft 3 respectively. Because the circle {f11 , f12 } inserting between f22 and f23 or after f34 will violates the constraint Z2 ≤ ε 2 for flight f11 . As aircraft 2 is disrupted, aircraft 2 has been attended in swapping passively even no inserting of {f11 , f12 } before f21 . While aircraft 3 is not disrupted, the insert activities before f31

294

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

Fig. 5. sample for crossing temporary cancelled circle with other flight circles covered by available aircraft.

Fig. 6. Sample of inserting delayed flight circles (flight sequence) into other available aircraft routes.

and between f32 and f33 will lead to the increase of the number of swapping aircraft by one, therefore, the inequality constraints Z3 ≤ε 3 should also be checked simultaneously. Furthermore, the three new routes should be checked to fulfill the airport curfew time as well. Method (2). Crossing one flight sequence (flight circle) with one temporary cancelled flight sequence (flight circle) Sometimes, inserting one flight circle may cause violation of airport curfew time limitation of the following flights. The violation phenomenon can be avoided by cancellation of other flight circles. This insert- cancel activity refers to crossing one flight circle with one temporary cancelled flight circle. Additionally, one flight subsequence in the temporary cancelled circle can also be crossed with one flight sequence with the same origination and terminations. Fig. 5 gives a sample for method (2). There are total four flight circles of aircraft 2 and aircraft 3, {f21 , f22 }, {f31 , f32 }, {f33 , f34 },{f31 , f32 , f33 , f34 },which can be crossed with the temporary

cancelled flight circle {f11 , f12 } respectively. At most four aircraft routes: route {f11 , f12 , f23 } covered by aircraft 2, route {f11 , f12 , f33 , f34 } covered by aircraft 3, route {f31 , f32 , f11 , f12 } covered by aircraft 3, route {f11 , f12 } covered by aircraft 3, can be feasible neighborhood solutions under the constraints of SP(ε 2 , ε 3 ). Fig. 5(b) gives the crossing between the flight subsequence in the temporary cancelled flight circle and one flight sequence. There are total three crossing possibilities between flight sequences: {f11 } and {f31 }, {f12 } and {f32 }, {f12 } and {f32 , f33 , f34 }. Fig. 5(c) gives the result derived from crossing of {f11 , f12 }and {f21 , f22 }. Method (3). Inserting delay flight circles (flight sequence) into other available aircraft routes Similarly with method (1), the flight circle in one disrupted aircraft can also be inserted into other available aircraft route (see Fig. 6(a)) if the disrupted aircraft originally assigned to the flight circle has smaller capacity than the available aircraft. It is noted that flight sequence (no circle) of one aircraft route can be inserted into the tail of the other

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

295

Fig. 7. Sample of crossing two flight subsequences (flight circles).

available aircraft route as long as the departure airport of the flight sequence is the same as arrival airport of the last flight in the available aircraft route (see Fig. 6(b)). Different with general insert activity, the inserting of the flight sequence occur only if two aircraft belongs to the same fleet type. It is reason for promising the aircraft fleet balance at the termination airports. There are at most four feasible neighborhood solutions obtained by method (4) for inserting one flight circle (flight sequence) into the route covered by aircraft 3. One of the solution results is shown in Fig. 6(c). Method (4). Crossing two flight subsequences (flight circles) Crossing between two aircraft routes can promote the decrease of flight delay time apparently. Generally speaking, crossing is limited within two aircraft with the same fleet type. There are at most 7 chances of crossing between aircraft 2 and aircraft 3. The crossing of two flight circles is displayed in Fig. 7(a), two of crossing chances of two flight sequences are shown in Fig. 7(b), and another two crossing chances are crossing {f21 , f22 , f23 }with {f33 }, and crossing {f23 } with {f33 }.one of the crossing result between {f21 , f22 } and {f31 , f32 } is shown in Fig. 6(c). Method 5. Canceling the temporary cancelled flight circles permanently The temporarily cancelled flights circles are permanently cancelled if there is no possibility to insert them into one of the other aircraft by any of the ways previously described. 4.3.2. Heuristic steps Fig. 8. illustrates procedure of the heuristic improving the initial solution based on the above neighborhood generation methods. First temporary cancelled flight circles are sorted by the total cancelled cost of the flight circle. Starting from the top of the list of the temporarily cancelled flights, each flight circle is tried to add to other available aircraft routes until the stopping criterion is met or there is no temporary cancelled flight circle need inserting. If there are still some temporary cancelled flight circles even the stopping criterion for method (1) is met, then the circles (or flight sub-sequences) will be crossed with other flight circles (sequences) in available aircraft routes until the stopping criterion for method (2) is met. If all temporary cancelled flight circles have been added into other available aircraft routes or the stopping criterion for method (2) is met, method (3) and method (4) will be implemented literately until the stopping criterion for method (3) and (4) is met. If the stopping criterion for the heuristic has been met, then

Fig. 8. Heuristic flow for obtaining the optimal solution of SP(ε 2 , ε 3 ).

some temporary cancelled flight circles will be cancelled permanently, and the optimal solution of SP(ε 2 , ε 3 ) is output; otherwise, method (1), (2), (3) and (4) will be operated repeatedly as stated above. 5. Computational experiments To evaluate the methodology designed in this paper, two types of examples are generated. First, a small example with 5 aircraft is shown to present the advantage of Pareto optimal solution of the multi-objective problem. Then, real-world empirical data for a Boeing 737 fleet from a major Chinese airline are used to show the efficiency of the method. The fleet type is composed of 104 aircraft covering approximately 410 domestic flights in daily operation. All computations were performed on a ASUS Y582L workstation with Intel i5-4200 U CPU and 4GB of RAM. All algorithms were coded in java.

296

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299 Table 1 information of the sample problem. plane_ID

available time

flight_ID

departure airport

arrival airport

scheduled departure time

scheduled arrival time

cancel cost

1

12:0 0:0 0

11 12 13 14

PEK XFN PEK SHE

XFN PEK SHE PEK

9:0 0:0 0 14:25:00 16:43:00 18:30:00

10:44:00 16:02:00 17:43:00 19:41:00

104,0 0 0 97,0 0 0 60,0 0 0 71,0 0 0

2

12:10:00

21 22 23 24

PEK HET PEK HET

HET PEK HET PEK

9:58:00 13:40:00 15:15:00 16:45:00

10:46:00 14:33:00 16:04:00 17:34:00

48,0 0 0 53,0 0 0 49,0 0 0 49,0 0 0

3

11:30:00

31 32 33 34

PEK WEH PEK HLH

WEH PEK HLH PEK

12:30:00 14:15:00 16:05:00 18:20:00

13:31:00 15:24:00 17:37:00 19:58:00

61,0 0 0 69,0 0 0 92,0 0 0 98,0 0 0

4

10:20:00

41 42 43 44 45 46

PEK SHA PEK CHG PEK YNT

SHA PEK CHG PEK YNT PEK

11:30:00 13:10:00 14:52:00 16:55:00 19:0 0:0 0 20:55:00

12:27:00 14:11:00 16:13:00 18:20:00 20:06:00 21:51:00

57,0 0 0 61,0 0 0 81,0 0 0 85,0 0 0 66,0 0 0 56,0 0 0

5

9:30:00

51 52 53 54 55 56

PEK TGO PEK YNZ PEK TGO

TGO PEK YNZ PEK TGO HET

11:0 0:0 0 12:52:00 14:40:00 17:42:00 20:50:00 22:42:00

12:11:00 13:59:00 17:02:00 20:08:00 21:52:00 23:42:00

71,0 0 0 67,0 0 0 142,0 0 0 146,0 0 0 62,0 0 0 60,0 0 0

5.1. One sample example A sample problem for 5 aircraft is introduced to evaluate the performance of our method. Table 1 illustrates the basic information of aircraft and flights. The first two columns give the plane ID and the corresponding available time in the beginning of the day. The following columns represent the flight scheduled information (flight_ID, departure airport, arrival airport, scheduled departure time from the departure airport and scheduled arrival time on the arrival airport). Based on the original flight schedules and available aircraft information represented in Table 1, one disruption scenarios were generated as follows. 2 aircraft are disrupted for almost 3 hours in airport PEK: aircraft 1 is not available until 12:00 and aircraft 2 is not available until 12:10. If there are no other recovery measures for the two aircraft routings, then two flights will be delayed: delay time of flight 11 is 180 mins and delay time of flight 21 is 132 mins. According to discussions with Airline, some cost data for recovery operations is shown as follows. Minimum turnaround time is assumed to be 40 minutes, Z2 _max is assumed to be 4 hours, and the recovery period varies from 7:00 AM-12:00 AM. Delay cost for each flight is ¥6/min. The cancellation cost of each flight is closely related to the average ticket price of the flight and its numerical value for this sample example is represented in the last column of Table 1. According to the common solution method (denoted as suboptimal solution method here) in literatures, model presented by (2) and (5)-(8) will be established and then some heuristics are used to solve the model. The methodology has also been simulated in our paper and the corresponding solution is illustrated in Fig. 9. From Fig. 9, we find that if there is no limitation for the maximal flight delay time limitation, extreme delays of flight 21 occurs for more than 2 hours, it may result in passengers’ dissatisfaction of flight 21 and even makes a mess of airline’s daily operations. For this small scale disruption, as all flights are implemented before curfew time, we need to find other available aircraft in airport PEK for swapping stage 1 in order to reduce the maximal flight delay time. In swapping stage 1, swapping intervals are found, |IN1 | = 0, |IN2 | = 0, |IN3 | = 2, and |IN4 | = 1. Table 2 gives the first flight information after the swapping intervals. The flight string after the first flight is implemented by the same aircraft

Table 2 The first flight information of IN2 , IN3 and IN4 . IN

ini

fi1

dep(fi1 )

ai

si1

IN3 IN3 IN4

in1 in2 in3

f31 f41 f51

PEK PEK PEK

11:30 10:20 9:30

12:30 11:30 11:00

Table 3 Comparison of solution results between no recovery method, sub-optimal solution method and BSMCF.

max delay t/min #swapping p total delay t/min #delayed fls #delayed fls > 30min #cancelled fls delay cost/¥ cancel cost/¥ totalcost/¥

initial

suboptimal

BSMCF

180 0 312 2 2 0 1872 0 1872

122 5 152 2 1 0 912 0 912

30 5 357 14 0 0 2142 0 2142

with that of the first flight. All flights strings assigned to any aircraft satisfy the curfew time limitation. We can use BSMCF (designed in Section 3.3) to resolve the flight rescheduling problem under the disruption of the two aircraft in airport PEK. The solution is displayed in Fig. 10. We find that the solution method of minimizing the maximal flight delay time can apparently decrease the maximal delay time. All flight delay time is not over than 30 minutes, which can be neglected and no penalty is for the slight delay time according to the regulation of Civil Aviation Administration of China. Comparison of solution values of no any recovery measure (“initial”), suboptimal solution method (“suboptimal”),and BSMCF (“BSMCF”) is shown in Table 3 to reflect completely the performance of suboptimal solution method and BSMCF. It includes the maximal flight delay time in the solution (“max delay t”), the number of aircraft which attend in the swapping (“#swapping p”), the total delay time of all flights in the solution (‘‘total delay t”), the number of delayed flights in the solution (‘‘#delayed fls”), the number of flights whose delay time is over then 30 minutes

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

297

Fig. 9. Solution derived from sub-optimal solution method.

Fig. 10. Solution derived from BSMCF. Table 4 Near Pareto optimal solution set of the sample problem. Z1 ∗ /¥

Z2 ∗ /min

Z3 ∗

total delay t/min

#delayed fls

#cancelled fls

delay cost/¥

cancel cost/¥

1872 201,732 302,0 0 0 1632 121,822 222,030 1212 101,480 912 101,180

180 122 0 150 132 110 122 80 122 30

2 2 2 3 3 3 4 4 5 5

312 122 0 272 637 505 202 80 152 30

2 1 0 2 2 1 2 1 2 1

0 2 4 0 2 4 0 2 0 2

1872 732 0 1632 3822 3030 1212 480 912 180

0 201,0 0 0 302,0 0 0 0 118,0 0 0 219,0 0 0 0 101,0 0 0 0 101,0 0 0

(‘‘#delayed fls > 30

), the number of cancelled flights in the solution (‘‘#cancelled fls”), flight delay cost (‘‘delay cost”), and flight cancellation cost (‘‘cancel cost”), and total cost (“total cost”) of the solution. From Table 3, we find that the initial solution before using any recovery measure is not satisfied due to its bigger values of maximal flight delay time, total flight delay time and the corresponding flight delay cost. Suboptimal solution method is simply concerned about minimizing the first objective Z1 (the total deviation from the original schedules). In the corresponding solution although the total deviation cost and the total flight delay time from the original schedule has been minimized, one flight has extremely long delay time. BSMCF is only concerned on minimizing the second objective Z2 (the maximal flight delay time). In the solution result derived from BSMCF, although the maximal flight delay time

is very small, the number of delay flights is large. The comparison between the two methods reflects the advantage of multi-objective optimal solution method (described in Section 4) for aircraft recovery problem. The corresponding Pareto solution result set of the sample problem is displayed in Table 4. It includes all possible optimal solutions with iteration of ε 2 and ε 3 . Dispatchers can compare different solutions visually and choose one suitable solution according to real world situations and available resources. 5.2. Real world experiment To further evaluate the utility and efficiency of multi-objective optimal solution method, we also test it on real-world empirical data provided by a major Chinese airline. The dataset is divided into multiple instances, with each corresponding to a randomly

298

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

Fig. 11. Solution result of one instance and relationship between objectives.

and the minimum time of solving all SP(ε 2 , ε 3 ) in each instance. In computational experiments, all test instances are solved around 20 min. CPU time is not over than 4 seconds for all SP(ε 2 , ε 3 ) and even fewer than 0.5 seconds for most SP(ε 2 , ε 3 ). The performance of the method shows that it can be suitable for the most requirements in the practical recovery problem. We also find that the complexity of the problem greatly depends on the number of swapping opportunities. CPU time for SP(ε 2 , ε 3 ) increases with rise of ε 3 , Therefore, when disruptions occur, suitable numbers of aircraft are chosen for swapping according to the scale of the disruption, which will promote an efficient recovery.

Table 5 Sub-fleet substitutions. Sub-fleet

733

73G

73D

738

733 73G 73D 738

– y y y

y – y y

n n – y

n n y –

generated disruption scenario. Sub-fleet substitution feasibility is listed in Table 5, where “y” indicates that the row sub-fleet can substitute the column sub-fleet to cover a flight as scheduled; otherwise, an “n” appears. Fig. 11 gives all solution results of one disruption instance and the relationship between the muti-objective values. There are totally three graphs in the figure. Fig. 11-(1) represents comprehensive relationship of three objectives Z1 , Z2 and Z3 . Although it is hard to show the relationship of the objectives for the three – dimensional graph, clear relationships between Z1 and Z2 , Z1 and Z3 are shown in Fig. 11-(2) and Fig. 11-(3) respectively. In Fig. 11(2), as the maximal flight delay time of the solution increases, the total deviation cost from the original schedules decreases simultaneously. It is understandable for the reason that flights can be delayed longer and it is not necessary to cancel some flights in case of longer delay time. Fig. 11-(3) illustrates the relationship between the numbers of aircraft available for actual swapping and the deviation cost. The deviation cost decreases significantly with an increase in the number of available swapping opportunities. Table 6 illustrates the CPU time of the 10 instances. The first three columns display the number of grounded aircraft (“#grd aircraft”), delayed aircraft (“#del aircraft”) and available aircraft for swapping (“#ava aircraft”). The fourth column represents total CPU time to obtain all near Pareto optimal solutions of each instance and the last three columns give the maximum time, average time

6. Conclusion This paper focuses on multi-objective aircraft recovery problem (MOARP) to normal daily operation after disruption. Three objectives include minimizing the total deviation from original schedules, minimizing the maximal flight delay time, and minimizing the number of aircraft attend in actual swapping. First according to the problem description of MOARP, a multiple objective integer programming formulation is established based on connection network. Then for small scale disruption which occur in the same airport, one polynomial algorithm (BSMCF) is applied to solve the problem, based on the optimal analysis with minimizing the maximal flight delay time and the number of aircraft attend in actual swapping successively. Finally, one heuristic (CNS) is designed, combined large neighborhood search algorithm and ε - constraints method, for large scale recovery problem with three objectives. Through a sample example and some computational experiments for different disruption scenarios, we show that multiobjective solution method is more suitable practical recovery in real world situation compared with traditional suboptimal solution method with only one objective of minimizing the total deviation

Table 6 CPU time of all instances to get near Pareto optimal solutions.

instance instance instance instance instance instance instance instance instance instance

1 2 3 4 5 6 7 8 9 10

#grd p

# del p

#ava p

CPU t/min

max t/sec

aver t/sec

min t/sec

3 3 3 6 6 6 3 3 8 8

27 24 21 21 24 27 30 33 14 17

74 77 80 77 74 71 71 68 82 79

12.2 14.4 14.3 19.9 18.5 19.8 16.9 18.3 20.1 21.4

1.81 2.91 1.99 2.68 2.64 2.89 3.77 3.02 2.66 3.57

0.58 0.60 0.60 0.99 0.93 1.04 0.83 1.01 0.63 0.68

0.19 0.20 0.19 0.19 0.18 0.19 0.19 0.2 0.21 0.21

Y. Hu et al. / Expert Systems With Applications 83 (2017) 283–299

from original schedules, because dispatchers can compare different solutions visually and choose one recovery solution according to real world situations and available resources. For future work, it would be interesting to examine several directions. First, we would like to extend our study to the optimal Polynomial-time algorithm for small scale problem of disruption within more than one airport simultaneously. When the number of disrupted airport is more than one, the swapping stages are more difficult to define and characters of the problem may also change. Second, if we develop an expanded model that integrates the recovery of both aircraft and passengers, improved results would likely to follow. Finally, other approaches, such as meta-heuristics or robust optimization might prove advantageous when multiple recovery options are called for. Although we were able to find optimal solutions, stable solutions are more often preferred by practitioners. Acknowledgements This research was partially supported by the project of Central University Basic Research Fund (HEUCF150903), the project of Heilongjiang province postdoctoral Fund (LBH-Z15047), the project of Chinese Postdoctoral Science Foundation (2016M590276), the Natural Science Foundation of Heilongjiang Province (QC2016095). References Abdelgahny, A., Ekollu, G., Narisimhan, R., & Abdelgahny, K. (2004). A proactive crew recovery decision support tool for commercial airlines during irregular operations. Annals of Operations Research, 127, 309–331. Abounacer, R., Rekik, M., & Renaud, J. (2014). An exact solution approach for multi-objective location–transportation problem for disaster response. Computers & Operations Research, 41, 83–93 2014. Andersson, T., & Varbrand, P. (2004). The flight perturbation problem. Transportation Planning and Technology, 27, 91–117. Ball, M., Barnhart, C., Nemhauser, G., & Odoni, A. (2007). In J. Jan Lenstra, & J. G. Dai (Eds.), Air transportation: Irregular operations and control (pp. 1–68). Atlanta: Handbook in Operations Research and Management Science. Bard, J. F., Yu, G., & Argüello, M. F. (2001). Optimizing aircraft routings in response to groundings and delays. IIE Transactions on Operations Engineering, 33(10), 931–947. Bisaillon, S., Cordeau, J.-F., Laporte, G., & Pasin, F. (2011). A large neighborhood search heuristic for the aircraft and passenger recovery problem. 4OR, A Quarterly Journal of Operations Research, 9(2), 139–157. Bratu, S., & Barnhart, C. (2006). Flight operations recovery: New approaches considering passenger recovery. Journal of Scheduling, 9(3), 279–298. Burke, E., Causmaecker, P., Maere, G., Mulder, J., Paelinck, M., & Berghe, G. (2010). A multi-objective approach for robust airline scheduling. Computers & Operations Research, 37, 822–832. Castro, A. J. M., & Oliveira, E. (2011). Airline operations control: A new concept for operations recovery, 61–97. Chou, T. Y., Liu, T. K., & Lee, C. (2008). Method of inequality-based multiobjective genetic algorithm for domestic daily aircraft routing. IEEE Transaction on System, Man, Cybernetic, Part A, 38, 299–308. Clausen, J., Larsen, A., Larsen, J., & Rezanova, N. J. (2010). Disruption management in the airline industry-concepts, models and methods. Computers & Operations Research, 37(5), 809–821. Dunbar, M., Froyland, G., & Wu, C-L. (2012). Robust airline schedule planning: Minimizing propagated delay in an integrated routing and crewing framework. Transportation Science, 46(2), 204–216.

299

Dunbar, M., Froyland, G., & Wu, C-L. (2014). An integrated scenario-based approach for robust aircraft routing, crew pairing and re-timing. Computers & Operations Research, 45, 68–86. Eggenberg, N., Bierlaire, M., & Salani, M. (2010). Constraint-specific recovery network for solving airline recovery problems. Computers & Operations Research, 37(6), 1014–1026. Etschmaier, M. M., & Mathaisel, D. F. X. (1985). Airline scheduling: An overview. Transportation Science, 19(2), 127–138. Fieldsend, J. E., & Singh, S. (2005). Pareto evolutionary neural networks. IEEE Transactions on Neural Networks, 16(2), 338–354. Filar, J. A., Manyem, P., & White, K. (2001). How airlines and airports recover from schedule perturbations: A survey. Annals of Operations Research, 108(1-4), 315–333. Gershkoff, I. (1987). Aircraft shortage evaluator. ORSA/TIMS joint national meeting. MO: St. Louis. Hu, Y., Song, Y., Zhao, K., & Xu, B. (2016). Integrated recovery of aircraft and passengers after airline operation disruption based on a GRASP algorithm. Transportation Research Part E: Logistics and Transportation Review, 87, 97–112. Hu, Y., Xu, B., Bard, J., Chi, H., & Gao, M. (2015). Optimization of multi-fleet aircraft routing considering passenger transiting under airline disruption. Computers & Industrial Engineering, 80, 132–144. Jafari, N., & Zegordi, S. H. (2010). Simultaneous recovery model for aircraft and passengers. Journal of the Franklin Institute, 348(7), 1638–1655. Kohl, N., Larsen, A., Larsen, J., Ross, A., & Tiourine, S. (2007). Airline disruption management- perspectives, experiences and outlook. Journal of Air Transport Management, 13(3), 149–162. Lan, S., Clarke, J-P., & Barnhart, C. (2006). Planning for robust airline operations: Optimizing aircraft routings and flight departure times to minimize passenger disruptions. Transportation Science, 40(1), 15–28. Lee, L., Lee, C., & Tan, Y. (2007). A multi-objective genetic algorithm for robust flight scheduling using simulation. European Journal of Operational Research, 177, 1948–1968. Liu, T. K., Chen, C. H., & Chou, J. H. (2010). Optimization of short-haul aircraft schedule recovery problems using a hybrid multiobjective genetic algorithm. Expert Systems with Applications, 37, 2307–2315. Maher, S. (2015). A novel passenger recovery approach for the integrated airline recovery problem. Computers & Operations Research, 57, 123–137. Maher, S. (2015). Solving the integrated airline recovery problem using column-and-row generation. Transportation Science, 50(1), 216–239. Rosenberger, J., Johnson, E., & Nemhauser, G. (2003). Rerouting aircraft for airline recovery. Transportation Science, 37, 408–421. Sinclair, K., Cordeau, J., & Laporte, G. (2014). Improvements to a large neighborhood search heuristic for an integrated aircraft and passenger recovery problem. European Journal of Operational Research, 233(1), 234–245. Sinclair, K., Cordeau, J., & Laporte, G. (2016). A column generation post-optimization heuristic for the integrated aircraft and passenger recovery problem. Computers & Operations Research, 65, 42–52. Teodorovic, D., & Guberinic, S. (1984). Optimal dispatching strategy on an airline network after a schedule perturbation. European Journal of Operational Research, 15(2), 178–182. Thengvall, B. G., Bard, J. F., & Yu, G. (20 0 0). Balancing user preferences for aircraft schedule recovery during irregular operations. IIE Transactions, 32(3), 181–193. Thengvall, B. G., Yu, G., & Bard, J. F. (2001). Multiple fleet aircraft schedule recovery following hub closures. Transportation Research Part A: Policy and Practice, 35(4), 289–308. Tomizawa, N. (1972). On some techniques useful for solution of transportation network problems. Networks, 1(2), 173–194. Yan, C., & Kung, J. (2016). Robust aircraft routing. Transportation Science Article in Advance. Yan, S., & Lin, C. G. (1997). Airline scheduling for the temporary closure of airports. Transportation Science, 31(1), 72–82. Yan, S., & Tu, Y. P. (1997). Multifleet routing and multistop flight scheduling for schedule perturbation. European Journal of Operational Research, 103(1), 155–169. Yan, S., & Yang D. , H. (1996). A decision support framework for handling schedule perturbation. Transportation Research Part B: Methodological, 30(6), 405–419. Zhang, D., Yu, C., Desai, J., & Lau, H. (2016). A math-heuristic algorithm for the integrated air service recovery. Transportation Research Part B: Methodological, 84, 211–236.