Abstract

After a disruption to its operations, an airline needs to dispatch its aircraft and accommodate affected passengers in order to resume normal operations. When facing a flight delay or cancelation, passengers usually have two options—switching to a different itinerary or receiving ticket refunds. This paper focuses on the integrated recovery of both aircraft and passengers by taking into consideration passengers’ preferences regarding the two options. Objectives include minimizing the financial loss of airlines and minimizing the utility loss of passengers. To solve the optimization problem, we present a loop-based multiobjective genetic algorithm (GA) by leveraging a special characteristic of flight networks. Experiments based on real-world data demonstrate the effectiveness and stability of the algorithm in various situations. The outcome has theoretical and practical implications for disruption management and airline operations.

1. Introduction

According to statistics from the International Air Transport Association (IATA), the global civil aviation industry transported 4.1 billion passengers in 2017, and its growth rate was 7.4%, which is the greatest one in history. However, the total punctuality rate in 2017 was 76.35% and this rate varies significantly among different airlines, which claims that different operational strategies may have great impacts on the punctuality rate. In China, with development of the Chinese economy, civil aviation transportation has unique advantages in convenience and time. With the support of various policies, the number of airline passengers of different countries as well as cities in China has rapidly grown. At the same time, however, complaints from civil aviation passengers have also risen sharply. According to statistics from the China Civil Aviation Administration (CAAC), the total number of complaints received from passengers reached more than twenty thousand in 2017, nearly twice as many as 2016, and flight delay problems accounted for 53% of the complaints. Why did this happen?

One of the most important reasons for the low customer satisfaction with airlines’ handling of disruptions is that airlines often prioritize their own operational convenience and economic impact over passengers’ preferences. Such prioritization of their own margins during flight rescheduling is very common for profit-seeking business like airlines. However, it is extremely unfair and intolerable for their disrupted passengers. The incident of a passenger dragged off a United flight, posted by CNN wire on Apr. 11th 2017, highlighted such ignorance of passengers’ preferences. Meanwhile, the core competitiveness of an airline is its ability to satisfy passengers’ needs. Chang et al. [1] evaluated airline competitiveness and determined that customer service quality is one of the primary considerations of passengers. Thus, during the postdisruption recovery process, airlines should pay more attention to passengers’ preferences when deciding how to accommodate them.

Common methods for airlines flight scheduling recovery from disruptions include flight delays and cancellation, aircraft swapping and ferrying, crew repair, and passenger reassigning; these have been considered in various prior studies, such as Clause et al. [2], Le et al. [3], Castro et al. [4], and Artigues et al. [5]. A topic of concern to an increasing number of researchers is integrated recovery of multiple resources, such as aircraft and crews. In contrast, recovery of passengers’ itineraries, especially their preferences regarding alternate itineraries, has gone unheeded by most studies. However, during the practical recovery process, there is a close connection between passenger recovery and aircraft recovery, and passengers’ satisfaction is paid significant attention by most airlines because of its potential impact on their own reputation and competition level in the aviation market. Therefore, it is critical to integrate interests of an airline and its passengers into a single recovery model, especially with consideration of passengers’ preferences regarding alternate itineraries.

In practice, the cost of recovery for airlines is still an important aspect for the operation of an airline so that a multicriteria model should be built to reflect these two considerations. In the field of integrated of flight, multicriteria model is a common modeling strategy, such as Lee et al. [6], Liu et al. [7], and Khaled [8]. The consideration of different kind of objects for the recovery and using functions separately to describe these objects is the key of this modeling style. After that, well-designed mechanisms are introduced to balance these objects. Although the most common way to handle multiobjective problems is to set weights of different objects so that they can be transferred into one single object, which is relatively easy to solve, we still try to keep two independent objective functions so that both objects can be fully considered. Moreover, the separate consideration of two objects avoid the complexity of how the weight of each object is set during the solving process. Thence, aiming to model the integrated postdisruption recoveries of aircraft and passengers, a multicriteria model is necessary.

The goal of this paper is to define, formulate, and efficiently solve the problem of satisfying passenger preferences in the integrated recovery of both aircraft and passengers. To achieve this goal, it was not only necessary to develop and extend mathematical formulation to consider the disrupted passengers’ preferences in the integrated recovery process but also to design an efficient solution method for the integrated recovery problem with consideration of passengers’ willingness. This paper provides an efficient and computationally manageable method for integrating the aircraft and passenger recovery problems considering passengers’ preferences by establishing a multiobjective optimization model formulation. The objectives of the formulation include clearly describing the airline recovery cost and passengers’ interest loss during disruption recovery. Additionally, a loop-based multiple objective genetic algorithm (GA) that leverages the characteristics of flight networks is also designed to obtain Pareto-optimal solutions for this problem extremely efficiently.

The remainder of the paper is organized as follows: after a review of relevant studies of airlines’ disruption management in Section 1, the problem and model are formally described in Section 2. Section 3 presents a heuristic method for the integrated recovery problem with passenger preferences. Computational experiments that evaluate our approach are presented in Section 4. Section 5 summarizes our work and discusses directions for future work.

2. Literature Review and Contributions

2.1. Postdisruption Recoveries of Aircraft and Passengers

