An Efficient Approximation Algorithm for Aircraft Arrival Sequencing and Scheduling Problem

The aircraft arrival sequencing and scheduling (ASS) problem is a salient problem in airports’ runway scheduling system, which proves to be nondeterministic polynomial (NP) hard. This paper formulates the ASS in the form of a constrained permutation problem and designs a new approximation algorithm to solve it. Then the numerical study is conducted, which validates that this new algorithm has much better performance than ant colony (AC) algorithm and CPLEX, especially when the aircraft types are not too many. In the end, some conclusions are summarized.


Introduction
With the rapid development of airline industry, serious congestions and frequent delays have been hitting most major airports in the world, especially in the United States and Europe [1]. How to enhance the air traffic capacity and reduce the delay becomes a severe problem [2,3].
In 1998, runway has been identified as the primary bottleneck in air traffic [4]; that is, even small enhancements to runway throughput will significantly reduce the delay. However, building more runways is often considered not a realistic option because of practical constraints and huge investment costs. Therefore, many researches and technologists resort to a promising approach, which is to more optimally schedule the aircraft arrival sequence so that the runway can land as many aircraft as possible within a period of time. The optimization process is formulated in this paper as the aircraft arrival sequencing and scheduling (ASS) problem (see in Section 2).
However, ASS is inherently hard to solve [5]; it is nondeterministic polynomial (NP) hard [6,7]. To cope with it, two methods are often adopted, which are mixed integer programming (MIP) and ant colony (AC) algorithm [5,[8][9][10]. ASS can be expressed by MIP formulations. In 1992, Brinton [11] has introduced, as far as we know, the first MIP formulations and designed an implicit enumeration (IE) algorithm to optimize it. In 1993, another MIP is presented by Abela et al. [12] for single-runway ASS. A branch and bound (B&B) algorithm is developed to solve it. Back in 1999, the MIP is presented not only in single but also in multiple runway [8]. Beasley et al. [5] give an improved B&B algorithm by employing linear programming (LP) based tree search. Then Bennell et al. [13] provide an extensive literature overview for ASS.
Ant colony (AC) algorithm [9,10,14] is another effective method for ASS. It is originally proposed by Dorigo in 1992 [15][16][17]. In 2002, Randall [10] presents its first application in ASS, which shows great advantages. Then AC is used to generate initial solutions and to incorporate local search heuristic for single-and multiple-runway ASS [9,18]. Back in 2010, AC is developed to tackle the real-time ASS based on the receding horizon control by Zhan et al. [14]. Experimental results validate that AC is robust, effective, and efficient for ASS.
In this paper, rather than using the above two methods, we develop a new approximation algorithm for ASS. The core idea is to find the lower bound solution of the ASS problem. Then the algorithm presents solutions infinitely approaching this bound. We compare the performance of this new algorithm with AC and MIP (by CPLEX). Computational results verify that this new algorithm returns much better solution and costs less time than AC and MIP, especially when there are several aircraft types. This paper is organized as follows. In Section 2, some constraints for ASS are introduced. The approximation algorithm is proposed in Section 3. In Section 4, ant colony (AC) and MIP are designed for ASS. In Section 5, numerical study is conducted to compare the performance of this new algorithm with AC and MIP (by CPLEX), while some conclusion is summarized in Section 6.

Aircraft Sequencing and Scheduling (ASS) Problem.
ASS aims to make the most use of runway, that is, to minimize the makespan of the landing sequence so that to land as many aircraft as possible within a period of time. The objective function is where end is the landing time of the last aircraft in sequence . For the th aircraft in sequence , its landing time is achieved by where and insure the time-window constraint and MST insures the minimum separation time constraint. These two constraints, as well as some other constraints, are described below. If −1 ≤ −1 is not satisfied, equals +∞.

Minimum Separation Time (MST) Constraints.
MST is a hard constraint to insure safety. When an aircraft flies in the air, it generates wake-vortex (WV). However, WV may result in the instability of the following aircraft (to shake or to lift) [19]. To avoid this, a MST is strictly kept between them.  after the landing of a heavy aircraft. However, when a heavy aircraft lands after a small aircraft, the MST is only 60 s. One reason is that larger aircraft commonly generates and tolerates more turbulent air, while smaller aircraft generates and tolerates less. The asymmetric nature of MST results in the feasibility and necessity of runway scheduling. Proper scheduling strategies can save a lot of landing time. For example in Figure 1, the makespan of sequence 1 equals 288 s; however, for sequence 2, the makespan is only 129 s, which saves more than 50% time.

Time-Window Constraints.
Time-window is a hard constraint to insure that aircraft lands between its earliest and latest possible landing time, which is in the interval time set [ , ]. The earliest possible landing time ( ) depends on the constraints such as maximum airspeed to speed up, runway availability, and possible manoeuvres, while the latest possible landing time ( ) depends on the fuel limitation, maximum allowed delay, minimum airspeed, and so on [13]. Actually, it is not necessary that the time-window for an aircraft is only one that continues interval set [19]. Though we only discuss the single-interval case in this paper, our approach is also applicable to handle the situation that time-window constraints are disjoint intervals.

Precedence Constraints.
Precedence constraints are pairwise requirements to insure whether one aircraft must land before another [19]. There are two reasons for these constraints. One is due to the airline company, which sometimes has strict restriction that one should land first for the reason of priority, banking operations, and so on. Another is due to the jet route, which does not allow two aircraft within the same jet route to overtake each other [20].

Approximation Algorithm for ASS
In this section, we design an approximation algorithm to solve ASS. The core idea is to find the lower bound solution of ASS. Then the approximation algorithm gives ASS solutions infinitely approaching this bound. to a new generated MST (denoted by MST ), which is an approximation of the actual MST. How to determine MST is shown in LP 1 (formulas (3)- (6)). The objective function (formula (3)) is to minimize the differences between MST and MST since ΔMST measures their differences (formula (4)).
Formula (5) calculates out MST . We decompose each element in MST into two parts (see in Figure 2). One is the ability of the leading aircraft to generate the WV (denoted by ), and the other is the ability of the following aircraft to bear the disturbance (denoted by ).
In formula (6), ΔMST , , and are restricted nonnegative. ΔMST ≥ 0 insures that the elements in MST are not bigger than the corresponding elements in MST, which insures Theorem 1.
It is worthy to note that we do not restrict MST ≥ 0, because it does not affect the result in Theorems 1, 2, and 4 from the perspective of mathematical calculation. In addition, it allows more elements in ΔMST to be 0: subject to
Actually, this lower bound is very close to the optimal solution when the time-window constraint is not in consideration since its optimality only depends on the first and last aircraft in the sequence (Theorem 2).  Proof. For the landing sequence ( ) with aircraft, the MST between aircraft and +1 is ( TP( ) − TP( +1 ) ). So we have which is equivalent to is a constant number since the type of each aircraft in is known. So the makespan ( ) is determined by ( TP( 1 ) − TP( ) ), which only concerns the first and last aircraft in .

Definition 3. A sequence is called strongly hypo
if and only if the change of the first and the last aircraft leads to its optimality.
So any sequence in ASS-MST is at least strongly hypooptimal (SHO). Figure 3 lists the relationship between ASS-MST curve and ASS-MST curve when the time-window constraint is not in consideration. ( 1 − 2 ) is very small. We wish to minimize Δ which measures how a solution is close to SHO solution (the lower bound).

Approximation Algorithm Approaches the Lower Bound.
The approximation algorithm gives ASS solutions infinitely approaching the lower bound; that is, it attempts to find a sequence with minimized Δ (see in Figure 3) since Δ measures how a solution is close to a SHO solution (the lower bound) when the time-window constraint is not in consideration (Δ is also a useful measure when considering time-window). The numerical result in Section 5 validates the statement and the efficiency of algorithm.

Contain Time-Window with Aircraft Not Too Many.
In ASS, each aircraft is associated with an estimated landing time (ELT), which is used to estimate the time when the aircraft lands. The ELT relates to the earliest and the latest landing time of the time-window. Generally, the earliest landing time ( ) is one minute less than the ELT, because more than oneminute forward move is often not economical [20]. The latest landing time ( ) can be 60 minutes after the ELT if aircraft have enough fuel and without any emergency accident.
Based on the above discussion, the earliest landing time is the main constraint in time-window, since the latest landing time is often inactive if the number of aircraft is not too many. So we firstly develop an algorithm for ASS without too many aircraft (often less than 30). And then we extend it to solve the case with many aircraft (often more than 30).
The core idea is to infinitely approach the lower bound, that is, minimizing additional makespan when comparing ASS-MST and ASS-MST . The additional makespan consists of ΔMST and the earliest landing time ( ). We give the approximation of the additional makespan.
For the FSFC sequence ϝ, all aircraft with the same type are put together according to their order in ϝ to construct subsequences. There are subsequences, which are ϝ 1 , ϝ 2 , . . . , ϝ , where TP(ϝ ) = , ∈ Λ, and ϝ is the th aircraft in ϝ .

Proof.
Step 2 costs ( ) to find the smallest Δ . Since there are aircraft, Step 2 repeats for times. So the complexity is ( ).

Contain Time-Window with Many
Aircraft. If the number of aircraft is too large (often more than 30), the "inactive" assumption of the latest landing time is often broken, since some aircraft may be assigned to land 1 hour later. To avoid this, we firstly reschedule the first 30 aircraft (from the 1st to the 30th) in FCFS by the above algorithm and then reschedule the next 30 aircraft (from the 31st to the 60th) and then the next 30 aircraft (from the 61st to the 90th) and go on until all aircraft have been considered. In the end, we connect all the rescheduled sequences one by one to construct the ultimate sequence. Here the length of 30 aircraft is based on the numerical result in Section 5. Of course we can set other lengths for subrescheduling but without such a good solution. These whole processes do not exceed the complexity of Theorem 4.

Contain Precedence.
To consider the precedence, we may freeze the subsequence. For example, aircraft in subsequence ϝ should land before aircraft in ϝ ; we can freeze subsequence ϝ when is in the first position of ϝ . The only sufficient condition to unfreeze ϝ is that aircraft has been transferred into the final sequence Π.

Ant Colony (AC) for ASS. Ant colony (AC) algorithm
is an intensively studied method [13]. Many papers adopt it for single-runway scheduling [9,10,14]. Here AC is used to find a solution for ASS, which is a valid comparison with the approximation algorithm. The following are some important settings of AC in the computation. To make sure most routes returned by ants satisfy the timewindow constraints, we restrict the allowable aircraft to the integer set {max{1, − }, . . . , min{ , + }} when an ant does its th visit. This principle is also called -CPS [19]. Under this principle, an aircraft should be visited before an ant does its ( + + 1)th visit. So if aircraft has not been visited when an ant does its ( + )th visit, it is forced to be visited by this ant. In the numerical study here, setting the integer number = 5 is good to satisfy time-window.

State Transition Rule. After finishing visiting an aircraft
, all ants choose the next allowable aircraft according to the probability of Pb( , ): 4.1.3. Pheromone Updating Rule. The pheromone updating operation is carried out on each arc in the completed tour as in formulas (12) and (13) for each ant: where (1 − ) ( ) models the evaporation of the pheromone and Δ is the lately released pheromone by ants. For an ant , Δ is defined as where is a given constant and is makespan of tour route by ant .

Mixed Integer Programming (MIP) for ASS. In 2000,
Beasley et al. [5] proposed a MIP for the single-runway scheduling problem. Here we extend it to ASS as a solution comparison of the approximation algorithm. The following are some important settings.
The objective function (formula (14)) is to minimize the makespan.
≥ (formula (18)) insures is the landing time of the last aircraft.
Formula (16) is the separation constraint. There are two cases.

Numerical Study of Three Methods
Here is the numerical result of the approximation algorithm for ASS. The AC and CPLEX are used as a valid comparison of this new algorithm.

Randomly Generate the MST.
In ASS, the concerning MST is important. Table 1 gives a typical MST table by a classical aircraft types classification. However, there are some other MST tables by other classifying principles. For example, Table 2 illustrates another MST table with more than 3 aircraft types. We believe that our algorithm applies almost all possible MST matrixes. So in the numerical study we randomly generate MST matrixes which have two common properties.
For MST matrixes with aircraft types, their general properties are summarized as follows.
In numerical study, MST can be generated by where MST 1 is a given number and (≥ 1) is a random number. In the following numerical study, if ≤ 4, is in the interval of [1,2] satisfying uniform distribution; if > 4, is in [1, 1 + 1/( − 3)] satisfying uniform distribution. This setting insures MST 1 is not too much bigger than MST 1 . It is in accord with the actual MST; for example, Tables 1 and 2 are basically satisfied.
It is worthy to note that even if a MST does not strictly satisfy the above settings, it can also be well coped with by the approximation algorithm.

Numerical Result.
We test the cases of = 3 and = 9. In each case, ten groups (aircraft number = 20, 40, 60, . . . , 200) of aircraft sequences are randomly generated. Each group contains 20 random sequences; the number of all types of aircraft is the same in each sequence. The MST is randomly generated for each sequence (see detailed setting and its parameters explanation in formula (20)).
The estimated landing time (ELT) of aircraft is simulated by a Poisson arrival process, which has been validated by Willemain et al. [21]. The expectation of interval time between neighbor aircraft in ELT is set to be (4/5) (MST), where (MST) is the mean value of all elements in MST. It models how crowded an airport is. The earliest landing time is set to be one minute less than the ELT, because more than a minute forward move is often not economical [20]. The latest possible landing time is set 60 minutes after the earliest landing time.
For approximation algorithm, long aircraft sequences (with more than 30 aircraft) are scheduled each time for 30 aircraft (see the explanation in Section 3.2.2). For AC algorithm, the termination criteria are the maximal generations of 100. The parameters in formulas (11)-(13) are set as = 1, = 5, = 0.1, and = 100. The allowable aircraft are restricted in set {max{1, − 5}, . . . , min{ , + 5}} when an ant does its th visit. For MIP, the traditional B&B search in CPLEX is used and the maximal tree nodes of B&B search are set to be 4000. The numerical result is shown in Table 3 with average runtime and average runway enhancement of the tested 20 sequences in each group. For each sequence, the "runway enhancement" is achieved by It is comparison of the makespan between the final sequence Π and the FSFC sequence since FCFS is a widely used method in current runway scheduling system.
When there are fewer aircraft types (3 types), the approximation algorithm enhances the runway throughput more than AC and CPLEX and costs less time. It validates the performance of the algorithm. Actually classifying the aircraft into 3 types is often enough. If there are too many aircraft types, we can classify the aircraft into several big categories and reconstruct the MST table concerning the big categories other than the aircraft types-each big category contains aircraft types with very close properties. So the ASS in this condition can also be well coped with by the approximation algorithm.
When there are 9 aircraft types, the approximation algorithm generally has higher runway enhancement when the number of aircraft is not too many (<100 aircraft). For too many aircraft (100-200 aircraft), AC enhances the runway throughput more. However, the approximation algorithm costs less time than the other two methods.
By comparing AC and MIP, we find that MIP (by CPLEX) has better performance when both the number of aircraft types and aircraft are small. For a large number of aircraft (or many aircraft types), it is better to use AC for scheduling other than CPLEX. It is also worthy to note that the approximation algorithm has better performance than CPLEX in almost all the enumerated cases in Table 3.

Conclusions
To cope with aircraft sequencing and scheduling (ASS) problem, we design an approximation algorithm, which is infinitely approaching the lower bound of the optimal solution. Numerical result validates that this algorithm, with less runtime, can get a smaller makespan than AC and CPLEX, especially when the number of aircraft types is not too many. This is a big increase of the runway.

Nomenclature and Notation
Nomenclature AC: Ant colony ASS: Aircraft sequencing and scheduling ASS-MST: ASS problem concerning MST; that is, in ASS the separation between aircraft is determined by MST matrix ASS-MST : In ASS the separation is determined by MST matrix FCFS: First come first served LP: Linear programming MIP: Mixed integer programming MST: Minimum separation time (matrix). MST denotes the MST of an aircraft of type followed by another aircraft of type . MST ⩽ MST + MST , ∀ , , ∈ Λ MST : A newly generated MST matrix. MST denotes the MST of an aircraft of type followed by another of type ΔMST: Differences between MST and MST . ΔMST = MST − MST . ΔMST denotes the ΔMST of an aircraft of type followed by another of type . All elements in ΔMST are nonnegative SHO: Strongly hypooptimal WV: Wake-vortex.
Notation allowed : A set of allowed aircraft for ants to visit followed by aircraft , : Parameters determining the relative importance of versus Themakespanofsequence in ASS-MST, that is, the landing time of the last aircraft in sequence in ASS-MST : Themakespanofsequence in ASS-MST Δ : Additional makespan from to ; that is, Δ = − Δ : Additional makespan of adding aircraft of type to a sequence : The landing time of aircraft : The ability of the leading aircraft of type to generate WV : The ability of the following aircraft of type to bear disturbance Λ: Theaircrafttypeset : Alargepositiveconstant ϝ: The FSFC sequence. ϝ is a subsequence where all aircraft of the same type are put together according to their order in ϝ : An aircraft landing sequence. is the th aircraft in sequence , and end is the last aircraft in Π: The final aircraft landing sequence. Π is the th aircraft in Π : Aircraft lands before aircraft ( = 1) and after ( = 0) ( , ): Pheromone in arc( , ) in AC. 0 ( , ) is initial pheromone ( , ): Heuristic information in arc( , ), where ( , ) = 1/MST TP( )TP( ) : Parameter modeling the pheromone evaporation ratio.