Airports throughout the world have seen a rapid increase in the number of flights and passengers over the past decade. Not withstanding the current decline due to COVID-19, the international Air Transport Association ( [17] ) expects 7.2 billion passengers to travel in 2035, a near doubling of the 3.8 billion air travelers in 2016.
Airport facilities, in particular, terminals were not sized to handle such traffic. Although in the long term, expansions in capacity have been planned, such expansions in the near future is unlikely. Consequently, airports have seen an increase in delays as well. The Federal Aviation Administration reports that the cost of delayed flights in 2018 increased by 6 percent, from $26.6 to $28.2 billion, an increase of $1.6 billion ( [16] ). As a result, better management of airport operations is needed.
One critical area of airport operations is the optimal and efficient assignment of arriving flights to gates.
The airport flight-to-gate assignment problem has been extensively studied by many researchers in both the operations research and aviation communities. Many models and algorithms have been proposed. For a detailed review of past work in this area, we refer the readers to [3] and [9] , but give a brief overview of the popular objectives and the common solution methods. There are three popular objectives to consider in the models.
The first one is the maximization of passenger satisfaction levels. In particular, researchers consider walking distances and waiting or transit time as proxies for the passenger satisfaction level. A more detailed approach involves a differentiation between the transfer passengers and destination passengers. In these models, the costs associated with the assignments of two flights to any two gates reflect the distance between the two gates and the connection time between the two flights. Interested readers are referred to [24] , [6] , [30] , and [21] . The second class of objective focuses on airport and airline operations. Some papers, such as [20] and [7] , consider the minimization of the number of ungated flights or equivalently the number of flights assigned to the remote gates or introduce costs only incurred when the flights are assigned to remote gates to differentiate the priorities of the flights. Moreover, [19] and [25] point out that airlines likely have preferences over which set of gates to park their flights leading to the consideration of maximizing total gate preference scores. The third popular objective is to improve the robustness of the solution to schedule variations. Due to the uncertain nature of the air transportation system, arrival times and departure times are likely to be stochastic subject to various factors such as weather and maintenance. Early arrival and late departure may cause the arriving flight to wait for the departing flight at the assigned gate to push back. This phenomenon is referred to as a gate conflict. Other studies consider robustness with respect to minimizing the expected gate conflict time. In the event of major disruptions such as severe weather or maintenance delays, airport has to be able to quickly recover from the deviations from the original arrival schedules and adjust the flight-to-gate assignments to reduce any potential delays. Readers are referred to the work of [13] and [23] . Since the above mentioned three objectives are all very important to take into account, some combinations of them are also explored. The weighted sum method is more common in such line of considerations, but a pareto-based approach has been utilized as well. [22] and [14] are examples of studies that use weighted sum approach while [8] approaches his problem with a pareto local search method.
Due to the general N P-hard nature of the problem, it can be extremely difficult to solve the problems directly and exactly. Many solution algorithms and heuristics have been proposed. In much of the work to date, an integer programming or mixed integer programming formulation is first presented and solver is called to obtain a solution and validate the model. These formulations are large and complex because of the many complicated constraints required to reflect real world phenomena. Consequently, it is challenging to improve the solution methodologies. Readers are referred to [23] and [27] for the details of these models. In addition to integer programming or mixed integer programming approach, [26] use a stochastic optimization model to capture the uncertain nature of the problem and study the occurrence of gate conflicts. [29] use a column generation method together with a branch-and-bound algorithm to obtain integral solutions. They give a description of the column generation method and a general strategy to apply the branch-and-bound algorithm. Besides exact algorithms, heuristics are extensively used. Popular heuristics such as genetic algorithms and tabu search algorithms have been explored. Tabu search algorithm performs a local search around the current solution and prohibits searching at previous search points by keeping a list of those points for a few iterations while genetic algorithm is inspired from the evolutionary idea where a chromosome of higher fitness is more likely to survive and it involves operations of selections, crossover, and mutations. For an implementation of these two metaheuristics, readers are referred to [21] . Moreover, [31] implemented a MIP-based heuristic which is adapted from many popular heuristic and they show that the adaptation gives a good performance. Lastly, there are some other heuristics based on neighborhood search, breakout local search, Fuzzy Bee Colony optimization (FBCO), and particle swarm optimization (PSO). Readers are referred to [32] , [2] , [2] , and [10] for these heuristics.
In this paper, we focus on the airport operations and aim to minimize the arrival delays over all the flights.
We propose an exact solution method that is capable of obtaining the flight-to-gate assignments at busy hub airports where large numbers of flights are seen on a daily basis. Our main contributions are summarized below.
1. Column generation formulation. We present an explicit column generation formulation with a set covering master problem and we decompose the pricing problem into independent problems to which we apply both heuristic and exact algorithms to obtain favorable assignment patterns. Such decomposition allows us to simplify the pricing problem formulations and reduce the size of the individual pricing problem significantly. Experiments show that this formulation is far more efficient than a conventional compact mixed integer programming formulation.
2. Approximation and exact algorithms. We explore a few approximation methods and exact algorithms to efficiently solve individual pricing problems after the decomposition. We apply transformations on the pricing problem formulations to understand the underlying structures, so we can devise new algorithms and adapt existing heuristics. In particular,
• We show and exploit the submodularity of the objective functions of the pricing problems to design an approximation algorithm. To the best of our knowledge, this is the first attempt to utilize submodularity to efficiently solve the pricing problems in a column generation setting.
• We derive an exact dynamic programming algorithm based on a recursive formula and show that the direct implementation of the algorithm works reasonably well in practice. We show the algorithm to be pseudo-polynomial in the special case of integer inputs and present a new heuristic based on this observation.
• We utilize a rolling horizon framework and propose another 2-approximation block decomposition heuristic for large-sized instances.
The structure of the rest of the paper is as follows. In Section 2, we provide a description of the flight-to-gate assignment problem. In Section 3, we provide the formulation for the column generation algorithm, and the exact algorithms and heuristics used for solving the pricing problems are discussed in Section 4. Specifically, a submodular maximization approximation algorithm is described in Subsection 4.2 and a dynamic programming algorithm with complexity analyses is presented in Subsection 4.3. In Subsection 4.4, we discuss the implementation of the rolling horizon framework and present a new block decomposition approximation for large-sized instances. We discuss how to retrieve integer feasible solutions in Section 5. Results of the computational experiments are presented in Section 6. Finally, we make some concluding remarks and discuss possible extensions of this work in Section 7 .
The flight-to-gate assignment problem may be described as follows. Let the set of flights be F and the set of gates be G. Each flight i ∈ F has an arrival time, a i , an airline, and an aircraft type which determines its minimum turn time, τ i , the minimum duration for which the flight must remain after parking at its assigned gate before being ready for departure. Each gate k ∈ G has a buffer time, b k , that has to be observed between consecutive flights, and a set of allowable aircraft types. The solution of the problem gives an assignment in which each flight is assigned to one gate. For each flight i, we define its arrival delay as the difference between its arrival time, a i and the park time at its assigned gate, t g i . The decisions are the gate assignment, x ik , the park time, t g i , and the push back time, t p i , for i ∈ F and k ∈ G. In particular, the decision variables x ik are equal to 1 if flight i is assigned to gate k and otherwise equal to 0. In addition, we also consider the notion of compatibility, α ik , where a flight can only be assigned to a compatible gate. In particular, α ik equal to 1 if flight i is compatible to gate k and equal to 0 otherwise. The compatibilities are determined by various factors including gate types (heavy or regular), flight types (heavy or regular), and airline preferences. Two additional conditions have to be met for an assignment to be valid: 1. A flight occupies its assigned gate for a duration of at least its minimum turn time, τ i , for i ∈ F.
A pair of flights that are assigned to the same gate cannot occupy the gate at the same time and the buffer time, b k , must be observed before the gate becomes ready for the next flight.
We further assume independence of the gates which exclude interferences between the gates. Based on the above descriptions, we give a compact mixed integer programming formulation below,
where M is a sufficiently large constant. The expression in (2.1) gives the objective function which minimizes the total arrival delays. Constraint (2.2) ensures that each flight is assigned to exactly one gate. Constraint Constraint (2.4) ensures that the park time, t g i , is no earlier than the arrival time, a i . We state constraint (2.4) for a mathematical convenience in the model as flights can arrive ahead of their schedules in real world. Note that we do not impose an upper bound on the park times, t g i , since it is possible that an arriving flight has to wait at its assigned gate for a long period of time for the departing flight that is currently occupying the gate to pushback due to various reasons, such as weather, technical difficulties, although the chance of such long wait is slim. Condition 2 above leads to constraint (2.5) which observes the buffer time, b k . It considers all pairs of flights i < j. When they are assigned to the same gate, the difference between t p i and t g j must be at least the buffer time of the gate, b k . Constraint (2.6) is the compatibility constraint to ensure the flights are only assigned to the compatible gates. Constraints (2.7) and (2.8) are binary and non-negativity requirements on the decision variables respectively. We will show in computational experiments that this compact formulation is not ideal to obtain the desired flight-to-gate assignment efficiently. We instead propose a column generation approach to solve the problem. Note also that this setup is in fact similar to what airlines do for their flight-to-gate assignments in the real world. Interested readers are referred to [28] .
Constraint (2.2) in the compact mixed integer programming formulation suggests a set partitioning master problem should be used in the column generation formulation. However, as many papers, such as [1] and [11] , have pointed out, a set covering master problem formulation is numerically more stable when solving its linear programming relaxation compared to a set partitioning master problem. Therefore, we present a set covering master problem formulation where the constraint in which each flight is assigned to exactly one gate is relaxed to a constraint in which each flight is assigned to at least one gate. Since this gives a relaxation and we are minimizing the total arrival delay that is non-increasing after removing a flight from a gate, the set covering master problem either produces a set partitioning assignment or a set covering assignment with a total arrival delay of zero. In the latter case, we can fix each flight with multiple assignments to one of its assigned gates and obtain a new assignment that has a total arrival delay at most as large as that of the original assignment.
We define the following parameters: Note that δ k ip and c k p are constants that can be computed for a given assignment pattern. The decision variables z k p (k ∈ G, p ∈ P k ) are equal to 1 if pattern p of gate k is used and 0 otherwise. Then the set covering master problem is given by
The objective (3.1) is an expression for total arrival delays by summing the arrival delays of all patterns over all gates. Constraint (3.2) is the cover constraint that ensures each flight is assigned to at least one gate. Constraint the binary requirements on the decision variables z k p (k ∈ G, p ∈ P k ) and relaxing these constraints to 0 ≤ z k p ≤ 1 gives the linear programming relaxation of the set covering master problem. The linear programming relaxation is then solved in each iteration.
In the pricing problem, we construct feasible assignment patterns for each gate. If we associate dual variables π i with (3.2) and µ k with (3.3), the reduced costc k p (k ∈ G, p ∈ P k ) of a pattern p of gate k is given bȳ
Given a feasible solution z to the linear programming relaxation of (3.2) -(3.4), we know that z is optimal if and only if the optimal assignment patterns generated from the pricing problem have non-negative reduced costs. Then the pricing problem is given by
Since we assume each gate is independent, we can further decompose the pricing problem into |G| independent pricing problems, one for each gate. The k th pricing problem is given by When the pricing problems are solved to generate a new assignment pattern p, the values of the constants, δ k ip , are determined using the values of the decision variables x ik and c k p is the sum of all individual arrival delay, t g i − a i . Let M denote a large constant, then the k th pricing problem is given by
14)
The objective function (3.8) is obtained by expressing c k p explicitly in decision variables t g i and parameters a i , and substituting the parameters δ k ip by the decision variables x ik . Constraints (3.9) -(3.14) follow directly from the constraints (2.3) -(2.8) in the compact mixed integer programming formulation. However, we do not need constraint (2.2) . For flights that are not assigned to this gate, we assume their arrival delays to be zero.
The solutions x ik from a pricing problem form a potential assignment pattern. If the corresponding reduced cost given by the objective value is negative, then this assignment pattern is favorable and added to the master problem.
We want to point out that µ k is a constant in the k th pricing problem and thus dropping this term we can equivalently solve (3.8) -(3.14) by solving the following maximization problem:
subject to (3.9) − (3.14).
This equivalent problem is intuitively easier to interpret as we may consider π i as the benefits of accepting flight i at the cost of an arrival delay t g i − a i . Consequently, we are attempting to maximize the total net benefits over all flights. We refer to this maximization problem as the pricing problem for any gate in all following sections and use the terms total net benefits and reduced cost interchangeably.
4 Solving the pricing problem 4 
The decomposition of the pricing problem into smaller pricing problems allows additional pre-processing of the set of flights F for each gate. Consider a particular gate k, we can construct a subset of flights F ⊆ F to include only compatible flights and flights with positive dual variables π i . To see this fixing is correct, the incompatible flights are not allowed to be assigned to this gate while accepting the flights with π i = 0 does not contribute any benefit at the cost of arrival delays of other flights. For the rest of the discussions in this section, we use the subset of flights F and assume |F | = n. Note that for different gates, the set F could be different.
Submodularity is a well studied property of set functions and there are many existing algorithms that can approximate the task of maximizing a submodular function very efficiently and effectively. We notice that, in the context of the flight-to-gate assignment problem, the submodularity of the objective function (3.15) is very intuitive. Consider any gate and given a set of accepted flights at this gate, the impact of accepting an additional flight on the total arrival delay of this set of flights is at least as large as that of accepting the additional flight with a subset of those flights.
Formally consider the pricing problem for gate k and recall that our objective after pre-processing the set of
For simplicity, we denote t g i − a i by t g i and clearly t g i ≥ 0. We define a function f : 2 F → R as follows,
We further denote the variables t g i under a set of assigned flights A by an additional superscript A and assume that the values of t g,A i for i / ∈ A to be zero. By the definition of f , we have that f (∅) = 0. For a given set of assigned flights A, the constraints (3.9) -(3.11) can be incorporated when evaluating f (A). Specifically, let the flights in the set A be 1, 2, . . . |A| and the total benefits i∈A π i can be computed. What remains is to compute the arrival delays. Note that the first flight parks at the gate at its arrival time, a 1 , and pushes back after a duration of its minimum turn time at a 1 + τ 1 . The gate becomes available for the second flight after a duration of the buffer time at a 1 + τ 1 + b k . Consequently, an arrival delay of a 1 + τ 1 + b k − a 2 on the second flight is incurred if the second flight arrives earlier than a 1 + τ 1 + b k and it parks at a 1 + τ 1 + b k . Otherwise, there is no arrival delay on the second flight and it parks at its arrival time, a 2 . The same procedure can be applied to the rest of flights in the set A to compute f (A).
Note that since the subset of flights F are all compatible to this gate and the above procedure incorporates constraints (3.9) -(3.11), we can pose the pricing problems as the following equivalent unconstrained maximization problem
In order to show the submodularity of the function f , we need the following observation about the variables
proof of Lemma 4.1. We can rewrite the left hand side expressions into
The last inequality is valid because A ⊆ B, for any i ∈ A, thus the arrival delay of i in the set of assigned flights B of larger cardinality is at least as large as its arrival delay in the set of assigned flights A ⊆ B. Additionally, since t g,B i ≥ 0, we have the desired inequality.
We next show that f is a submodular function. 
Similarly, we have that
Since A ⊆ B, the arrival delay of u in the set of assigned flights B is at least as large as that of u in the set of assigned flights A and then t
. For the first term in the expressions (4.2) and (4.3), we consider the following cases when accepting an additional flight u to A and B:
1. The park times remain unchanged for all i ∈ B and thus i ∈ A. Both terms are zero. 
The two cases above verify that when we accept an additional flight u,
and this shows that f is a submodular function.
We want to point out that f is not monotonic in general. Suppose that A ⊆ B ⊆ F , then
is always non-positive by Lemma 4.1. However, since π i ≥ 0 for
is not necessarily always non-negative or non-positive and thus f is not monotone in general.
Now we can utilize an existing submodular maximization algorithm by [4] which is shown below as algorithm 1.
In this algorithm, we select each flight in the set of flights F with a probability. Let X i and Y i be two random sets to be updated in each iteration and we initialize X 0 = ∅ and Y 0 = F . We also denote the elements in F by 1, 2, . . . , n.
end 11: end for 12: return X n (or equivalently Y n )
The following theorem establishes a theoretical approximation guarantee of Algorithm 1. (4.1), and OP T be the true optimal solution of the problem (subMax) and X n (or equivalently Y n ) is the set returned by the algorithm.
proof of Theorem 4.1. The proof is a very minor modification of what is given in [4] and attached in Appendix A for a reference and completeness. The modification is since we require
Although it is possible that the condition on f to obtain the theoretical guarantee is not satisfied, the assignment patterns that are generated can still have favorable total net benefits. We discuss later in the numerical experiment section about the implementation of this algorithm and its performance.
Approximation algorithm can be efficient in generating favorable assignment patterns, however, exact optimal solutions to the pricing problems have to be obtained to prove optimality at each node. This would typically be done using a general integer programming solver. However, this can be slow. To reduce overall solution time,
we have developed and now present a much more effective dynamic programming algorithm to achieve this task.
For any gate k, we define an auxiliary function g i (t) to be the maximum total net benefits from optimally accepting flights in the set {i, i + 1, · · · , n} where the smallest indexed accepted flight from {i, i + 1, . . . , n} can be accepted any time after time t. To give a formal definition,
As an illustration, we present an example of a set of 8 arrivals and show in Figure 1 the functions g i (t). In particular, the left plot shows the function g 1 (t) while the right plot shows the function g 8 (t).
We observe that the optimal value of the k th pricing problem is g 1 (0). We further assume the buffer time of 
Based on the definition of g i (t), we can derive the following recursive formula, We consider each case separately.
1. This is the terminating case. If i ≥ n + 1, we do not have any flights and thus receive a zero net benefit.
2. Since a flight cannot park earlier than its arrival time, the total net benefits for any t ≤ a i are equal to
3. If the current time t > a i + π i , accepting the flight i does not contribute to the total net benefits as t is equivalent to the procedure EVALORACLE(i, t) with input n and t which further calls the procedure EVALORACLE(i, t) at most 2 times with input n + 1 and t or n + 1 and t + p n . This gives the base case.
Suppose that the statement is true for i = n, n − 1, . . . , j + 1. Now consider i = j and g j (t) for any t, the procedure EVALORACLE(i, t) with input j and t further calls the procedure EVALORACLE(i, t) at most 2 times with input j + 1 and t or j + 1 and t + p i . The former is equivalent to the value of g j+1 (t) and the latter is equivalent to the value of g j+1 (t + p i ). The inductive step shows that both call requires at most 2 n−j further recursive calls to the procedure EVALORACLE(i, t). Therefore, g i (t) for any t requires a total of 2 · 2 n−j = 2 n−j+1 recursive calls, which completes the proof. if t ≤ a i then 6: return EvalOracle(i, a i ) 7:
else if t ≥ a i + π i then 8: return EvalOracle(i + 1, t), S ← S ∪ {0} 9: Recall a direct implementation of the dynamic programming algorithm terminates with a worst case 2 n recursive calls as proven in Proposition 4.1. We can improve the complexity significantly if we further assume all input data are rationals and run a backward version of the dynamic programming algorithm. Input data includes all input parameters to the pricing problems, namely the arrival times, a i , the dual information, π i , and the processing time, p i . Without loss of generality, we assume the input data are expressed as rationals with a common denominator. We introduce a new notion of last time to consider for the pricing problems, denoted by c. This is the time beyond which no further flights are accepted. It suffices to choose c as a n + max 1≤i≤n π i .
For any flight j, if we consider the function g j (t) with t > c, we get that t > a l + π l for j ≤ l ≤ n. Consequently the function value g j (t) = g l (t) = g n+1 (t) = 0 and x lk = 0 for j ≤ l ≤ n.
For every i ∈ [n], we see that the function g i (t) is piecewise linear and we label the points where the slope of the function changes as the breakpoints of the function. When we have integral input data, each g i (t) does not have too many breakpoints. The following theorem gives a formal statement about this observation. To prove Theorem 4.2, we first consider the possible slopes for a particular function g n−i (t). π n if t ≤ a n a n + π n − t if a n < t ≤ a n + π n 0 if t > a n + π n .
(4.6)
The set of slopes of function g n (t) is {0, −1}. This is the base case of n. Suppose the statement is true for n, n − 1, · · · , n − i. Then the set of possible slopes of g n−i (t) is {0, −1, . . . , −i − 1}. Now consider the function g n−i−1 (t), based on the recursive formula (4.5),
(4.7)
In the second case, the set of possible slopes of g n−i−1 (t) is the same as that of g n−i (t). In the third case, proof of Theorem 4.2. Consider a function g n−i (t), from the recursive formula (4.5), we observe that all the potential breakpoints can be computed by equating the two components in the max expression. The solution of t may constitute a breakpoint of g n−i . As a result, we consider the following linear equation (a n−i + π n−i − t) + g n−i+1 (t + p i ) = g n−i+1 (t).
(4.8)
Functions g n−i+1 (t + p i ) and g n−i+1 (t) are piecewise linear in variable t. We can write
where c 1 and c 2 are rational constants and k 1 and k 2 are integer constants. We only consider cases where k 1 + 1 = k 2 , otherwise the two functions on the left hand side and right hand side of (4.8) have the same slope and do not constitute a breakpoint. Now equation (4.8) can be rewritten as The total number of breakpoints is then at most 1dc + 2dc + 3dc + · · · + (i + 1)dc = O(n 2 dc), which completes the proof.
Note that since each of the functions g i (t) has at most O(n 2 dc) breakpoints, we only need to evaluate the function at these potential breakpoints to construct the function and consequently we need to run in O(n 3 dc)
to construct all functions g i (t) in the interval [0, c]. In particular, note that only the function values at the points in the list of the potential breakpoints are needed since the functions are all piece-wise linear. We can first construct the function g n (t) in the interval [0, c] by evaluating the function at the points in the list of the potential breakpoints given in (4.11) in reverse order and then we construct the function g j (t) for 1 ≤ j ≤ n − 1 using previously constructed functions g l (t) for j + 1 ≤ l ≤ n based on the recursive formula (4.5).
Although the above analyses provide a simpler complexity in the case of rational input data, it is not very effective in practice as the common denominator d can be large after scaling the input data.
We now propose an implementation of the dynamic programming algorithm with a running time of O(nc) in the case of integral input data. Recall that when the input data are rationals, we can construct the each function Next, we optimally assign the flights 1, 2, . . . , n in the set F based on the recursive formula (4.5). For any flight j, if we accept j, it can only be accepted at an integral time that is either determined by the accepted flights before j or a j . If we do not accept j, we move on to optimally assigning the flights j + 1, j + 2, . . . , n and the arguments for j can be applied. Therefore, decisions on whether to accept the flights all occur at integral times and we only need to evaluate g i (t) at the integral times in [0, c] to construct that g i (t). for t ← c to 0 do 6: if t ≤ a i then 7: 
Alternatively, the pricing problem can be formulated as a variant of the shortest path problem with time windows known as the shortest path problem with time windows and time costs (SPPTWTC) in which there is an additional linear cost associated with the service start time at each node. Such a problem is studied in details in [18] . To formulate the pricing problems as a SPPTWTC, a source and a sink need to be introduced and each of the flights in the set F represents a node in the network. The cost associated with each arc between nodes i and j is given by the negative of the corresponding dual variable value, −π ij . To solve the SPPTWTC problem, a dynamic programming algorithm is proposed in [18] . This dynamic programming algorithm is derived based on the general labelling algorithm for the shortest path problems with time windows. For a detailed discussion of the labeling algorithm and some relevant extensions, we refer the readers to [12] and [15] . One of the key reasons for the effectiveness of this dynamic programming algorithm is the use of upper bounds on the service start times to eliminate labels that are infeasible with respect to the time windows, so that the number of labels created is significantly reduced. However, if we reformulate our pricing problems as SPPTWTC, the lack of upper bounds on the park times leads to an exponential number of labels rendering the dynamic programming algorithm proposed in [18] ineffective. Nonetheless, if one wanted to enforce upper bounds on the park times, one would definitely be able to leverage the work in [12] , [15] , and [18] . Therefore, we did not pursue this reformulation further.
The total number of arrivals into a busy international airport can be very large and the proposed dynamic programming algorithm can be slow to obtain the optimal assignments. We present two ways to further decompose the pricing problem of each gate to tackle large-sized instances.
When dealing with a large number of arrival flights, the flights can be usually divided into blocks based on the interactions between the flights to reduce the problem sizes. We introduce the idea of an adjacency parameter, denoted by σ. A formal definition of σ is as follows.
Definition 4.1. For each flight i, let j be the earliest flight after i such that a j > a i + π i + p i , and denote
Note that this adjacency parameter is likely to differ across the gates and change after a new iteration of the master problem is solved with the updated set P k . The use of the adjacency parameter also appears in the work of [5] . Now we are ready to give the block decomposition approximation algorithm.
1: Input: adjacency parameter σ 2: Divide the set F into n/σ +1 blocks such that each of block has σ flights and the last block has n− n/σ ·σ flights.
3: Solve for the optimal assignment in each of the blocks 4: Compute the sum of the total net benefits of the odd-indexed blocks and even-indexed blocks and let B be the collection of blocks with a larger sum of the total net benefits 5: Construct an assignment, S, by cascading the optimal assignments from the blocks in B 6: Compute the total net benefits of the assignment S, f (S) 7 : return S, f (S)
Note that we refer to the discussions in Subsection 4.2 for the computation of the total net benefits f (S).
We now provide an approximation guarantee for Algorithm 4. Now without loss of generality, we assume that the sum of the total net benefits of the odd-indexed blocks is larger than that of the even-indexed blocks, then we have that
Now we show that we can construct a valid assignment by just cascading the optimal assignments from the odd-indexed blocks. Let I 1 , I 2 be any two consecutive odd-indexed blocks and let f 1 be the last flight in block I 1 and f 2 be the first flight in I 2 respectively. From the construction of the blocks, we have that f 2 − f 1 > σ since they are from two different blocks. In addition, by the definition of the adjacency parameter, we have a f2 > a f1 + π f1 + p f1 and, in fact, a f2 > a i + π i + p i for any flight i in block I 1 . Recall that a i + π i is the end point of flight i's acceptance window beyond which i is not accepted. Therefore, the optimal assignment of the flights in the block I 1 does not conflict with whether to accept f 2 in the optimal assignment of the flights in the block I 2 and cascading the optimal assignments of I 1 and I 2 forms a valid assignment. Repeating this argument for all the odd-indexed blocks, we see cascading the optimal assignments of those blocks form a valid assignment pattern.
Therefore, we have that
which completes the proof.
Note that Theorem 4.4 implies that if there exists an assignment pattern with positive total net benefits, the block decomposition approximation algorithm has to return an assignment pattern, S, with a positive OP T S and thus S is a favorable assignment pattern. Furthermore, Algorithm 4 can be improved as follows. Suppose the odd-indexed blocks have a larger sum of the total net benefits and we construct a valid assignment pattern from the odd-indexed blocks. We can still potentially accept more flights from the even-indexed blocks to contribute more net benefits as long as accepting those flights does not violate constraints (3.9) -(3.11).
Lastly, as the adjacency parameter, σ, is computed based on the interactions between flights, it is not likely to be very large. Most flights park at the gate to get ready for the next flight leg with a relatively short turn time and thus they have limited interactions with other flights that arrive hours later.
The rolling horizon framework can also be utilized to decompose the pricing problems to reduce their sizes. As with any implementation of the rolling horizon method, we must decide on two parameters. One is the horizon size and the other is the window size. The horizon size can be either fixed or dynamically determined for each gate and each iteration. The adjacency parameter given in Definition 4.1 is an example of possible dynamically determined horizon sizes. We give a brief implementation of the rolling horizon method in Algorithm 5 that we have implemented in our experiments. Solve for the optimal assignment for the flights {s + 1, s + 2, . . . , s + l} at time t
Fix the assignments of first w flights and add them to S 
The set of feasible patterns P k for all gate k have to be initialized to contain at least one pattern to satisfy the availability constraint (3.3) in the master propblem. Moreover, the union of P k has to satisfy the covering constraint (3.2). Since we start with empty P k , we provide each gate a feasible assignment pattern by randomly assigning each flight to a compatible gate to have an overall feasible assignment.
After optimality is achieved in the linear programming relaxation of the master problem (3.1) -(3.4), we solve the master problem with the binary requirements of the decision variables reinstated to obtain a feasible solution.
We refer to this problem as the binary program. The objective value of the linear programming relaxation provides a lower bound on the optimal objective value while the binary program provides an upper bound on the optimal objective value. Subsequently, we can study the gap between these two bounds. Formally, we compute the gap by
where U B is the upper bound from the binary program and LB is the lower bound from the linear programming relaxation. In the cases where LB is zero, we use the absolute gap of U B −LB which is denoted by an additional notation of "(a)" after the numericals. We will use the abbreviations of UB and LB to represent the upper bound from the binary program and lower bound from the linear programming relaxation respectively in the computational experiments section. Table 1 . We test multiple instances of the same size in the moderate and large size categories to see if there are variations in the performances.
For the implementation, the pricing problem algorithms are implemented in Python. Gurobi (version 9.0) is used whenever a commercial solver is needed. The computations are performed on an Unix system with a 12-core CPU and 16GB RAM. We also set a limit of 3 hours for each instance.
For a complete comparison, the performance of the compact mixed integer programming (MIP) formulation (2.1) -(2.8) is considered using the small-sized instances and shown in Table 2 . The formulation is solved using the Gurobi solver. As we have noted before, the MIP formulation is not ideal for this problem as evidenced in Table 2 . It is important to note that we can combine the pricing problem algorithms presented in Section 4 in our implementation. For the small-sized instances, we test the following ways: The detailed computational results for the small-sized instances are shown in Table 3 and Figure 2 . We see that our proposed approaches out-perform the Gurobi solver in all six instances. While the approximation algorithms are designed to provide additional improvements relative to the dynamic programming algorithm, however, they do not offer any substantial boost. In general, if we increase the size of the instance, it becomes more difficult to solve. We observe that the flight to gate ratio has an impact on the time. In particular, although instance 1 is smaller in size than instances 2 and 4 are, it takes longer than both instances 2 and 4 do to be solved. The flight to gate ratio determines how congested the gates are. The higher the ratio, the larger the number of flights each gate on average has to accept. Instances with higher ratio require more delicate assignments. Furthermore, in contrast to the worse case theoretical analysis, the dynamic programming algorithm performs well. This is likely due to the fact that the inter-arrival times are much shorter than the processing times and consequently many recursive calls to the evaluation oracle have input t beyond the corresponding flights' acceptance windows. Therefore, the dynamic programming algorithm effectively only considers a smaller subset of F .
Additional experiments were performed to understand if there exists an optimal rule that is independent of the instances about the number of iterations of approximations to use. We take instances 1 and 3 and vary the number of iterations of approximations used and compare the time taken to solve the instances. The plots are shown in Figure 3 . In both plots, the dotted grey line gives a baseline for comparison. It seems such a rule does not exist for either approximation algorithm. The optimal number of iterations is highly dependent on the instance.
Next we show the computational results for the moderate-sized instances. As we have seen that the solver option is not as effective as other options are in solving the small-sized instances, only the options of DP, SM+DP, and ADP+DP are considered in this set of experiments.
The results are shown in Table 4 . We also include part of the results from the small-sized instances here for much less time compared to each iteration of the dynamic programming algorithm. We also observe that the submodular maximization out-performs the approximative dynamic programming in most of the cases as the assignment patterns that are generated by the submodular maximization are usually of better total net benefits than those that are generated by the approximative dynamic programming. We believe that the submodular maximization is more total net benefits oriented as the probability of accepting a flight is computed directly based on the change in the total net benefits accepting and not accepting the flight. Moreover, we observe small variations in the times taken across instances of the same size and this suggests that besides the sizes of the set of flights and the set of gates, other factors associated with the flights and the gates can complicate the problem.
Following the experiments on the moderate-sized instances, we move on to the large-sized instances. As we discussed in the Subsection 4.4, we can utilize the block decomposition approximation and rolling horizon framework to tackle the large-sized instances. The horizon size and window size are important parameters in the rolling horizon framework. We first perform parametric studies to understand how the horizon size and For the large-sized instances, we again have a few possible ways to solve the pricing problems because we can either use a fixed or a dynamically determined horizon size for the rolling horizon method and also utilize the approximation algorithms. Here is a detailed breakdown: instances. It can potentially improve the performances and serve as a benchmark for the block decomposition approximation. Note that whenever an algorithm is needed to evaluate the optimal assignment during the rolling horizon process or the block decomposition approximation, the direct implementation of the dynamic programming algorithm (Algorithm 2) is used. We also initialize the set P k to be identical for the comparisons.
The results of the computations are shown in Table 5 . Note that the last column in the SM+RHM option reports the constant threshold c used (option 4 above). After extensive experiments, we observe that the adjacency parameters, σ, can be very large at beginning and, consequently, the options of RHD and BD+RHD that involve decompositions based on these parameters become very ineffective and unable to obtain a solution within the allocated time limit. Therefore, we do not report the results of these two options in Table 5 . From the results, we see the rolling horizon method with a fixed horizon size performs much better than the same method with a horizon size of σ which struggles to solve these large-sized instances. If we further keep the horizon size at most 10/20 in the RHM option, we obtain better performances in all instances compared to the rolling horizon method with a fixed horizon size. In addition, submodular maximization algorithm efficiently provides assignment patterns of good quality when σ is large and improves the performances compared to both RHF and RHM options in all instances. Lastly, like what we get for the moderate-sized instances, we still see some variations in the performances across the instances of the same size. We next show how the adjacency parameters vary with the number of iterations. Figure 4 plots the adjacency parameters for each iteration when instance 15 is solved by the rolling horizon method with a fixed horizon size of 10 and window size of 1. The adjacency parameters during the early iterations have a larger variance than they do during the later iterations and the average of the adjacency parameters for each iteration also displays a decreasing trend and eventually reaches zero. 
In summary, it seems that the dynamic programming algorithm together with the submodular maximization offer the best performance for the small-sized and moderate-sized problems. For the large-sized instances, it seems that running the submodular maximization algorithm when σ is large and using the rolling horizon method while keeping the horizon sizes at most 10/20 otherwise lead to the best performances. Moreover, the binary program produces feasible solutions of good quality for almost all instances as the gaps are usually very small. The amount of time taken to obtain those feasible solutions is also reasonable in practice.
Finally note that all instance data can be accessed online at the authors' websites.
We use a column generation approach to solve the flight-to-gate assignment problem aimed at minimizing the total arrival delays. Instead of using the integer program solver for the pricing problem, more efficient approximation algorithms and dynamic programming algorithms are proposed. In addition, the rolling horizon framework is used together with approximation algorithms to tackle the instances with larger size. Multiple computational experiments are performed to validate the approach and show significant improvements over the integer program solver.
There are few possible extensions to consider. One of them is to consider the interferences between different gates. With different airport layouts, this can be realized by linking constraints in the master problem or a decomposition of the pricing problem by gate groups instead of by each individual gate. Another extension is to take the uncertainties in the arrival times and turn times into consideration. A stochastic programming approach may be suitable to tackle this version of the problem.
This proof is given in [4] and we modified it to match our assumption.
Let OP T denote an optimal solution and OP T i := (OP T ∪ X i ) ∩ Y i , then we have that OP T 0 = OP T and OP T n = X n = Y n . Here are two useful lemmas to be used later in the proof of the theorem. 
Lemma 7.2. For every 1 ≤ i ≤ n,
Proof. Proof of Lemma 7.2. It is sufficient to prove (7.3) conditioned on any event of the form X i−1 = S i−1 , when S i−1 = {u 1 , · · · , u i−1 } and the probability X i−1 = S i−1 is non-zero. Hence fix such an event for a given S i−1 . The rest of the proof implicitly assumes everything is conditioned on this event. Note due to conditioning, the following quantities become constants,
• a i and b i . From Lemma 7.1, we have a i + b i ≥ 0, so we only need to consider three cases.
1. a i ≥ 0 and b i < 0. In this case, a i /(a i + b i ) = 1. Then we have Y i = Y i−1 and
Then we are left to prove that
If u i ∈ OP T i−1 , then the left hand side is zero and this inequality is satisfied. If u i / ∈ OP T i−1 , then we note that
Next, by the definition of submodularity of f , we have now
2. a i < 0 and b i ≥ 0. This case is analogous to the previous case.
3. a i ≥ 0 and b i > 0. In this case, we have a i = a i and b i = b i . Then we can compute the left hand side of (7.3) by
Last inequality follows from the following cases. Note that u i ∈ Y i−1 and u i / ∈ X i−1 , (a) u i / ∈ OP T i−1 , then we note that the second term of the left hand side the last inequality is zero and
Next, by the definition of submodularity of f , we have now
(b) u i ∈ OP T i−1 , then we note that the first term of the left hand side the last inequality is zero and
Next, by the definition of submodularity of f , we have now
Then with (7.8) and (7.10), we are left to verify
which can be easily verified.
With the two Lemmas, we are ready to show the approximation guarantee of the algorithm.
Proof. Proof of Theorem 4.1. If we sum all inequalities from lemma 7.2 over 1 ≤ i ≤ n, we have that
]. (7.17) Note that X 0 = ∅, so we have f (X 0 ) = 0 and by our assumption we have f (Y 0 ) = f (F ) ≥ 0. In addition, we have OP T 0 = OP T and OP T n = X n = Y n , so (7.17) gives that