Airline disruption recovery, which is an important guarantee of airline daily operation, includes total aircraft, crew recovery, and passenger recovery such as Clause et al. [2]. Moreover, Hu et al. [9] think passenger itinerary recovery is the most important measure to determine the quality of airline service. The initial research topic focuses on aircraft recovery. Teodorovic et al. [37 performed the first aircraft recovery study using an optimization method. Subsequently, many scholars have primarily focused on the problem of rearranging aircraft after disruptions using minimum–cost network flow theory and corresponding algorithms such as Gershkoff [11], Jarrah et al. [12], Yan et al. [13], and Thengvall et al. [14]. A time-band network was constructed by Bard et al. [15] for the aircraft rerouting optimization problem; it is often combined with the Column Generation algorithm and Bender Decomposition algorithm, to generate feasible aircraft routes in the airline disruption recovery problem by Eggenberg et al. [16]. For the large-scale flight delay recovery problem, metaheuristic algorithms are preferred to obtain near-optimal solutions efficiently, such as Argüello et al. [17] and Løve et al. [18]. However, most studies only focus on single objectives concerning the cost of deviations from original flight schedules. Only a few attempts, such as Liu et al. [19], have employed a hybrid multiobjective GA to find solutions for a daily short-haul aircraft schedule recovery problem.

Besides recovering aircraft, it is also important to recover several resources simultaneously in airline disruption management. Since integrating recovery of the total resources is a difficult task due to the size of the resultant problem, only a few attempts, such as Lettovsky [20], Peterse et al. [21], and Arıkan et al. (2013), have appeared in the literature.

Most studies try to describe the integrated recovery problem of aircraft and passengers in a single model formulation. Bratu et al. [22] were the first to model passenger recovery using a passenger delay metric model (PDM) and a disrupted passenger metric model (DPM). Both models consider the passengers’ disrupted cost and airlines’ operation cost in a single objective. The latter one is validated in three different levels of scenarios and three different degrees of disruption. Jafari et al. [23] presented a single-objective model that simultaneously recovered aircraft and passengers. Their dataset contained 13 aircraft of 2 fleets, and they were able to address disruptions in a small-scale airline.

Bisaillon et al. [24] developed a large neighborhood search heuristic for an airline recovery problem combining aircraft routing and passenger reassignment with one objective of minimizing operating costs and impacts on passengers. Sinclair et al. [25] improved the heuristic by adding some additional steps in each phase for the same problem. Then, Sinclair et al. [26] solved the similar problem with a column generation postoptimization heuristic, which slightly increased the allotted computing time. The heuristic is based on a mixed-integer programming mode and forms various hierarchies of passengers, flights, and other elements to recursively solve the problem. Zhang et al. [27] proposed a three-stage sequential math-heuristic framework to solve the integrated airline service recovery problem. Time-space network flow representations, mixed-integer programming formulations, and algorithms that take advantage of the underlying problem structures were proposed for each of the three stages. However, the computational results of the above studies were employed for the 2009 ROADEF Challenge that was organized by Artigues et al. [5] instead of for real-life situations.

Hu et al. [28] and Hu et al. [9] aimed to find the optimal trade-off between passenger delay costs, passenger reassignment costs, and the cost of refunding tickets in a single-objective model for the integrated recovery problem of both aircraft and passengers. They both consider the airline’s point of view instead of those of the passengers. The former solved the model using CPLEX Solver, and the latter design a meta-heuristic based on GRASP to obtain the near-optimal solutions. In the former research, the average time for solving a real example, including 319 flights and 59 airports, is 49 seconds and the average optimization efficiency is about 10%. Therefore, the computational time and optimization efficiency still have chances to improve.

From the literature review analysis, we observe that most of these papers combine utility losses of passengers with airlines’ financial losses to form a total loss as the single objective of the optimization model. This combined object simplifies the problem-solving process with a certain degree of relevance. However, unlike aircraft or crew members, passengers do not belong to the airline. They are two juxtaposed participants during flight production. Faced with the unforeseen disruption and the following recovery process, an airline cannot cover and compensate for the total loss of disrupted passengers, especially the utility cost and dissatisfaction of the passengers with the airline’s service quality. Therefore, the disruption and recovery costs of the airline and passengers cannot be simply combined. They should be subdivided in order to obtain a recovery solution that is more suitable for both participants in long-term plans.

2.2. Passengers Preferences

In recent years, research on passengers’ travel preferences has gained popularity. Most studies use a logit model to analyze the choice of passengers between flights or railways by Moeckel et al. [29], revealing factors influencing passengers’ flight choices by Algers et al. [30], Yan et al. [31], Hagmann et al. [32], and Fleischer et al. [33], predicting the total number of air travelers [34], or optimizing passenger flows in a flight network, Dou et al. [35] and Yang et al.(2017). The above models are only suitable for normal transportation networks, not analysis of passengers’ preferences in the case of an accident or disruptions. Only a few attempts have discussed the possible behavior of railway passengers during emergency evacuation, and the analysis concluded that gender may result in significant difference in evacuation behavior in the ordered logit model.

However, although an airline’s handling of disruptions is a major source of its passengers’ dissatisfaction, Bratu et al. [22], no study on airline postdisruption recovery has attempted to incorporate passenger preferences. This paper tries to address this gap by quantifying passenger preference via utility losses and considering such losses along with airlines’ financial losses in a multiobjective optimization model. To solve such a multiobjective optimization problem, we choose to use GAs.

2.3. Multiobjective Model and Optimization

In most studies, as stated above, traditional airline recovery problems are only described as single-objective models. Even when faced with multiobjective problems, they often propagate a linear sum on objectives to change them into a single one in the mathematical formulations. However, it is stated that the single-objective solution cannot overall reflect the true interests of both participants, airlines and passengers, Fieldsend et al. [36]. Therefore, a multiobjective model is naturally in-need to describe these two considerations. On one hand, the introduction of multiobjective model is used on the modeling, which means different kinds of object functions are built. Liu et al. [7] used a five-object model to manage a disruption and different weights were set to continue the solving process. On the other hand, the multiobjective idea can be used in the process of solving, such as Burke et al. [37], Chou et al. [38], and Lee et al. [6]. However, there are very few studies that consider the multiobjective optimization problems of integrated recovery of aircraft and passengers. Liu et al. [19] 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 inequality-based multiobjective GA. A simulated disturbance experiment involving daily domestic airline plans in Taiwan with only 7 aircraft and 39 flights was performed. The efficiency was relatively low. Moreover, the objectives lack formal descriptions of flight cancellation and passenger reassignment, which are common in most airline recovery situations. The short-haul recovery approach is not suitable for general-scale disruptions, especially in larger airlines.

2.4. Genetic Algorithm

GAs represent a widely used method for intricate multiobjective optimization, especially for airlines’ integrated recovery problems. This method, based on a Pareto noninferior solution, is a powerful multiobjective tool, Goldberg [39], that exhibits the superiority and inferiority of each solution. These solutions are then classified to choose the parent individuals during the solution process, leading the whole population to the forefront of Pareto solutions. GAs are especially suitable for problems with large feasible regions and complex constraints.

Especially in the research field of flight rescheduling, the aircraft flow being continuous in time and space increases the difficulty of the integrated recovery problem with multiple objectives, Hu et al. [40]. It is also the reason why few studies have investigated multiobjective optimization for airline recovery. However, the characters of flight loops, especially in the Chinese flight network, can package and protect the aircraft flow as continuous in time and space. Therefore, a loop-based multiple objective GA can be designed by combining an efficient coding strategy and the loop characteristics of the flight networks, in order to obtain near-optimal solutions of the integrated recovery problem.

2.5. Contributions

The main contributions of our research are as follows.

(1) An integer programming model with two objectives constructed in this study reflects the practical interests of two major participants: airline and the disrupted passengers during the disruption recovery process. One objective refers to reducing the financial costs of airlines, and the other one focuses on reducing passengers’ utility losses.

(2) Different from the previous research topic, passengers’ preferences between refunding or endorsing are first considered in the integrated recovery model formulation of aircraft and passengers in order to improve passengers’ satisfaction. It is the first attempt in which passengers’ psychological behaviors are quantified in a probabilistic form in the integrated recovery optimization model.

(3) It is the first time that a recovery solution approach for resolving multiobjective airline rescheduling problems that combines a characteristic flight loop network and a GA is developed. The combination results in distinctly improved efficiency and better solution performance for the integrated airline recovery problem.

3. Problem Description

3.1. Passenger Preferences Analysis

An airline must arrange its aircraft to cover its flights each day. In the case of inclement weather, emergency repairs, and aviation control, the availability of aircraft could change, causing disruptions to the airline’s operations. The airport operation control center (AOCC) will collect all information and start the recovery period. A typical recovery period involves rearranging the assignments of aircraft to flights, mainly with the goal to reduce the financial loss for airlines. However, such a recovery process often ignores the preferences of passengers whose itineraries are affected by the disruption. In other words, it is assumed that passengers will accept the new itinerary provided by the airline.

Figure 1 shows the decision process of a passenger who belongs to a disrupted flight. If the flight is not canceled, which means the flight will land off later, we assume this passenger will still wait for the delayed flight. However, if the flight is canceled, the problem for the passenger of this flight is whether to apply for refund. If so, this passenger will receive a refund from the airline and his/her journey will end. If not, this passenger insists his/her journey and wants to move to another flight. In this paper, a proportion factor is used to illustrate how many passengers want to choose endorsement when the cancellation of flight happened. This factor is external and based on the operation information of the airline. In this situation, the problem waves back to the airline: are there enough available seats for the passenger who wants to endorse? The answer relies on the recovery decision the airline made. If the passenger’s demand can be satisfied, he/she can continue the journey. If not, the passenger will be annoyed about the recovering process so that the utility loss of passenger increases. In conclusion, in the whole recovering process, passengers only need to show their preference to which extend they choose refunding and this is the information they need to give to the airline when adjusted flight schedule is announced.

According to the statement above, in the recovery process, although passengers do not participate in creating a new flight schedule, they can pick between two options. One is endorsing, which means the passenger chooses to transfer to another flight with the same airline. The other one is ticket refunding, which means passenger cancels their current itinerary with this airline, and the airline does not have to provide any extra transformation service for them but rather helps them apply for a refund. In this paper, passengers do not passively accept the adjustment but actively choose to endorse or refund when flight cancellations occur. The airline reassigns the aircraft based on the passengers’ willingness during the recovery period. For the endorsing choice, some passengers’ willingness cannot be satisfied due to the limited available aircraft capabilities and other reasons, which causes a decline in satisfaction. Thus, the utility loss is used to describe passenger dissatisfaction with the airlines. For ticket refunding choice of passengers, the service task seems easy for the airline; however, due to the failure of the service, the airline has to absorb both the economic cost and reputation cost due to the disruption.

Considering that the utility loss in the delay situation has a high degree of consistency with the economic loss, which may produce multiple collinearities, and that a significant amount of computation is needed to calculate the endorsed passengers’ utility loss, this research considers only the utility loss of the passengers’ willingness to endorse but cannot be satisfied when a cancellation occurs.

3.2. Recovery Operation Hypothesis

The options of integrated recovery in this paper are as follows: cancel—no aircraft will be arranged for this flight, and the passengers will be reassigned according to their willingness and aircraft capacity; naturally delay the same aircraft—a new timetable of flights will be provided according to the available times of the delayed aircraft; swapping aircraft—the matches between aircraft and flights will be rearranged.

During recovery operations, it is assumed that the crew requirements are not considered, so modifications to aircraft routes are not limited by crew availability. When an aircraft route is selected, it is acceptable as long as appropriate crew members can be found to fly the aircraft. Because the time window of the recovery is only a single day, the arrival and departure times of the involved flight must be the same day to prevent interference with the flight schedule of the following day. The time length of a single day does not exceed the maximal duty time of most crews according to the Airline Operation Manual.

In reality, a fleet is widely used in civil aviation to describe a type of aircraft, such as the Boeing 787 fleet. The main reason for this use is that within the same fleet, the crew members, maintenance team, and other resources can be shared. Some flights can be flown only by a specific type of aircraft due to distance limitations, arrival airport situations, and other reasons. For example, remote aircraft are usually used for intercontinental flights. Therefore, in this paper, we do consider multiple fleet types and the fact that each fleet has different characteristics that limit cross-assignments. Substitution for flights is allowed between some fleets in accordance with airline policies and aircraft configurations.

When flight interference occurs, passengers’ itineraries are forced to be adjusted. We claim that only the flights are cancelled, and the corresponding passengers can choose new itineraries. Otherwise, all passengers should travel on their original flights. The passengers of the cancelled flight make their decisions based on their own characteristics. As described above, faced with flight disruption, passengers have two options: endorsing or ticket refunding. Endorsing is the first choice of most passengers except when there are not enough seats to accommodate them. So, some assumption about passengers’ endorsing should be presented clearly. Since different fleets have different passenger carrying capacities, this paper considers only aircraft that belong to the same fleet for passenger exchanges. In this circumstance, the class and the capacity are the same, which provides more feasibility for airlines and passengers. Additionally, it is assumed that a passenger has only one flight in an itinerary due to the low connecting flight rates and that most domestic flights will not change flight numbers during a voyage, which has substantial equivalence to a single flight. This paper focuses on Chinese domestic flights, for which the majority of the passengers can complete their trip in one flight.

3.3. Problem Objectives

The problem raised in this paper reflects the practical operation of an airline when accidents force its operating plan to be altered. When some aircraft are delayed or grounded due to unforeseen events, the flights covered by these aircraft and passengers’ itineraries will be disrupted. To recover the schedule and passengers’ itineraries as soon as possible, available aircraft need to be rerouted and passengers reassigned according to their willingness. Several objectives have been established to guide the recovery and evaluate the total losses of the airlines. In this paper, the objectives include the following two aspects.

(1) Economic loss. The objective is to minimize the airline financial cost, including compensation for passenger delays, the transiting cost of transferring them to other flights, and ticket refunding cost when the passengers cannot be reassigned. The flight delay cost refers to the direct cost due to flight delays, such as additional fuel cost, overtime payment, and additional meal cost for each passenger. It is apparently associated not only with delay time but also with the number of passengers. The transiting cost is the service cost for each passenger transiting from a cancelled flight to the other flight. It is related to the difference between the actual departure time of the new flight and the scheduled departure time of the cancelled flight. The refunding cost mainly refers to airline compensation for the passenger original ticket payment due to the airline service discontinuing.

(2) Utility loss. The utility cost is mainly to describe the dissatisfaction of passengers on airline recovery operations and their potential economic loss. Faced with the itinerary disruption, most passengers prefer endorsing. However, whether endorsement is feasible depends on the capacity of the aircraft, which covers the passengers’ new itinerary. If the new arrangement does not have sufficient capacity to accommodate all passengers who endorse, some of the passengers will need to have their tickets refunded, which will generate a utility cost for these passengers. Compared with passengers who prefer refunding, the passengers who want to be endorsed but cannot be satisfied need to be emphasized because they suffer more from the cancellation; thus, the utility loss of this type of passenger is relatively high. The more dissatisfied the passengers become with the airline, the less likely they are to choose airline’s flights in the future, which can generate potential economic loss and loss of market share for the airline.

3.4. Problem Constraints

As the integrated recovery problem in this paper must be solved in practice, the flight rescheduling, aircraft rerouting, and passenger recovery must also obey the following constraints.

First, each flight can only be adjusted once or cancelled in the feasible aircraft route, its adjusted departure time must not be earlier than the original time, and its delay time must be less than the maximum delayed time. Moreover, for every couple of neighborhood flights in one flight sequence, the former flight’s arrival airport must be the same as the latter flight’s departure airport, each neighborhood flight’s interval time must be smaller than the minimum turnaround time (the minimum time for the aircraft to rearrange for the next flight after landing), and the arrival time of the last flight must be earlier than the curfew time.

Each aircraft can cover at most one route to form one feasible aircraft route that remains continuous in the time and space connection network. Additionally, no aircraft scheduled for maintenance service during the recovery period will have its original route altered. Finally, at the end of the recovery period, there should be a sufficient number of appropriate aircraft types positioned at each airport to ensure that the published schedule beyond the recovery period can be executed to avoid disruption continuation into the following day.

For every endorsement willingness, the available flight’s departure time must be later than the departure time of the cancelled flight, and the aircraft type between two switched flights must be the same. Moreover, the routes must be the same between two switched flights, i.e., passengers who choose to endorse should be promised to arrive at their original destination. Finally, the capacity of the aircraft must allow for the extra passengers; otherwise, they will have to refund their tickets.

3.5. Mathematical Formulation

To accurately describe the problem, an integer programming formulation is constructed. This model aims to allocate limited aircraft to satisfy the needs of all types of flights and the constraints of aircraft flow and airport curfew and to reallocate the passenger into available flights with sufficient capacity according to their desire.

Index: aircraft fleet index;: aircraft route index;: aircraft index;: flight index;: airport index.

Sets: set of aircraft fleet;: set of aircraft routes;: set of aircraft;: set of flights;: set of airports;: set of aircraft that belong to fleet ; ;: set of aircraft routes that include flight ; ;: set of aircraft routes that end in airport ; ;: set of aircraft routes that can be covered by aircraft ; ;: set of aircraft routes that can be covered by aircraft fleet ; ;: set of flights that can receive passengers from flight ; ;: set of flights that can send passengers to flight ; .

Parameters: cost of each passenger who belongs to the delay flight covered by route and aircraft ; , , ;: cost of each passenger who belongs to the cancelled flight ; ;: cost of each passenger who belongs to the cancelled flight and transfers to a flight covered by the route and aircraft;: delay time of flight covered by route and aircraft ;: the maximum amount of time that a flight can be delayed;: number of aircraft that belong to fleet and should terminate at airport ; , ;: number of passengers who are initially in flight ; ;: number of seats in aircraft fleet ; ;: the rate of passengers who want to endorse and belong to the cancelled flight ; ;: the coefficient of utility loss of each passenger who belongs to the cancelled flight and want to transfer but cannot be accommodated; .

Decision Variables: if route is covered by one aircraft; if not, ; ;: the actual number of passengers who transfer from flight to a flight covered by route and aircraft; , , ;: the actual number of passengers who refund tickets and belong to flight .

Main Formulas.t.

The objective of function (1) is to minimize the sum of the delayed cost, the transferring cost, and the refunding cost. The objective of function (2) is to minimize the sum of utility losses. Constraint (3) ensures that every flight is either cancelled or covered by a flight. Constraint (4) ensures that a sufficient number of aircraft are arranged for the needed airport prior to the curfew time. Constraint (5) ensures that every route should be covered by at most one aircraft. Constraint (6) ensures that the maximum delay time satisfies the demand. Constraint (7) ensures the balance of the transferring stream. Constraint (8) ensures that the adjustment satisfies the volume of aircraft. Constraint (9) emphasizes that some desires of passengers cannot be satisfied. Constraints (10) - (12) define the decision variables as integer or binary.

The proposed model (1)-(12) is an integer programming model based on a set covering model formulation with side constraints about flights and passenger. In some studies, an integrated recovery model of aircraft and passengers has been directly solved using business software, such as CPLEX or LINGO [28] Arıkan et al. (2013) and Arıkan et al. (2016). However, considering the exponential relationship between the aircraft route and aircraft and the extreme complexity of passenger transiting process with their willingness, plenty of time and storage must be used for business software to generate all of possible transiting relationship and all possible aircraft routes. Some other studies have introduced the column-generation algorithm and bender composition algorithm to solve this problem. However, the biggest resistance for the above algorithm is that there are two objectives in the proposed model (1)-(12), and the complicated passenger transiting relationship between flights may negatively affect the efficiency of the above algorithm.

In contrast, faced with this complicated biobjective integrated recovery problem, it is critical to design an efficient heuristic algorithm to obtain a series of nondominant recovery solutions that meet the selection requirements of the AOCC staff according to the real-time environment. In this paper, we provide a recovery solution method with more competitive performance and more efficiency by designing a heuristic based on a multiple objective GA and flight loop network characteristics. Detailed information about the proposed solution method is described in Section 4.

4. Solution Method

Under the consideration of the proposed biobjective model (1)-(12), referring to the characteristics of flights in China and the rules for actual flight operations, a novel recovery solution approach (denoted as LBMGA) is designed, based on a multiple objective GA and flight loop network characteristics, to solve the problem efficiently. The significant efficiency is mainly derived from excellent parallel searching ability of GAs in complex conditions and excellent structural compatibility with a multiobjective mathematical model of the proposed approach. Additionally, since the clash between the mechanical application of a general GA processing strategy and the space continuity characteristics of the flight sequence may result in a negative effect on the approach efficiency, a clever encoding method is designed according to loop features of the flight network, which can perfectly avoid the negative effect without blemishes. Moreover, the path of the searching for better solutions between two objects is also significant for the total performance of the algorithm. Thus, we also chiefly introduce the multiple objective choosing strategy in the following subsection.

Overall, the subsections are organized as follows. Based on graphical analysis of the flight loop features in Section 4.1, a novel coding strategy of the GA is presented in Section 4.2, and then the description of the fitness calculation method follows immediately. Finally, the choosing strategy of multiple objectives of the proposed model in the recovery solution method is introduced in detail.

4.1. Feature Analysis of Flight Schedules

Flight schedules indicate how flights are organized by an airline. In this paper, we analyze current practical flight schedules and find a loop-based flight network character. In the flight network, one flight sequence connects each air spoke via air hubs; no direct connections exist among air spokes. Simultaneously, all available airports are separated into groups; each group contains one air hub and some air spokes. Thus, each air hub only services the same group of air spokes. In most cases, an aircraft that flies from an air hub to an air spoke must return to the same air hub to create a flight loop (FL). In other cases, the aircraft flies from an air hub to another air hub and then immediately returns to the original hub; this situation also makes an FL.

An FL consists of two adjacent flights performed by the same aircraft, including a contrast origin, a terminal point, and at least one air hub.

Figure 2 and Table 1 show a small example of this special system. A and B are air hubs, and the remaining letters represent air spokes. Flight ① and flight ② constitute a flight loop, and flight ③ and flight ④ constitute a flight loop. If a flight is b1→b2, the passengers of the flight need to fly from b1 to B and then from B to b2. However, if a flight is b1→a1, the passengers need to transit by B and A to reach the destination.

Referring to some existing flight schedules in China, numerous Chinese airways choose this system to arrange their routine. The most significant characteristic of this special system is that the percentage of FL is at an extremely and remarkably high level.

4.2. Gene Coding

The traditional theory of a chromosome is to encode a partial schedule of all flights. If a sample consists of n flights and m aircraft, the length of the chromosome is n, which changes from 0 to m for each gene point. For example, the gene sequence (2, 0, 3) indicates that the first flight is performed by aircraft 2, the second flight is cancelled, and the third flight is performed by aircraft 3. This strategy is intuitionistic, and the algorithm can be easily realized. Concerning massive restrictions that need to be obeyed during the recovery problem, especially the continuity constraint, the traditional choice will be very inefficient and slow.

In this paper, a segmented encoding theory based on FLs is proposed to improve efficiency and accelerate the problem-solving processes. According to the features of an FL, if a recovery case contains n air hubs, the chromosome will be divided into n sections, and the length of each section is determined by the number of specific sets of FLs of the air hub. Thus, each gene point represents a specific FL. Using the information in Figure 2 as an example, supposed the flights belonging to air hub A is A→a1,a2,a3; a1,a2,a3 →A so there are three flight loops for air hub A and the corresponding aircraft is shown in Table 2. It means that aircraft 1 flights the loops (A→a1,a1→A) and (A→a2,a2→A). The information about air hub B is also shown in Table 2. Therefore, the gene sequence of this example is (1,1,2) (3,4,4). Particularly, if the number of aircraft is 0, it means there is no aircraft is arranged to the flight loop so that both flights of the loop will be canceled

The major concern of this coding theory is to cope with the continuity constraint. Using an FL, every crossover, mutation, and other operations for each gene section will not break the continuity constraint to reduce the number of calculations. In this special strategy, only the FL can be involved in the recovery process. Because the percentage of FLs is very high in the Chinese flight network, the optimization efficiency is guaranteed.

4.3. Fitness Calculation

Fitness one aims to calculate the total economic cost for each flight schedule. Using the previous method, when every flight is confirmed to be performed by an aircraft, the flight schedule and the cancel table is confirmable. By matching every original departure time with the present departure time and matching every original arrival time with the present arrival time of each flight, the delay time of every noncancelled flight can be calculated to obtain the delay cost. By adding the cancellation cost for the cancelled fights, the total economic cost can be calculated.

Fitness two aims to calculate the utility losses. Every flight schedule corresponds to a cancellation table. For every cancelled flight f, there exists a set of flights () that can receive passengers from flight . is a dynamic set. For every flight , the number of available seats () has been given as background information. represents the number of empty seats after an airline stops selling tickets for flight . Thus, for every cancelled flight , the total number of available endorsing seats can be calculated based on a specific flight schedule.

Flight fe needs to obey the available constraint, the airport constraint, the time constraint, and the aircraft type constraint. The airport and the aircraft type of flight fe is relatively fixed, whereas the available situation, the present departure time, and the present arrival time of flight fe is dynamic based on every flight schedule.

Based on the actual situation, some passengers of the cancelled flight choose to endorse. Thus, the percentage αi is set to reflect the degree to which passengers will choose endorsing when a cancellation occurs for flight i. Based on the cancellation table, the number of passengers endorsing willingness can be calculated, and the number of dissatisfied passengers who want to endorse can be calculated based on the real situation.

4.4. Multiple Objective Choosing Strategy

Two types of fitness exist for every member of the population. Thus, the method for weighting two objects to optimize the population is a critical process for the LBMGA. The main idea of coping with multiobjects in this paper is to rank every solution to layer the population; within each rank, solutions need to be ordered.

The population p, which includes the present population and the filial population that is generated during the crossover operation, is assumed. In this paper, the degree of adaptation is used to distinguish rank, and the relative density is used to order the solutions in every rank. The degree of adaptation is described by the domination: if solution p dominates solution q, both objects of p are better than both objects of solution q. The solutions within the same rank cannot dominate each other.

Assume , the relative density of is , and the definition of is given in formula (13):

In formula (13), and are both sufficiently large numbers, and , where represents the value of object of and represents the number of solutions in .

The multiple target choosing operation is divided into two parts. The first part is to layer the population into different ranks to confirm the dominating relationships among solutions. The larger is the number of solutions dominated by the rank, and the smaller is the ordinal of this rank. The second part is to calculate the relative density of each solution in every rank. The solution is filled with the rank ordinal ascending and the relative density ascending in every rank until the next population is filled.

In Figure 3, a small case is shown to reveal the multiple target choosing operation. The present generation can be divided into four parts: R1, R21, R22, and others. R1 dominates R21 & R22, and R21 & R22 dominates others. Regarding the constraint of volume, only some of the solutions can be inherited by the next generation. For R21 & R22, although they are in the same rank, R22 is better than the other part of R21 based on the relative density; thus, R22 will be given priority during choosing. According to the volume, R1 and R22 are determined in order to form the next generation.

4.5. Total Algorithm Flow

A flow chart of the entire process of the loop-based GA is shown in Figure 4. At the beginning point, the original population is formed based on the natural delayed flight schedule. Then, a crossover operation is performed to form the filial population. The object values of these two populations are calculated and ranked to form the next population.

After the next population is formed, mutation and cancellation operations are performed. Then, the population undergoes a stop inspection. If the number of generations has reached the limit, circling will stop, and the final population will be calculated. If not, circling will continue.

In total, the special structural design of the algorithm ensures the improvement of the efficiency, and the settings of two types of objects make airlines try to find a balance between the economic profit and the passengers’ willingness. In other words, in this process, the passengers’ willingness is quantized and deeply considered such that the operator of airlines can clearly receive the changes of passengers’ willingness among different recovery plans and choose a preferred one among them. Thus, compared with the traditional recovery process, the passengers have more chances to receive the information of the flight adjustment early and accept it when then new recovery process mentioned in this paper is performed. The reasons can be separated into two aspects: first, the earlier passengers are informed by the airline about the recovery process means more time for the passengers to adjust their journey, and second, the preconsideration of the passenger’s willingness means more passengers can be satisfied. Hence, the overall service level of the airline can be significantly improved, as can the market competitiveness.

5. Computational Experiments

To explain and evaluate the effectiveness of the proposed recovery solution method, two types of examples are generated. First, a small example with 5 aircraft is shown to explain how our method works. Then, real-world empirical data regarding the Boeing 737 fleet of a major Chinese airline are used to evaluate the effectiveness of the recovery solution method. In this section, we will discuss the construction of problem instances and the experimental results. All computations were performed on a HASEE HYTIK5AIDIC0SE0 with i7-4270HQ CPU and 7.89 gigabits of RAM. The proposed algorithm was coded in C++.

5.1. Experiments with a Simple Case

A simple case in Table 3 for 5 aircraft and 24 flights is first introduced to analyze the performance of the proposed solution method compared with the benchmark solution method (traditional GA). In Figures 710, FN, DA, AA, DT, AT, AM and OA refer to the flight number, the departure airport, the arrival airport, the scheduled departure time, the scheduled arrival time, the aircraft fleet type, and the original aircraft ID that performs the flight, respectively.

The disruption information is set as follows: Flight 4 and Flight 5 cannot be available until 15:00 due to emergency repairs. The remaining aircraft have been available since 09:00. During the recovery process, the minimum turnaround time is assumed to be 40 min, the maximum delay time is 5 h, and the curfew time is 24:00 on the current day. For reflecting the comparison results fairly, the proposed recovery solution method and the traditional GA share the same crossover rate, mutation rate, cancellation rate, and generation limitation (40). The values of the coefficients refer to the present operating plans of the airlines and the relevant laws and regulations in China.

Based on the solving process mentioned before, LBMGA and GA both run several times and the optimal example of both algorithms is given to be analyzed in depth. In the example of GA, there is only one solution ranking first in the final generation. In the example of LBMGA, two solutions rank first together in the final generation. In order to make an intuitive comparison, the circumstance of natural delay is also given. Figures 58 show the different performances of the proposed method and the benchmark method and their contradistinctions are displayed in Table 3. Figure 5 shows the naturally delayed arrangement at the beginning of the recovery process, which indicates that the chromosome is formed based on Table 3. The time route auto-generation strategy is employed to determine which flight is cancelled, and it obeys the rule that no aircraft is rearranged.

The comparison of Table 3 and Figure 5 concludes that six flights are delayed, and six flights are cancelled. Flights CZ1849 and CZ1850 are cancelled because the AT of CZ1850 will break the curfew time. CZ1849 is in the same flight loop with CZ1850 and thus is also cancelled. The flight loops CZ1457-CZ1458 and CZ1143-CZ1144 are cancelled due to the curfew constraint. Consequently, 72 passengers who want to endorse cannot be satisfied, which results in utility losses of the disrupted passengers. Due to the high coefficient of utility loss of flight loops CZ1849-CZ1850 and CZ1143-CZ1144, which are both the end flight loops of the routes, few substance flights are available, and the total utility losses of this arrangement are relatively high.

The remaining flights of aircraft 4 are delayed by 3.5 hours, whereas the remaining flights of aircraft 5 are delayed by 4 hours. Six flights are cancelled. As a result, the economic loss of this arrangement is very high.

Figures 6 and 7 list the Rank-1 Pareto solutions of the final generation by the proposed recovery solution method. The total solution time is 1.1379 s.

The arrangement shown in Figure 6 enables CZ1109-CZ1110 to be cancelled and ten flights to be delayed, which only causes 20 endorsing passengers to be dissatisfied. The arrangement shown in Figure 6 enables CZ1109-CZ1110 and CZ1149-1150 to be cancelled and eight flights to be delayed. Compared with the arrangement in Figure 5, the total delayed time in this timetable is significantly smaller, and the cancellation fee of CZ1149-CZ1150 is relatively lower, which produces a relatively low economic loss. However, two more flights are cancelled so that 48 passengers endorsing willingness cannot be satisfied. Due to the increase in cancelled flights, the utility cost of the arrangement in Figure 7 is higher than the utility cost of the arrangement in Figure 6.

Figure 8 shows the Rank-1 Pareto solutions of the final generation by the traditional GA. First, a coding strategy is directly applied. The length of the chromosome is 24, and for each point, the number changes from zero to five. Second, a flight sequence generation strategy is employed to obtain a suitable solution. Last, all chromosome operations are randomly performed, and a punishment strategy is employed to constrain the solution.

The total solution time is 4.5532 s. The arrangement in the table enables CZ1325-CZ1326 and CZ1149-CZ1150 to be cancelled and eight flights to be delayed. Consequently, 40 passengers who want to endorse cannot be satisfied. Although the total delayed time is smaller than that in Figure 6, the economic loss is larger because the case is not such as that in Figure 6; the cancellation fee fills the gap and overflow. Compared with the arrangement in Figure 6, even though the number of cancelled flights is the same, the coefficient of utility loss of CZ1325-CZ1326 is smaller than that of CZ1109-CZ1110, which produces a smaller utility loss in Figure 8.

Table 4 provides a comparison between traditional GA and the proposed method LBMGA. Although the optimization efficiency of the proposed method is slightly better than that of traditional GA, the speed of the proposed method is substantially faster. More importantly, via the proposed method, a relatively even result can be generated, which means the willingness from the passengers can be fully considered so that the rearrangement is more likely able to satisfy the passengers

Additional analysis of the solution process indicates that approximately 65% of the time is allocated to forming a timetable in the case of the proposed method. Conversely, approximately 80% of the time is allocated to generating suitable flight sequences and calculating the punishment value in the case of the traditional GA.

According to the previous analysis, we know that the characteristics of the flight, especially the relationship between the arrival time and the curfew time, and the coefficient of utility loss have a great impact on the total benefit, on either the economy side or the utility side. In detail, the later the initial arrival time of the canceled flight is, the less the possibility that the passenger can move to another flight. Moreover, it is obvious that a high coefficient of utility loss of a flight results in high total utility loss if this flight is canceled. Therefore, the airline should pay more attention to these two factors when rearranging flights, which means that it should try to not cancel the latest flight of a day or high utility-costing flights.

Moving to the perspective of coding, this simple case shows the importance of data structure. Based on a well-designed data structure fitting the structural features of the mathematical model, plenty of time and space can be saved when solving these problems, and the additional analysis indicates that the feasibility and the punishment value process are two of the key points to improve the efficiency.

5.2. Experiments Based on Real-World Data

In this study, a large case based on a present flight schedule in China is employed to verify the effectiveness and efficiency of the proposed algorithm. Table 5 lists the scheduled information about the real data case.

The delayed aircraft are randomly chosen. For each specific delayed table, the program runs 50 times to retrieve an average. The optimization efficiency is defined in formula (14):

In formula (14), original object1 means the value of object 1 in the circumstance natural delay recovery process happens, so as the original2. The final object1 means the value of object1 under the recovery process given by the LBMGA, so as the original2.

The functionality parameters are chosen in prior experiments for excellent performance. The operating environment is an Intel Core i7-4720HQ 2.60 GHz CPU, 8 G memory, Win 10 OS and Dev-C++ 5.11 compiler. This program uses only a single core of the CPU.

In a real situation, some flights belong to flight loops. To accelerate the recovery process, only flight loops become involved with the aircraft switch, which indicates that these non-flight-loop flights will be naturally delayed if accidents occur. For the case in this section, the chromosome will be divided into three sections, which correspond to three air hubs. Only 144 flights (72 flight loops) will participate in the aircraft switch. The remaining flights will remain in their original aircraft and naturally perform a delay.

In the practical application of GA, dynamic mutation is commonly introduced to expand the search areas and avoid the solution premature. In LBMGA, the mutation rate is the same for every section of a solution, but the rate among solutions may differ.

Mutation rate for solution i:

In formula (15), means the optimization efficiency of solution . For every generation, the denominator is solid but for every solution, the molecule is different. Moreover, the molecule, , means that if the solution is more different than the best one, this solution may have more chance to mutate so that the search areas expend. More importantly, this action effectively avoids the solutions concentrating on the average so that the premature less likely appearing. Meanwhile, the dynamic mutation rate only affects the solutions above the average level and the wave of mutation is controlled to a reasonable extent so that the searching process can be more effective.

In Figure 9, the number of delayed aircraft is 15; the data are obtained from a specific running result. The outside line represents the Rank-1 Pareto solutions of generation 5, and the inside line represents the situation of final generation (100). To vividly display the data, the population size in Figure 9 is 32. From the two groups of solutions, it can be concluded that the quality of the solution has been significantly improved. Moreover, the distribution of the Rank-1 Pareto solutions has also been noticeably changed from sparse to dense. This result shows that the proposed algorithm has good convergence to determine an area of good solutions

Figure 10 shows the relationship between operation time and the number of delayed aircraft. An increase in the number of delayed aircraft causes a steady increase in operation times. As shown in Figure 10, as the number of delayed aircraft increases, the optimization efficiency gradually fluctuates and shows a rising trend. This result was observed because for every additional delayed aircraft, the disturbed flights are limited, and the number of the total disturbed flights is increased; thus, extra computation is needed to handle the new flights involved in the recovery process, which results in an increase in the time to obtain a solution. Moreover, from the analysis of the proposed algorithm, the generation of the new route timetable for every aircraft occupies most of the operations. Therefore, each newly added delayed aircraft means a new route timetable must be calculated, which increases the computational costs and time needed to obtain a solution significantly. Moreover, in the situation in which the economic loss is the only object, the relationship between the two factors above shows a trend similar to that in Figure 10, which means that the consideration of utility loss does not have significant impact on the operation speed when the scale of delayed aircraft becomes larger. Hence, considering the utility loss is reasonable in terms of in-time reaction for airlines.

It can be concluded from Figures 10 and 11 that the stability of optimization efficiency decreases as the number of delayed aircraft increases, and the complexity of the situation increases, which significantly guarantees the degree of effectiveness of the algorithm. This phenomenon can be explained by two aspects. First, as more aircraft are delayed, the feasible solutions become fewer; thus, the possible adjustment operations also become fewer and the optimization efficiency decreases. Second, a gene algorithm is based on stochastic processes; thus, the optimization efficiency shows some volatility.

In Figure 9, the data describe how the optimization efficiency changes for different minimum turnaround times and delayed aircraft; the delayed aircraft are randomly chosen for each group. When the number of delayed aircraft is 15, 20, and 30, the optimization efficiency shows a tendency to fluctuate, and the fluctuation is mainly caused by the nature of the GA, which is related to the randomness of the solution. Further analysis in the process data unveils the reason why the optimization efficiency does not change significantly for different minimum turnaround times; when the time increases, the total delayed time of naturally delayed flights in the initial solution also increases. Therefore, although the absolute quantity of the final good solutions increases with increasing maximum delay time, the optimization fluctuates. However, there is a dramatic fluctuation in the optimization when the number of delayed aircraft is 25, especially when the minimum turnaround time is 60. Analysis of the delayed aircraft shows that there are 3 aircraft only being delayed in this circumstance but not in others, and these aircraft are in the same fleet. Thus, the main extra cost is from the flights performed by this fleet. However, the number of aircraft in this fleet is very small, and only a few flights can be performed by this fleet. When the minimum turnaround time is less than 60, these flights are not severely affected. However, when the turnaround time reaches 60, they are severely disturbed. Therefore, the importance of each aircraft in the flight recovery system is not equal; the available time of some aircraft can have a large influence on the recovery process. According to the lines in Figure 9, it is obvious that the smaller turnaround time always brings the higher optimization efficiency. However, cutting down the turnaround time is a tough challenge for airlines because it highly depends on the ability level of the ground staff and the situation of the airport, which means that it is quite a challenge for the airline to behave much better than its rivals. Therefore, trying to find a different idea for improvement, such as considering the utility cost mentioned in this paper, is more feasible for the airline to improve the competitiveness. An interesting point is that if the airline only chases the minimize turnaround time, it may have a negative effect on the optimization efficiency. The line 25 mins in Figure 12 is a good example. Wiser advice is that the turnaround time should be well-designed with the time sensitivity of the total fight arrangement. This paper suggests that the airline can set their turnaround time based on the average gap-time between the flights.

In Figure 13, all of the curves reach their peaks when the maximum delay time is 5 and then show overall declines. The rising section can be interpreted as that with a relatively small maximum delay time; many flights are forced to be canceled due to the maximum delay time constraint. Thus, the higher cancellation costs and the dissatisfaction of the passengers on these canceled flights are inevitable. However, with more relaxed maximum delay time restrictions, the recovery arrangement has more flexibility to determine whether a flight is canceled or not. However, the process data show that in this circumstance, the total delayed time is likely to increase and more flights may surpass the curfew time. More importantly, the optimization process consists of the optimization of two objects. When the maximum delay time increases, object 2 shows a more pronounced downward because the coefficient in object 2 is more sensitive than object 1 when cancellation happens; moreover, as the arrival time of the flight approaches the curfew time, the utility loss will also increase. Data analysis also supports the sharp decline in the optimization of object 2. Therefore, the maximum delay time should be set wisely. Moreover, the setting process should depend on the flights in the system, especially the later flights.

Figure 14 compares the results of the traditional GA and the proposed algorithm for different maximum generations. The operation time of the proposed algorithm is significantly smaller than that of the traditional GA for all maximum generation levels. Theoretically, a larger search area increases the likelihood of finding a good solution, which is obviously the advantage of the normal gene algorithm because its searching process includes many infeasible solutions. However, the curves in Figure 13 show that when the maximum generation is relatively low, the algorithm in this paper clearly performs better, although the two algorithms convergence to a roughly equal level. Further research into the solution of the iterative process reveals that the distribution of feasible solutions is sparse but locally concentrated due to the massive constraints and the 0-1 encoding strategy. In this condition, the number of infeasible solutions is much larger than that of feasible solutions, and the searching process jumps from one area to another. Therefore, although a given sufficient maximum generation level, the normal gene algorithm can convergence to a satisfactory result, an algorithm that is based on the feature of the solutions and running within the feasible solutions can improve the search efficiency and the convergence speed significantly. Moreover, the global optimality can be ensured at the same time.

6. Conclusions

This study focused on the problem of integrated recovery of aircraft and passengers with passengers’ choice willingness between itinerary changes and the ticket refunding to put the concept of passenger-oriented service into practice and make the airlines better serve the construction of the civil aviation industry in China. The following objectives were proposed: minimizing the total economic cost of recovering flights, aircraft, and passengers and minimizing the total utility loss of passengers. Based on these objectives, an integer programming model formulation was built for describing and analyzing the problem.

The LBMGA was introduced to lead the solution process and ensure that the willingness of passengers was on the same level as the airlines’ profits. Then, a special time route auto-generation operation was designed to establish a new flight schedule to fit various constraints. A rank-based multiobject strategy was used to balance the two objectives in order to find better solutions.

The results of a small-scale example used to present the changes in a flight schedule imply that the LBMGA is more available and efficient than the normal GA. Then, a large-scale example based on a real-world situation was performed to prove the robustness and efficiency of LBMGA.

Considering the importance of passengers’ willingness, future research needs to focus on how to describe the utility more accurately and consider more circumstances, such as using statistical methods to determine the utility parameters of the passengers and expanding one-stroke to multiple-stroke. Additionally, future research will generalize this algorithm to solve the general flight recovery problem, for example, by expanding the LBMGA to embrace more types of flight networks. The pertinent improvement of the GA can increase the efficiency to a new level, such as the application of a new data structure and a new crossover or mutation operation.

Data Availability

The data, used for validating the efficiency and effectiveness of loop-based multiple objective genetic algorithm (LBMGA), can be divided into two parts. One part is a sample case. The basic data of the sample case is displayed in Table 2 in the article. The other part is real data case. The flight schedules are derived from Air China. The aircraft delay information is randomly generated according to analysis on real situation. As the scheduled database is not publicly available because of the commercial sensitivity, we could not provide the basic data. In terms of verifying the results of an article, we think that the sample case is enough for verifying the results of an article. We also have given detailed analysis for validating the results of the sample case and the real data case in the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was partially supported by the project of Central University Basic Research Fund (3072019CFM0901), the Natural Science Foundation of Heilongjiang Province (QC2016095), Social Science Foundation of Heilongjiang Province (18GLC208, 18JYC257), the project of Heilongjiang Province Postdoctoral Fund (LBH-Z15047), the project of Chinese Postdoctoral Science Foundation (2016M590276), Heilongjiang Post-Doctoral Research Initiation Fund (LBH-Q18047), the National Natural Science Foundation of China (71801061, 71401162, 71273072, 71841054), and the Major Project of Philosophy and Social Sciences Research, Ministry of Education (17JZD026).