Abstract
Aiming at the phenomenon of a large number of flight delays in the terminal area makes a reasonable scheduling for the approach and departure flights, which will minimize flight delay losses and improve runway utilization. This paper considered factors such as operating conditions and safety interval of multi runways; the maximum throughput and minimum flight delay losses as well as robustness were taken as objective functions; the model of optimization scheduling of approach and departure flights was established. Finally, the genetic algorithm was introduced to solve the model. The results showed that, in the program whose advance is not counted as a loss, its runway throughput is improved by 18.4%, the delay losses are reduced by 85.8%, and the robustness is increased by 20% compared with the results of FCFS (first come first served) algorithm, while, compared with the program whose advance is counted as a loss, the runway throughput is improved by 15.16%, flight delay losses are decreased by 75.64%, and the robustness is also increased by 20%. The algorithm can improve the efficiency and reduce delay losses effectively and reduce the workload of controllers, thereby improving economic results.
1. Introduction
Flight delay is one of the problems to restrict the development of the world aviation industry and also the main source for the airline passengers’ dissatisfaction with the service. At present, countries around the world have taken various measures to reduce flight delay; many airports in our country improve capacity by increasing the number of runways and other methods, thereby alleviating the flight delay. Because the airport runways are the resources of the air traffic system, within a certain time, making good use of runway resources can effectively alleviate the flight delay. In addition, the approach and departure flight scheduling of terminal area also plays an important role in reducing the flight delay, so reasonable scheduling of flights has an influence on ensuring flight safety, improving resource utilization, and reducing the loss of delay as well as improving airline credibility and so forth. In our country, FCFS (first come first served) was used for aircraft terminal controllers to make a sorting for the arrival aircraft. While the fact shows that FCFS is not the best strategy to maximize the use of existing airport capacity so is reducing the average delay losses [1].
Starting in early twentieth Century, domestic and foreign scholars have launched the research on scheduling flight optimization problem. Dear and Sherif proposed a methodology for sequencing and scheduling of aircraft in high density terminal areas. Termed constrained position shifting (CPS), this methodology was examined and its effectiveness was tested. Potential capacity improvements were noted over the first-come, first-served, runway (FCFS-RW) strategy, especially during peak periods [2]. Brinton designed an algorithm called the implicit enumeration (IE) scheduling algorithm to optimize arrival aircraft sequences and schedules. The approach taken in the algorithm allows the planning process to be expanded to include runway assignment. Initial results from using the IE algorithm for runway assignment indicated that significant performance enhancements were possible when both runway and sequence assignments were considered in the scheduling process [3]. Abela et al. described two approaches for solving the problem of scheduling aircraft landing times. They also formulated this problem as a mixed integer program (MIP) and developed a branch and bound algorithm for its solution [4]. Ronbinson designed a fuzzy reasoning-based method for scheduling air traffic in the terminal area. This method base applied fuzzy reasoning to evaluate propositions that consider both performance criteria and workload criteria, such as delay reduction and conflict avoidance [5]. Andersson et al. described a novel optimization-based approach to benefits analysis that allowed us to analyze new collaborative air traffic management tools. They applied this method to analysis of the collaborative arrival planner (CAP), a concept under development by NASA as part of a suite of decision-support tools for improved air traffic management [6]. Saraf and Slater developed a new efficient scheduling algorithm, which used combinatorial optimization techniques to find the optimal arrival aircraft sequence and the optimal STAs for the aircraft at a certain reference point, given that the maximum number of position switches from—the FCFS order was specified. The algorithm also enabled prioritization of aircraft according to airline preferences and distributed the allocated delays optimally among all the sectors that lie along the route of the aircraft [7]. Lee and Balakrishnan presented a dynamic programming algorithm for determining the minimum cost arrival schedule at an airport, given aircraft-dependent delay costs. They also proposed an approach which made it possible to evaluate tradeoffs in terminal-area schedules, through the comparison of throughput maximizing and delay minimizing schedules. A comprehensive analysis of the tradeoffs between average delay and fuel costs was also conducted, accounting for both the cost of delay and the increased rate of fuel burn incurred by an accelerating aircraft [8]. Hancerliogullari et al. examined the aircraft sequencing problem (ASP) over multiple runways, under mixed mode operations with the objective of minimizing the total weighted tardiness of aircraft landings and departures simultaneously. It can be modeled as a parallel machine scheduling problem with unequal ready-times, target times, and deadlines. Due to the problem being NP-hard, greedy heuristics and metaheuristics are applied in this paper to obtain solutions in reasonable computation times [9].
Compared to other developed countries, the research of air traffic flow management started late in China. Ye and Tao established a dynamic model of the aircraft sequencing in terminal area and dynamic scheduling which was based on genetic algorithm and made a reasonable arrangement for the flight landing sequence, thereby improving the runway utilization, reducing flight delay losses [10]. Cheng et al. put forward adaptive genetic algorithm in arrival aircraft scheduling, solving the problem of landing aircraft scheduling on a single runway, which was significantly better than FCFS algorithm [11]. Wang et al. proposed a method of approach and departure flight scheduling based on time colored Petri network (TCPN). It can obviously reduce the flight delay losses and to a certain extent help the ATC personnel making scheduling decision [12]. Li et al. proposed a dynamic optimization model and an optimal flow allocation method of flight scheduling for terminal area, on the basis of traffic optimization, integrated the type of aircrafts, wake separation time, and other factors. This method can effectively optimize allocation of airport approach and departure traffic, while ensuring minimum flight separation among aircraft, to determine the optimal landing sequence, making a combination of the aircraft approach and departure traffic optimization allocation and the plan of arriving and landing [13]; Cao and Zhang built a model of airport’s servicing queue by using Chapman-Kolmogorov equation theory of Markov chain. It was proved that this model could reflect the delays of flight intuitively. And the flight ordering strategy would be formulated according to the delaying condition [14]. Zhang et al. introduced the receding horizon control (RHC) strategy and established a flight sequencing model. It was solved by using the genetic algorithm based on RHC strategy (RHC-GA). The method was proved effectively [15].
Although the scholars had done a lot of research on flights sorting optimization methods, still they had disadvantages. For example, some models were established in the ideal situations and did not take into account the actual situation of the weather, air traffic control, airport surrounding environment, and so forth. Moreover most models were for single runway, rather than the multirunways, and did not meet with the trend of current development of the airport.In addition, the problem of taking off and landing of flight scheduling on multirunways is more complex with the increasement of runways, and the influence of human factors will be deepened. So it is needed to further improve the algorithm, so as to run closer to the real-time control. Considering the weather, air traffic control, route and other factors, the model was established in this paper to make optimization of approach and departure flight scheduling. The model was solved by genetic algorithm, thus improving the runway throughput and reducing flight delay losses. At last, the economic benefit of the airlines is improved.
Robust optimization solution is obtained in every possible case (but the probability is unknown in each case). It is the solution which is with small deviation from the optimal solution. For different problems, the robust optimization methods are also different. In order to solve the robust optimization problem in this paper, make use of the advantage of genetic algorithm, combine with the robustness, and find out the flight sequence with minimal changes, then make simulation. This method is proved to be correct.
2. Model and Methods
2.1. Problem Description
In recent years, with the rapid development of China’s civil aviation, busy airport passenger throughput increased significantly, resulting in the aviation market demand for sustained growth and the increasement of flight traffic, coming with it is the problem of flight delay, and,thus, leading to the dissatisfaction of passengers, or even making a confliction. Therefore, considering a variety of factors, such as weather, traffic control, air traffic control, route, and other restrictions, make a reasonable scheduling of approach and departure flights under the premise of ensuring safety, which will minimize the flight delay and alleviate the dissatisfaction of passengers caused by flight delay.
Air traffic controllers make a reasonable allocation of the approach and departure flights waiting in the queue. Distribute taxiways and runways for the incoming flights to make them a smooth landing, reaching the apron safely; while for the departure flights, arrange appropriate taxiways and runways, making it leave the airport safely and smoothly take off. Approach and departure flights schematic is showed in Figure 1. The flight scheduling problem considered in this paper, mainly is on account of the basic principle—first come first served (FCFS), from different types of aircraft must maintain different “minimum safe separation standards,” by rearranging the order of the queue for flights, searching for all the possible flight queues, finding a finally flight queue with maximum runway throughput, minimum flight delay losses, and minimal shift.
2.2. Model Introduction
Flight sorting is a dynamic continuous process. It is needed to make corresponding adjustments according to the change of real-time information. Assume that there are a lot of flights coming into the terminal area, waiting for the traffic controllers to make arrangements for them. It is required that meeting the restrictions of minimum safety interval among aircraft, arranging the order of flights, thus minimizing the time of the completion of approach and departure flights scheduling.
The robust optimization method of approach and departure flight scheduling is studied in this paper, which assumes the following:(1)there are flights of airlines waiting for scheduling. Airline controllers make a dynamic sorting for flights and distribute them into taxiways. Under the premise of ensuring the safety, making a reasonable arrangement of the time of approach and departure flights, the runways and the order, which will minimize the total delay of the whole flight queue;(2)the estimated time and actual time of each flight are not the same and they can be determined at the time 0;(3)the airport studied in this paper contains multiple parallel runways, and each of them meets independent operation standards;(4)the information of approach and departure flight (including flight number, flight model, and the estimated arriving time of flights) is known;(5)the capacity of the airport meets the assumption. That is the number of flight and the time distributed within the range of airport capacity licensing;(6)the approach flights do not delay when they take-off at the departure airport, and arrive on time at the terminal area waiting for landing.
2.3. Parameters and Symbols
: is the collection of approach flights, ; : is the collection of departure flights, ; : is the collection of all the flights , , , ; : is the collection of independent parallel runways, ; : is the collection of taxiways, ; : The maximum delay time of approach flight during the time of arriving at the airport; : The maximum delay time of departure flight during the time of leaving the airport; : is the actual time of flight arriving at the airport; : is the estimated time of flight arriving at the airport; : is the actual time of flight leaving the airport; : is the estimated time of fight leaving the airport; : is the arrival time of flight after being optimized; : is the leave time of flight after being optimized; : is the minimum safety interval time between two successive approaching flights (they are landing on the same runway and ); different flights have different safety interval, meanwhile, on the same runway, if a same type of flight has different types of approach and departure, its safety interval is not the same; : is the unit time delay cost of approach flight ; : is the unit time delay cost of departure flight ; different types of flights have different unit time delay costs; the larger flight has higher unit time delay cost.
2.4. Modeling
2.4.1. The Objective Function of Runway Throughput and Its Fairness
Suppose that the actual scheduling time of the first flight on the runway is and the actual scheduling time of the last flight is , so all the time for scheduling the flight queue is . In order to maximize the throughput of the runway, it is required that the time finishing the scheduling of the flight queue is minimal. So the smaller is, the larger the runway throughput is. At the same time, the number of flights on each runway should be close to each other, thus ensuring the fairness among multirunways.
The objective function of the maximum throughput of the runway is
2.4.2. The Objective Function of Flight Delay Losses
Taking into account the actual situation, in most cases, either the approach or departure flights, the actual arrival time is almost impossible to be consistent with the estimated time, so we make the following provision. The following two situations are not treated as delay: the actual time is earlier or later than the estimated time within ; all the flights which are in advance of the estimated time landing or taking off. The model is established as follows.
The delay time of the flight is The delay losses of the flight queue.
Suppose that there are flights whose actual time of arriving or leaving is later than their estimated time, so their total delay losses is ; meanwhile, there are flights whose actual time of arriving or leaving is earlier than their estimated time, because their delay time is 0, so their total delay losses are ;
That is, The total delay losses of the entire flight queue are Assume that there are kinds of sorting programs, and the delay losses of each program is , so the objective function of flight delay losses are
2.4.3. The Objective Function of Robustness
The objective of robustness can be achieved by adjusting the minimal flight sorties, that is, it is needed to try to make the flight schedule time consistent with the estimated time, So that, the smaller adjustments, the better robustness. Meanwhile, reduce the workload of controllers.
Therefore, the objective function robustness can be measured by the number of flights which the schedule time is not the same as the estimated time; the smaller the robustness the better.
We introduce variables 0, 1: The objective function of robustness is It also can be expressed as a percentage:
2.5. Constraints
Equation (9) indicates that adjacent flights should meet certain wake separation and vortex separation. is the number of flight, is the scheduling time allocated to the flight , is the vortex separation, and is the wake separation.
Equation (10) is to restrain the delay time of the approach flight which should not be larger than the maximum approach delay time.
Equation (11) is to restrain the delay time of the departure flight which should not be larger than the maximum departure delay time.
Equations (12) and (13) are the constraint of runways; each flight only has one runway for arriving or leaving.
Equation (14) indicates that if there are two continuous flights arriving or leaving on the same runway (suppose that the flight is the former and the flight is the latter) the time interval they need to meet.
Equation (15) indicates that for two approach flights, and , either is the leading, is the tailing, or is the tailing, is the leading.
Equation (16) indicates that any approach flight only has one leading.
Equation (17) indicates that any approach flight only has one tailing.
For the mixed types of aircraft, the International Civil Aviation Organization (ICAO) specifies the minimal interval standards, denoted as ; the distance interval and its corresponding time interval is shown in Table 1. (The table represents the condition of windless).
2.6. Genetic Algorithm
The basic idea of the genetic algorithm is to simulate the natural process of genetic mechanism and biological evolution that form a process to search for the optimal solution. The characteristic of it is that the treatment objects of it are the parameters of the code collection rather than the problem parameters themselves. Besides, its searching process is not influenced by the constraint of connection of the optimization functions, also the optimization functions do not need to be differentiable. Therefore it has better ability of searching [16].
2.6.1. Chromosome Coding Scheme
According to the real characteristic of the flight scheduling problem, in this paper, real-coded schema is adopted, that is, the digital serial number encoding. Every flight queue is a chromosome, and each flight in the queue is a gene value. That is, every chromosome represents a scheduling scheme. For example, 6 3 4 1 2 5 7 is a chromosome, in it, 6 represents the first flight to be dispatched, while 3 4 1 2 5 7 is sequentially scheduled. The genotype of the coding scheme is a real number and the phenotype of it is the flight queue. The advantage of it is that making the genotype and phenotype correspondence.
2.6.2. Fitness Function
The fitness function is also known as the evaluation function. It is a symbol to judge the individual is good or bad based on the objective functions. It is also the driving force of the evolution process. Because the fitness function is always nonnegative, so under any circumstances hoping its value the bigger the better.
Real-coded schema is used in this paper, so by using the sorting method according to performance degree to determine its fitness [17]. That is making descending order of the individual according to their quality of the performance. As for each objective , all the individuals will generate a feasible sequence basis for the value of the objective function. Sorting for each objective, we can get overall performance of the individual for all objective functions. The calculation of fitness is as follows: Equation (19), ; ; is the number of objective functions; is all individuals of the population; is the individual of the population; is the serial number of the individual after all individual sorting in view of the objective ; is the fitness value for the individual on account of the objective ; is the comprehensive fitness value of the individual on account of all objectives; and is the constant between , for increasing the fitness value when the performance of individual is optimal. As can be seen from the equations above, for optimum overall performance of individuals can get better fitness to gain more opportunities to participate in evolution.
2.6.3. Genetic Operators
Make genetic operations of selection, crossover, and mutation of individual to produce more new individuals in the genetic algorithm. Although with the evolution of populations, it will produce more and more excellent individuals, because selection, crossover, mutation, and other operations are random, so it may also destroy the current individual with the best fitness. In order to choose the best to retain to the next generation, using the optimal preservation strategy to make elimination or survival; namely, the individual with best fitness does not involve in the operation of crossover and mutation, instead of replacing the individual with the worst fitness produced by crossover and mutation in the population.
(1) Selection Operator. The selection operation is that selecting the excellent individuals with a certain probability from the data of the upper level, then reproduces next generation. The purpose of it is to copy the good genes of the individual who has better fitness to the next generation.
The roulette selection operator is used in this paper; namely, the probability of the fitness in proportion decides the possibility of its descendants going or staying. If a certain individual is , the fitness is , the probability to be selected can be expressed as
(2) Crossover Operator. Crossover means selecting two individuals randomly from the population, by making exchange and combination of the chromosomes to inherit outstanding feature from the parents, thus producing new excellent individuals. Single-point crossover mapping method is used in this paper [18]; this method is an improvement of the partially matched crossover (PMX) proposed by Goldberg and Lingle. Now make an example: write two flight queues and as parent queues and randomly assign multirunways to the flight sequence. Here, randomly select crossover point, shown as follows: After performing a simple single-point crossover, we can get: Check traverse repetition before the intersection point of the code. Make exchange one by one according to the position mapping relation.
For For So
In the process of crossover, the runway of the flight does not make chiasma but allocates runways randomly to the new fight queue after chiasma. That is allocating multirunways for as well as . Therefore, it will expand the range of feasible solutions.
(3) Mutation Operator. Mutation operator is to improve the local search ability of genetic algorithm and avoid falling into local optimal solutions at the same time. It is also an important means of maintaining the diversity of population. Uniform mutation is used in this paper; it is that, random functions which are uniform distribution within a certain range replace original gene on each gene locus with a certain small possibility. The method is particularly applicable to the primary operational phase of the genetic algorithm and makes the search points able to move throughout the search space freely, thereby increasing the diversity of the population, so the algorithm can deal with more patterns.
Suppose that is a gene, is the change point, its value ranges from to . After making uniform mutation of the individual at this point, we can get a new individual , the new gene value of the change point is
In Equation (26), is a random number which is in accordance with uniform distribution within a range.
2.6.4. Algorithm Steps
The process of the algorithm is showed in Figure 2; the specific steps are as follows.
Step 1. Determine the genetic strategy, including population size , the methods of selection, crossover, and mutation. At the same time, determine the probability of crossover and the probability of mutation , as well as other genetic parameters;
Step 2. Define fitness function , as shown in (18);
Step 3. Generate initialization population randomly;
Step 4. Calculate the objective functions of the chromosome which is corresponding with flight queue;
Step 5. Calculate the fitness value of the individuals in the population, as shown in (19);
Step 6. Find the best individual in the population under the current conditions;
Step 7. Making a judgment of whether the evolution algebra meets the condition of smaller than the maximal algebra. If it is, algebra plus 1, then keep it in order. If it is not, turn to the Step 11;
Step 8 (roulette selection). Make crossover operation of the chromosome by Single-point crossover mapping method; make mutation operation of the chromosome by uniform mutation method, as shown in (20), (26);
Step 9. Execute Step 4~Step 6;
Step 10. Assess the effects of the genetic algorithm;
Step 11. Output the optimal function value then get the optimized flight sequence.
2.6.5. Processing Constraints
The constraints involved in the model are more and complicated, so they need to be processed. For the constraints such as , treat them as heritable variations. As for other types of constraints, such as , deal with them with penalty function. That is adding a penalty function to the fitness value on the individuals which violate constraints. Consider In Equation (27), is the original fitness value; is the new fitness function after adjustment; is a penalty function factor which is more than 0; and is the penalty function; it is a positive value when is does not meet the constraints.
As for , the penalty function can be written as . If , it does not need to be punished; If , the penalty is . Here, penalty factor needs to be set according to specific conditions.
3. Simulation and Analysis
3.1. Simulation
According to the standard of the International Civil Aviation Organization (ICAO), the aircraft is divided into 3 categories in accordance with the wake intensity. The study showed that the unit time delay losses of the approach flight are larger than those of the departure flight. The operating costs are showed in Table 2 [19, 20]. Unit time delay costs of approach and departure flights are different from different airlines, but the difference was not significant.
30 flight data of a continuous period base on Chengdu Shuangliu International Airport. There are 15 approach flights and 15 departure flights, 2 mutually independent parallel runway. s. The initial flight data is showed in Tables 3 and 4. Table 3 shows the departure flights and Table 4 shows the approach flights, including airlines, flight number, type, unit time delay loss, estimated time and actual time of departure or approach flights, runway and delay losses.
Running the above multiobjective genetic algorithm, the crossover probability is 0.9, the mutation probability is 0.1, the generation gap is 0.9, the elimination rate is 0.2, penalty factor is 0.2, the size of population is 100, and the evolution algebra is 2500. In the fitness function, , the selection criteria of the three objective functions is compared with the FCFS, that is lower than FCFS. The range of change point is , .
Select 7 schemes after optimization, the sequence results of the flights are as shown in Table 5.
Calculate objective function value of each program, and the results are compared with the FCFS algorithm, as shown in Table 6.
According to Table 6, Figure 3 is painted. It shows the comparison of FCFS and 7 kinds of optimization schemes, including three objectives: delay losses, runway capacity, and robustness (which is represented by the shift of flights).
3.2. Analysis
According to the data in the Table 4, the average delay losses of the 7 kinds of schemes which are optimized are 45810.69. Compared with FCFS, it has declined by 78.97%, the runway capacity improved by 18.00%, and the robustness also improved by 13.14%. It showed that the delay losses, runway capacity, and robustness all had been optimized compared with the FCFS algorithm.
There are three objectives in this paper, namely, minimum delay losses, maximal runway capacity as well as best robustness. Using multiobjectives genetic algorithm to solve the problem, then one can get a group of Pareto- optimal solutions. In this way, the decision makers can choose one of the schemes as the final scheduling program according to their own preferences.
In order to get a better program, design an evaluation function to make an assessment of the optimization schemes. Suppose that are the changing percentage of the runway capacity, delay losses, and robustness, respectively. Meantime, is the weight of each objective function. After calculation and comparison, find the percentage of maximum scheme as the final program. The evaluation function is as follows: Set constraints of the evaluation function to select the best program. The constraints are as follows: Make a comprehensive evaluation of the three objective functions, . Now define the optimization program 1; it is the program that does not regard the in advance as a loss and makes a comparison with FCFS. The delay losses could be calculated by (4). At the same time, in some cases, treating the flights which are arriving or leaving in advance as abnormal flights that need to be punished also results in the loss of the flight. But the unit time loss is half of delay loss. So define the optimization program 2; it is the program that regarding the in advance as a loss. Thus, the flight sorting results of these three schemes are showed in Table 7.
Moreover, making a comparison of the three target values of these three schemes, get the data showed in Table 8. And Figure 4 is painted according to Table 8, showing the comparative results.
As can be seen from the Table 8 and Figure 4:(1) Optimization Program 1. The delay losses is declined by 85.8% , the delay losses of 30 flights of each program is showed in Figure 5. In addition to some individual flights, the majority of flight delay losses are significantly reduced. Runway capacity is improved by 18.4% , the time of entire flight queue scheduling has saved 290 s, increasing the runway capacity. Besides, there is a more even distribution of these two runways, thus increasing the utilization of system resources. The robustness is improved by 24% , that is, the workload of the air traffic controllers declined by 24%, the number of flight shift is declined from 25 of FCFS to 20. In the flight queue after optimization, there are 4 flights arriving in advance, thus saving the fuel consumption for air waiting and reducing operating costs. It is a great advantage compared with the FCFS algorithm.(2) Optimization Program 2. The delay losses are declined by 75.64%, runway capacity is improved by 15.16%, and the robustness is improved by 20%.In addition, compare it with the genetic algorithm used in other papers. For example, compared with the results of literature 10, the delay losses are declined by 31.42%. Similarly, compared with the results of literature 11, the delay losses are also declined by 20.35%. It is obvious that the results of this paper are largely reduced flight delay losses, the results are optimal. In this paper, consider three objectives: delay losses, runway capacity, and robustness comprehensively, use multiobjective genetic algorithm to obtain multiple Pareto optimal solutions, and provide a variety of options for decision-makers, so decision makers can make the decision according to their preferences.
In summary, the results got by using multiobjective genetic algorithm in this paper have improved a lot compared with the FCFS algorithm. Decision makers can select the appropriate scheduling solution according to their needs based on the method of this paper, in order to obtain the satisfactory results.
4. Conclusions
In this paper the optimization problem of multirunway approach and departure flight sorting based on genetic algorithm is discussed, the targets of runway capacity, delay losses, and robustness are made; then it presents a multiobjective simulation model to solve it by using genetic algorithms. The simulation results show that the model and algorithm established in this paper on the flight scheduling problem have been greatly improved, not only runway capacity is improved, the delay losses is reduced, but it also reduces the workload of controllers and enhances the robustness of flight, which has high optimization efficiency. The method is effective and feasible in solving scheduling problems of flight at terminal area in reality; it can meet the requirements of operation controllers in real-time. The multiobjective algorithm is used in this paper, thus the decision makers can choose the final solution in the Pareto solutions. In addition, the optimization algorithm does not consider the impact on the airlines caused by flight scheduling, namely, lacking of the research on the fairness of the airlines. Therefore, the next job is to study the fairness of the airlines, consider the influence of airlines on flight sorting queue, make the improvement to the algorithm, and put forward a more perfect optimization program.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors thank the reviewers for helping them to improve this paper. This work is supported by National Natural Science Foundation of China, no. 61262002 and the Fundamental Research Funds for the Central Universities, no. NS2014064.