CN109544998B - Flight time slot allocation multi-objective optimization method based on distribution estimation algorithm - Google Patents

Abstract

The invention discloses a flight time slot allocation multi-objective optimization method based on a distribution estimation algorithm, and provides a multi-objective flight time slot allocation heuristic optimization method based on the distribution estimation algorithm aiming at the current situations that the existing method is insufficient in research on the time slot allocation multi-objective optimization algorithm and few in research on flight relevance under resource constraint. The invention considers the time slot distribution problem in the ground waiting strategy, establishes a multi-objective optimization problem under the resource constraint including delay cost and fairness indexes, and also adds the capacity of a takeoff airport in the constraint condition to prevent the capacity flow imbalance of other busy airports. On the basis of a distribution estimation algorithm, a Markov network is introduced to model a problem decision variable, the correlation among multiple variables is described and mined, and the search optimization efficiency is increased; meanwhile, a fitness function with an exponential weight parameter is provided for a multi-objective optimization process, and the distribution of the pareto solution is improved under the conditions of low calculation complexity and flexible parameter configuration.

Description

Flight time slot allocation multi-objective optimization method based on distribution estimation algorithm

Technical Field

The invention belongs to the field of air traffic flow management, and particularly relates to a flight time slot allocation multi-objective optimization method based on a distribution estimation algorithm.

Background

In traffic management, the ground wait strategy is an actively implemented, very important method of traffic management. The essence is to convert the expensive and high risk air waiting into a safe and economical ground waiting. The traffic management department knows in advance that the capacity of an air traffic flow is about to decline due to the fact that an air route and a sector where an aircraft passes or a destination airport where the aircraft arrives is congested, and then the aircraft is required to wait on the ground of a take-off airport to adjust the flow of an air traffic network, so that capacity and demand balance are achieved, and the purposes of optimizing flight delay cost, improving utilization rates of airports and airspaces and reducing flight safety risks are achieved.

In the ground waiting strategy, one core technology is the flight time slot allocation. That is, an airport generates multiple time slot resources (which can be regarded as an approach time slot) due to limited capacity, and a plurality of approach flights exist at the same time, each time slot is required to be allocated to at most one flight, and each flight must be configured with a time slot later than the scheduled approach time of the flight. In general, delay cost is an objective that must be considered, and fairness indicators are important factors to consider. Therefore, the slot allocation problem is often a multi-objective optimization problem under (slot) resource constraints.

At present, on one hand, the research of the time slot allocation problem focuses on problem modeling and comprises the problems of determinacy, uncertainty, dynamics, multiple targets and the like; on the other hand, research is carried out on optimization methods, such as integer programming, stochastic programming, robust optimization and the like.

Aiming at the problem of time slot allocation, the conventional method is insufficient in research on a multi-objective optimization algorithm, and few in research on flight relevance under resource constraint. Aiming at the defects, the invention provides a multi-objective optimization method aiming at flight time slot allocation, which is established by modeling flight relevance by using a Markov network, searching and optimizing through a distribution estimation algorithm based on probability statistics and giving a new fitness function to multi-objective optimization.

Disclosure of Invention

The purpose of the invention is as follows: on the premise of ensuring airport capacity constraint, the time slot is fully utilized, multi-objective optimization of delay cost and fairness indexes is realized, and the requirement of computation time complexity is met. In the optimization process, a time slot allocation solution is given, meanwhile, the relevance among flights is mined, and data support is provided for auxiliary decision making.

The technical scheme is as follows: a flight time slot allocation multi-objective optimization method based on a distribution estimation algorithm comprises the following steps:

step 1, establishing a time slot allocation multi-target model based on efficiency and fairness;

step 2, establishing a corresponding Markov network based on a time slot distribution multi-target model;

step 3, establishing a probability matrix;

step 4, calculating and finding out a contemporary pareto solution by adopting a new fitness function;

step 5, updating a trust solution set for filing according to the current pareto solution;

step 6, learning a Markov network and updating a probability matrix;

step 7, sampling according to the probability matrix and the Markov network to generate a new solution;

and 8, reaching a termination condition to obtain a final result.

In the present invention, step 1 comprises:

step 1-1, establishing the following multi-objective function:

where N denotes the total number of flights, M denotes the total number of time slots, c_ijRepresenting the delay cost, x, when flight i is assigned to slot j_ijRepresenting a decision variable, F (x)_ij) Representing a fairness index when flight i is assigned to slot j; equation (1) represents a multi-objective function that minimizes the cost to delay and the fairness discrepancy indicator.

Step 1-2, setting constraint conditions:

c_ij＝ω_i×(t_j-SLDT_i) (5)

t_j≥SLDT_iwhen x is_ijWhen becoming 1 (6)

Wherein, ω is_iRepresenting the importance weight of the flight i in the aspect of delay time, wherein the value range is 0.0-1.0; q denotes the total number of airlines, N_qRepresenting the total number of flights, t, for the airline q_jRepresents the start time of slot j; x is the number of_ijFor decision variables, representing whether flight i is assigned to slot j, only the value 0 or 1, S L DT_iSchedule plan landing time, c, representing flight i_i ^*Delay, d (x), representing the optimum of flight i_ij) The difference, nd (x), between when flight i is assigned slot j and the optimum_ij) Representing d (x) after regularization_ij) Difference of (f)_q(x_ij) Representing the average normalized difference of the airline q; formula (2) is a decision variable of a problem, formulas (3) and (4) respectively represent a constraint condition of flight allocation and a constraint condition of time slot allocation, formula (5) is used for calculating flight delay cost, formula (6) is constraint of a flight schedule, formulas (7) and (8) calculate a difference value between actual time slot allocation and optimal allocation and perform regularization, and formulas (9) and (10) calculate variance of all participating time slot allocation companies as fairness indexes of airlines.

In the invention, the step 2 comprises the following steps: establishing a Markov network corresponding to a time slot distribution multi-target model, wherein in the established Markov network, each node represents a decision variable, namely time slot distribution selection of flights; node X_iSlot selection representing flight i, which takes on a range D (X)_i)＝{s₁,s₂…,s_MWhere M is the total number of available slots, s_MIndicating the mth available slot; in the established Markov network, the connection line between two nodes indicates that stronger coupling exists between two flights and flight time slot contention exists.

In the present invention, step 3 comprises: the following probability matrix is established:

p (t) represents the probability matrix of the t generation in the iterative computation process;

P_NMa value representing the Nth row and the Mth column of the matrix;

row i of the matrix represents the probability of flight i being assigned to each time slot j, totaling 100%;

column j of the matrix represents the probability for slot j that each flight i is assigned to that slot, totaling 100%;

when initializing the matrix, assignment is carried out by adopting equal probability, and aiming at the specific problem of time slot allocation, time constraint in a formula (6) exists, and a rule is formulated that for each row i, in M time slots, the time constraint is earlier than S L DT_iIs set to 0, will be later than S L DT_iIs set to 1/z, where z is later than S L DT_iThe number of time slots.

In the invention, the step 4 comprises the following steps:

in the step 4-1, the pareto solution is summarized into two parts: the method comprises a boundary part and a middle part, wherein the boundary part is a single-target optimization part, the middle part is a multi-target optimization part, and a boundary fitness function in a VEGA algorithm is adopted to calculate aiming at the boundary part:

EdgeFitness_f(X)＝obj_value_f(X) (13)

wherein X represents an arbitrary solution, edgeFitness_f(X) represents the fitness function of solution X at target f, obj _ value_f(X) represents the value of the function of solution X on target f;

for the middle part, the following intermediate fitness function is used for calculation:

wherein X represents an arbitrary solution, centralfiltness (X) represents the fitness function of solution X in the middle part, p (X) and q (X) represent the number dominated (dominate) by X and the number dominated X, respectively;

the fitness function values of the middle parts of all pareto solutions are more than or equal to 1, and meanwhile, all non-pareto solutions are less than 1 and more than 0;

respectively calculating the boundary fitness and the intermediate fitness according to formulas (13) and (14);

step 4-2, combining the boundary fitness and the intermediate fitness function by adopting the following composite fitness function with exponential weight parameters:

wherein D-Fitness (X) is the comprehensive fitness of X considering the composite middle part fitness and the boundary part fitness, m is the total number of targets, N _ POP is the total number of population, and Ranking_EdgeFitnessi(X) Ranking of solution X based on ith target, Ranking_{CentralFitness}(X) is the ranking of X on the intermediate fitness, ω_iAnd ω_m+1An exponential weight parameter (e.g. a two-objective optimization problem, m is 2, so there are 3 of the exponential weight parameters ω₁，ω₂，ω₃) And the scores are used for distinguishing different ranks, and meanwhile, the preference of one target can be realized.

In the invention, the step 5 comprises establishing a file in the database for storing the trust solution as a trust solution set, wherein the capacity of the file is the same as or less than the N _ POP of the population number in each iteration, and for example, the value of the N _ POP × 1/2 is taken.

And calculating by using a formula (15) to obtain fitness function values of all solutions of the current generation, selecting α% of solutions from high to low according to the fitness function values, storing the solutions into a trust solution set archive, and selecting 1- α% of trust solutions from the previous generation trust solution set archive according to the fitness function values from high to low, and storing the solutions into the trust solution set archive of the current generation.

In the present invention, step 6 comprises:

step 6-1, learning the Markov network:

by the formula (16), the variable X is calculated_iAnd Y_jMutual information amount MI (X) of_i；Y_j)：

Wherein, X_iAnd Y_jRepresenting two random variables, wherein in the time slot allocation problem, the variable specifically represents the time slot allocation of one flight; p (x)_iI D) and p (y)_j| D) respectively represent variables X based on a set of trusted solutions D_i＝x_iProbability sum Y of_j＝y_jProbability of p (x)_i,y_jI D) is X_i＝x_iAnd Y_j＝y_jA joint probability of (a);

if the mutual information quantity of the two variables is greater than a given threshold (for example, the value can be 0.5), judging that the two variables have strong correlation, and establishing a connecting line in the Markov network and mutually calling the connecting line as a neighbor; the threshold is set by the user, or dynamically updated according to the known mutual information amount in the calculation process, as shown in formula (17):

where MI represents two variables X_iAnd Y_jNDV represents the total amount of variables, and parameter β is used to control the structural complexity (which can be generally set toBetween 0.5 and 1.5);

step 6-2, updating the probability matrix:

variable X_iIn its neighbor N_iAs shown in equation (18):

wherein P (X)_i|N_i) Represents variable X_iIn all its neighbors N_iConditional probability matrix of { x }_i,1，…，x_i,j，…，x_i,YiRespectively represent variables X_iA value that can be assigned in the slot assignment, i.e. for each slot number to be assigned, x_i,YiRepresents variable X_iThe Y th_iA value that can be assigned, Y_iRepresents change X_iMaximum number of values, x_i,YiIs namely X_iThe maximum value which can be taken as a value is the latest time slot in the time slot allocation problem; p (x)_i,Valuej|n_i,Neighbork) Represents variable X_iValue_jIn time, the combined value scene of all the neighbors is Neighbor_kA conditional probability of; z_iRepresenting a set of neighbors N_iThe maximum number of the combined value-taking scene;

for the probabilities in equations (12) and (18), the boundary probability P of the current generation is first calculated_t ^*(X＝x)：

Where N (X ═ X) denotes the number of variables X that take on value X in the trusted solution, and N (d) denotes the total number of trusted solutions; the probabilistic model is updated by equation (20):

P_t(X＝x)＝(1-λ)×P_t-1(X＝x)+λ×P_t ^*(X＝x) (20)

wherein P is_t(X ═ X) and P_t-1(X ═ X) denotes the probability that the time variable X takes the value X in the t-th generation and the probability that the time variable X takes the value X in the t-1 th generation in the iterative calculation process, respectivelyλ is the learning rate from the current generation, if λ ═ 1, it means that the probabilistic model will learn completely from the current generation belief solution set archive;

step 6-3, after updating the probability model by equation (20), performing mutation operation of equation (21) by mutation probability (e.g. 5%):

P_t(X＝x)＝min(P_t(X＝x)+θ,1-) (21)

where θ represents a variance (e.g., takes 0.2), which is a small positive number (e.g., takes 0.01).

In the present invention, step 7 comprises:

step 7-1, randomly generating a feasible solution s ═ s(s)₁,s₂,…,s_n)，s_nIs the last variable X_nValue of (1), s in the slot allocation problem₁,s₂,…,s_nTime slot selection, s, on behalf of each flight_nAnd (3) representing the time slot selection of the nth flight, wherein the feasible solution s is generated by the following method: arrange all flights according to flight schedule, from the first flight x₁Initially, a time slot s is randomly selected from all the remaining time slots₁Is assigned to x₁Requires s₁Not earlier than x₁And time slot s₁Delete slot s without selection by other flights₁Move to the next flight x₂Randomly selecting by the same method until all flights are selected;

step 7-2, converting the solution s to(s)₁,s₂,…,s_n) The sequence is randomly disordered to obtain a new sequence s₁ ^*,s₂ ^*,…,s_n ^*Wherein s is_n ^*A slot selection indicating the last flight in the new flight sequence; it should be noted that: at this time, the correspondence between the flights and the time slots is not changed, but the sequence of arrangement is changed;

step 7-3, from s₁ ^*At first, set s₁ ^*Corresponding to the flight h, finding out all neighbors of the flight h variable according to the established Markov network, and solving s-s(s) according to the moment₁,s₂,…,s_n) Value of middle neighbor and corresponding conditional probability p (x)_h|N_h) Determining the probability of selecting each time slot by the flight h, and generating s by adopting a Roulette wheel selection algorithm (i.e. randomly selecting according to the probability of various values, wherein the probability of a certain value is higher, the probability of being selected is higher) according to the probability distribution₁ ^*A new value of (d); then move to s₂ ^*Recalculating s by the same method₂ ^*Until all variables are completed, the sampling of a solution is completed;

and 7-4, repeating the methods from the step 7-1 to the step 7-3 to complete the sampling of N _ POP solutions of a group, and finally forming a solution set of the current generation.

Has the advantages that: compared with the prior art, the method has the advantages that:

compared with heuristic algorithms such as a genetic algorithm and the like, the distribution estimation algorithm is adopted, so that the global search capability and the convergence speed are improved;

a Markov network structure is adopted, the variable relation is better described and explored, and more auxiliary decision making capabilities are provided;

and the low calculation complexity and the searching flexibility are realized by adopting an exponential-based weight parameter fitness function.

The method fully utilizes the time slot on the premise of ensuring the airport capacity constraint, realizes the multi-objective optimization of delay cost and fairness indexes, and meets the requirement of computation time complexity. In the optimization process, a time slot allocation solution is given, meanwhile, the relevance among flights is mined, and data support is provided for auxiliary decision making.

Drawings

The foregoing and other advantages of the invention will become more apparent from the following detailed description of the invention when taken in conjunction with the accompanying drawings.

Fig. 1 is a flow of a multi-target distribution estimation algorithm based on a markov network.

Fig. 2 is a flow of a distribution estimation algorithm.

Figure 3 is an example of a markov network architecture.

Fig. 4 is a pareto solution comprising a boundary portion and a middle portion.

Fig. 5 is a fitness function for the middle portion.

Fig. 6a is an example of a minimization problem.

Fig. 6b is a function with an exponential weighting parameter.

FIG. 7 is an update archive candidate solution set.

Fig. 8 is experimental validation results (section).

FIG. 9 is a flow chart of a flight time slot allocation multi-objective optimization method based on a distribution estimation algorithm.

Detailed Description

The invention is further explained below with reference to the drawings and the embodiments.

As shown in fig. 1, it is a schematic flow chart of a multi-target distribution estimation algorithm based on a markov network. The algorithm is based on Population-to-Population theory, each generation generates a new candidate solution, and is connected with the trust solution selected by the previous generation, and the multi-target pareto selection is carried out through the fitness function, so as to update the trust solution. And then based on the newly obtained trust solution, estimating and sampling through a probability model in a distribution estimation algorithm to obtain a candidate solution of the next generation, and circulating until a termination condition appears. It should be noted that the markov network introduction is combined with a distribution estimation algorithm to enhance the estimation and sampling process of the probability model.

Fig. 2 is a schematic flow chart of the distribution estimation algorithm employed. The emphasis is on both sampling to generate new solutions and estimating the probability model.

As shown in fig. 9, the specific process of the present invention is as follows:

step 1: and establishing a time slot allocation multi-target model based on efficiency and fairness.

(1) Parameter(s)

Parameter name parameter description

Total number of N flights

Total number of M time slots

Total number of airlines Q

i flight i, i ═ 1, …, N

j time slot j, j equals 1, …, M

Q airlines Q, Q ═ 1, …, Q

N_qTotal number of flights for airline q

c_ijDelayed cost when flight i is assigned to slot j

x_ijWhether flight i is assigned to slot j

F(x_ij) Fairness index when flight i is assigned to slot j

t_jThe start time of time slot j

SLDT_iSchedule of flight i schedules landing time

c_i ^*Delay in optimum situation for flight i

d(x_ij) Difference from optimal when flight i is assigned time slot j

nd(x_ij) Normalized d (x)_ij) Difference value

f_q(x_ij) Average normalized difference of airline q

ω_iImportance weighting of flight i in terms of delay time

(2) Multiple objective function

(3) Constraint conditions

c_ij＝ω_i×(t_j-SLDT_i) (5)

t_j≥SLDT_iWhen x is_ijWhen becoming 1 (6)

Wherein, the formula (1) represents a multi-objective function of minimizing delay cost and fairness difference index. X in formula (2)_ijTo decide the variable, representing whether flight i is assigned to slot j, only values 0 or 1 can be taken. Equations (3) and (4) ensure that for any flight to be allocated, there is and only one slot allocated to it, and for any slot to be allocated, at most, only one flight can be allocated. Equation (5) represents calculating the delay cost based on the difference between the start time of the assigned slot and the schedule time. Equation (6) ensures that for any flight, the assigned time slot must be equal to or later than the schedule time. Equations (7) and (8) calculate the difference between the delay cost under the allocation scheme and the delay cost of the optimal case and regularize. Equations (9) and (10) calculate the variance of all participating timeslot allocation companies as an airline fairness index.

Step 2: and establishing a corresponding Markov network based on the time slot distribution problem.

The markov network is a mesh undirected graph (or called a bipartite graph) and consists of a structure G and parameters Ψ. The structure G includes nodes Node and connecting edges E, so a markov network can be represented as: MN ═ G, Ψ ═ Node, E, Ψ }. As shown in fig. 3, the markov network contains 5 variables (x)₁-x₅) The line between the variables represents that there is some strong association between the two variables, which are called neighbors. For example, variable X₁Containing 2 neighbors X₂And X₃；X₂Containing 3 neighbors, X₁，X₃And X₄。

Variable X_i(D(X_i)＝{x_i ¹,…,x_i ⁿⁱValue of (x) }_i ⁿⁱRepresents X_iThe maximum value that can be taken) is calculated from the conditional probabilities of all its neighbors as shown in equation (11).

p(x_i|x-{x_i})＝p(x_i|N_i) (11)

Wherein N is_iIs a variable X_iOf p (x)_i|x-{x_i}) represents a variable X_iConditional probability, p (x), under the joint value of all other variables_i|N_i) Represents variable X_iConditional probability under the joint value of all its neighbors.

Aiming at the problem of time slot allocation optimization, the step is responsible for establishing a corresponding Markov network, including the definition and value range of each node in the network.

In a markov network, each node represents a decision variable, specifically a slot allocation selection for an airline flight. Node X_iSlot selection representing flight i, which takes on a range D (X)_i)＝{s₁,s₂…,s_MI.e. all timeslots to be allocated. The connection line between the two nodes indicates that stronger coupling and resource (flight time slot) contention exist between the two flights, and the functions of the connection line are that on one hand, neighbor information is provided for the evolution process and probability distribution is calculated; and on the other hand, certain basis is provided for subsequent operations such as time slot exchange and the like.

And step 3: and establishing a probability matrix.

The probability matrix is based on the markov network in the above step, the probability model is shown in formula (12), p (t) represents the name of the probability matrix at the t th generation in the iterative calculation process, and is a matrix of N × M, row i represents the probability of flight i being allocated to each time slot j, and the total is 100%.

In the initialization, equal probability is generally adopted for assignment, the time constraint in the formula (6) exists for the specific problem of time slot allocation, and the rule is formulated that for each row i, in M time slots, the time constraint is earlier than S L DT_iIs set to 0, will be later than S L DT_iIs set to 1/Q, where Q is later than S L DT_iThe number of time slots.

And 4, step 4: and calculating and finding out the current pareto solution by adopting a new fitness function.

The mathematical model contains two targets which need to be optimized simultaneously, and it is very important to design a set of high-efficiency fitness function to complete the solution search. In different multi-objective optimization algorithms, the fitness function plays a considerable role in search efficiency and effect. The invention provides a fitness function with an exponential weight parameter, which can improve the search efficiency and bring higher flexibility.

As shown in fig. 4, taking the minimization of the multi-objective problem as an example, the pareto solution is summarized into two parts: a boundary part (single-objective optimization part) and a middle part (multi-objective optimization part), wherein the boundary part is calculated by adopting a boundary fitness function in a VEGA algorithm:

EdgeFitness_f(X)＝obj_value_f(X) (13)

wherein X represents an arbitrary solution, edgeFitness_f(X) represents the fitness function of solution X at target f, obj _ value_f(X) represents the value of the solution X on the target f.

For the middle part, the following intermediate fitness function is used for calculation:

wherein X represents an arbitrary solution, centralfiltness (X) represents the fitness function of solution X in the middle part, p (X) and q (X) represent the number dominated (dominate) by X and the number dominated X, respectively;

as can be seen from fig. 5, the fitness function given by equation (14) can distinguish well whether it belongs to the pareto set. The centralfiltness values of all dominant solutions (pareto solutions) are 1 or more, and the more dominant the other solutions, the higher the values. Meanwhile, all non-dominant solutions are less than 1 and greater than 0. Thus, a value of 1 can clearly distinguish whether a solution is the dominant one.

According to the formulas (13) and (14), the boundary fitness and the intermediate fitness can be calculated respectively, and the calculation complexity is very low.

In order to obtain the final solution, the two fitness functions need to be considered comprehensively. There are two schemes, one is to adopt Multi-island (Multi-island) optimization, i.e. part of population is optimized according to the boundary fitness function, and the other part of population is optimized according to the intermediate fitness. This method is efficient to perform, but poses problems in that it is difficult to finally form a unified solution and the pareto solution is not highly optimal.

The algorithm uses a function with an exponential weighting parameter to combine the boundary fitness with the intermediate fitness function, as shown in equation (15):

wherein D-fitness (X) is a comprehensive fitness obtained by solving X and considering the composite intermediate part fitness and the boundary part fitness, m is a total number of targets (a double target problem, m is 2), N _ POP is a total number of people, and Ranking is performed_EdgeFitnessi(X) is the Ranking of X based on the ith target (higher EdgeFitness, higher Ranking), and Ranking_{CentralFitness}(X) is for the ranking of X at intermediate fitness (higher Central Fitness, higher ranking), ω_iAnd ω_m+1Exponentially weighted parameters (e.g. dual target optimization)Problem, m is 2, so there are 3 indexing weight parameters: omega₁，ω₂，ω₃) And the method is used for distinguishing scores of different ranks, and meanwhile, the preference of a certain target can be realized.

The indexing weight parameters are illustrated by way of example in fig. 6a and 6 b. Setting the problem as a minimization problem, in fig. 6a, consider a black solution and a white solution, the black solution being a pareto solution, and the white solution not. The black solution is characterized by being excellent in one dimension (boundary, f)₂) White solution though at boundary f₂This dimension is slightly worse than the black solution, but is better than the black solution in other dimensions. Pareto solutions need only be very "outstanding" in one dimension, rather than "good" in every dimension. Since, in order to be able to distinguish such pareto solutions clearly and quickly, fig. 6b is a function with an exponential weighting parameter, i.e. the more one solution is ranked earlier in a dimension than the other solutions, the higher its score (fitness function value) is compared to the common linear relationship.

Meanwhile, the adoption of the exponential weighting parameters has a very important role: flexibility in searching is achieved. Similar to the concept of weights, if a certain dimension does not appear in the pareto solution in sufficient quantity, the exploration of the dimension needs to be increased. By increasing the corresponding exponential weight parameter, the dimension is emphasized from the calculation of the final fitness function value, more excellent solutions can be presented in the dimension, and the distribution of the pareto solutions is also improved.

And 5: and updating the trust solution set archive according to the current pareto solution.

As shown in FIG. 7, based on the Candidate Solutions (candidates) of the present generation and the trust Solutions (recommending Data) of the previous generation, one of the two Solutions is decomposed into pareto Solutions through the complex fitness function in the previous step, and the trust solution set Archive (Archive) is used for storing the updated trust Solutions.

Step 6: and learning the Markov network and updating the probability matrix.

The algorithm is different from the traditional distribution estimation algorithm in that a Markov network is introduced, so that when the distribution estimation is carried out, besides the probability matrix needs to be updated, the structure of the Markov network also needs to be learnt again. The following are described in two aspects:

(1) learning Markov networks

The association between two nodes in the markov network means that there is a strong association (non-causal relationship) between two variables, and the assignment of a variable is also based on the "neighbors" with which it is associated, so the learning process of the markov network is an important link for the algorithm.

One method of building a network structure is by problem experts, but for most problems it is difficult to find an expert that can give a credible opinion about the association between variables. The algorithm estimates the Markov network structure using a conditional independent Test (ConditionIndreference Test), Mutual Information (MI). Mutual information has the advantage of being simple to use and fast to compute.

By the formula (16), the variable X is calculated_iAnd Y_jMutual information amount MI (X) of_i；Y_j)：

Wherein, X_iAnd Y_jTwo random variables are represented (in the slot allocation problem, the variable specifically represents the slot allocation of a flight), p (x)_iI D) and p (y)_j| D) is based onVariable X of solution set D_i＝x_iProbability sum Y of_j＝y_jProbability of p (x)_i,y_jI D) is X_i＝x_iAnd Y_j＝y_jA joint probability of (a);

if the mutual information quantity of the two variables is larger than a given threshold (for example, the value is 0.5), judging that the two variables have strong correlation, and establishing a connecting line in the Markov network and mutually calling the connecting line as a neighbor; the threshold is set by the user, or dynamically updated according to the known mutual information amount in the calculation process, as shown in formula (17):

wherein, X_iAnd Y_jRepresenting two random variables, MI representing two variables X_iAnd Y_jThe NDV represents the total number of variables, and the parameter β is used to control the structural complexity (which can typically be set between 0.5 and 1.5). if β is set to a higher value, the threshold is higher, meaning that there are fewer links and shorter computation time in the Markov network, but some of the associated information is lost.

The structure learning tends to take a long time, and therefore learning is not required every generation for search efficiency. For example, structure learning may be set to occur once every 10 generations (the probability matrix needs to be updated every generation).

(2) Updating a probability matrix

In equation (12), the probability matrix is expressed as an independent probability for each variable. Meanwhile, due to the existence of relevance in the Markov network, a neighbor-based conditional probability matrix needs to be established and updated. Variable X_iIn its neighbor N_iAs shown in equation (18),

wherein P (X)_i|N_i) Represents variable X_iIn all its neighbors N_iConditional probability matrix of { x }_i,1，…，x_i,j，…，x_i,YiRespectively represent variables X_iA value that can be assigned in the slot assignment, i.e. for each slot number to be assigned, Y_iRepresents change X_iMaximum number of values, x_i,YiRepresents X_iThe maximum value that can be taken (in the slot allocation problem, i.e. the latest slot);

represents variable X_iValue of_jIn time, the combined value scene of all the neighbors is Neighbor_kConditional probability of, Z_iRepresenting a set of neighbors N_iThe maximum number of the combined value scene.

The probabilities in equations (12) and (18) are updated according to the computed trusted solutions or the contents of the archive of the set of trusted solutions. Since the search is iterative, the boundary probability of the current generation is first calculated:

where N (X ═ X) represents the number of variables X that take the value X in the trusted solution, N (d) represents the total number of trusted solutions, and 1/| X | is set to prevent the possibility of 0 appearing, and does not substantially affect the probability calculation.

After the probability of the current generation is calculated, the probability model needs to be updated, and in order to ensure the stability of the model, the calculation is performed through the formula (20):

P_t(X＝x)＝(1-λ)×P_t-1(X＝x)+λ×P_t ^*(X＝x) (20)

wherein P is_t(X ═ X) and P_t-1(X ═ X) respectively indicates the number of iterations in the meterAnd in the calculation process, the probability that the t-th generation time variable X takes the value of X and the probability that the t-1-th generation time variable X takes the value of X are calculated, wherein lambda is the learning rate of the current generation, and if lambda is 1, the probability model can completely learn from the credible solution obtained by the current generation. Considering that an archiving mechanism has been designed in the algorithm, the value of λ can be set to be greater than 0.5.

In order to further increase the diversity of the solution and avoid falling into local optima, the mutation operation like GA is also added to the calculation process of probability update. After the probability is updated by equation (20), the mutation operation of equation (21) is performed by the mutation probability (e.g., 5%):

P_t(X＝x)＝min(P_t(X＝x)+θ,1-) (21)

where θ represents a variance (e.g., 0.2) and is a small positive number (e.g., 0.01), in order to ensure that all probabilities are less than 100%.

And 7: and sampling according to the probability matrix and the Markov network to generate a new solution.

The probability model and the network model are obtained through the last step, and the next generation of solution needs to be generated. The markov network and the Bayesian Network (BN) are very different and have no causal relationship, so there is no clear order in the process of generating a new solution by sampling. To obtain the new solution, the algorithm uses a Markov chain-based Monte Carlo method, Gibbs sampling (Gibbs Sampler). The specific process is as follows:

(1) randomly generating a feasible solution s ═ s(s)₁,s₂,…,s_n)，s_nIs the last variable X_nValue of (time slot selection for each flight in the time slot allocation problem, s)_nSlot selection indicating the nth flight). The generation method comprises the following steps: arrange all flights according to flight schedule, from the first flight x₁Initially, a time slot s is randomly selected from all the remaining time slots₁(requirement s)₁Not earlier than x₁And time slot s₁Not selected by other flights) to x₁Deleting time slot s₁Move to the next flight x₂And randomly selecting in the same method until all flight selections are finished.

(2) Converting the solution s to(s)₁,s₂,…,s_n) The sequence is randomly disordered to obtain a new sequence s₁ ^*,s₂ ^*,…,s_n ^*Wherein s is_n ^*Indicating the slot selection for the last flight in the new flight sequence. (note that, the correspondence between flights and time slots at this time is not changed, but the order of arrangement is changed).

(3) From s₁ ^*(assuming corresponding flight h) and finding all neighbors of the flight h variable according to the established Markov network, and solving s to(s) according to the moment₁,s₂,…,s_n) Value of middle neighbor and corresponding conditional probability p (x)_h|N_h) Determining the probability of selecting each time slot by the flight h, and generating s by adopting a Roulette wheel selection algorithm (i.e. randomly selecting according to the probability of various values, wherein the probability of a certain value is higher, the probability of being selected is higher) according to the probability distribution₁ ^*A new value of (d); then move to s₂ ^*Recalculating s by the same method₂ ^*Until all variables are completed, the sampling of a solution is completed.

And (4) repeating the steps (1) to (3) to finish the sampling of N _ POP solutions of a group, and finally forming a solution set of the current generation.

And 8: and reaching a termination condition to obtain a final result.

And repeating the steps until the termination condition is met. The termination condition may be set to the number of iterations, e.g., 1000 generations. When the end condition is reached, the pareto solution obtained at this point is the final result.

In order to verify the correctness of the model and the algorithm in the method, based on the flow management environment in the north China area, 18: 00-24: 00 total 6 hours of simulated flight data is used for verification and assume 22: 00-23: between 00, capacity decline and over capacity occurred in the Beijing capital airport (ZBAA).

The result of the flight slot assignment is shown in fig. 8 (partial result). The results in the figure show the estimated departure time and the calculated departure time of the affected flights, and the suggestion of the ground delay time. From the results, it can be seen intuitively that, on one hand, part of flights of a company with more flights do not need to be delayed, which indicates that the efficiency (delay time) index in the multi-objective optimization is met, on the other hand, the average delay time is relatively similar, the flight of a certain airline company is not delayed too much, and the fairness index is ensured. Through experimental simulation verification, the method can effectively solve the time slot allocation problem and obtain a satisfactory result according to multi-objective optimization calculation of problem modeling.

The present invention provides a flight time slot allocation multi-objective optimization method based on a distribution estimation algorithm, and a plurality of methods and approaches for implementing the technical scheme, and the above description is only a preferred embodiment of the present invention, it should be noted that, for those skilled in the art, a plurality of improvements and modifications may be made without departing from the principle of the present invention, and these improvements and modifications should also be regarded as the protection scope of the present invention. All the components not specified in the present embodiment can be realized by the prior art.

Claims (7)

1. A flight time slot allocation multi-objective optimization method based on a distribution estimation algorithm is characterized by comprising the following steps:

step 1, establishing a time slot allocation multi-target model based on efficiency and fairness;

step 2, establishing a corresponding Markov network based on a time slot distribution multi-target model;

step 3, establishing a probability matrix;

step 4, calculating and finding out a contemporary pareto solution by adopting a new fitness function;

step 5, updating a trust solution set for filing according to the current pareto solution;

step 6, learning a Markov network and updating a probability matrix;

step 7, sampling according to the probability matrix and the Markov network to generate a new solution;

step 8, reaching a termination condition to obtain a final result;

the step 1 comprises the following steps:

step 1-1, establishing the following multi-objective function:

where N denotes the total number of flights, M denotes the total number of time slots, c_ijRepresenting the delay cost, x, when flight i is assigned to slot j_ijRepresenting a decision variable, F (x)_ij) Representing a fairness index when flight i is assigned to slot j; formula (1) represents a multi-objective function of minimizing delay cost and fairness difference index;

step 1-2, setting constraint conditions:

c_ij＝ω_i×(t_j-SLDT_i) (5)

t_j≥SLDT_iwhen x is_ijWhen becoming 1 (6)

Wherein, ω is_iRepresenting the importance weight of the flight i in the aspect of delay time, wherein the value range is 0.0-1.0; q denotes the total number of airlines, N_qRepresenting the total number of flights, t, for the airline q_jRepresents the start time of slot j; x is the number of_ijFor decision variables, representing whether flight i is assigned to slot j, only the value 0 or 1, S L DT_iSchedule plan landing time, c, representing flight i_i ^*Delay, d (x), representing the optimum of flight i_ij) The difference, nd (x), between when flight i is assigned slot j and the optimum_ij) Representing d (x) after regularization_ij) Difference of (f)_q(x_ij) Representing the average normalized difference of the airline q; formula (2) is a decision variable of a problem, formulas (3) and (4) respectively represent a constraint condition of flight allocation and a constraint condition of time slot allocation, formula (5) is used for calculating flight delay cost, formula (6) is constraint of a flight schedule, formulas (7) and (8) calculate a difference value between actual time slot allocation and optimal allocation and perform regularization, and formulas (9) and (10) calculate variance of all participating time slot allocation companies as fairness indexes of airlines.

2. The method of claim 1, wherein step 2 comprises: establishing a Markov network corresponding to a time slot distribution multi-target model, wherein in the established Markov network, each node represents a decision variable, namely time slot distribution selection of flights; node X_iSlot selection representing flight i, which takes on a range D (X)_i)＝{s₁,s₂…,s_MWhere M is the total number of available slots, s_MIndicating the mth available slot; in the established Markov network, the connection line between two nodes indicates that stronger coupling exists between two flights and flight time slot contention exists.

3. The method of claim 2, wherein step 3 comprises: the following probability matrix is established:

p (t) represents the probability matrix of the t generation in the iterative computation process;

P_NMa value representing the Nth row and the Mth column of the matrix;

row i of the matrix represents the probability of flight i being assigned to each time slot j, totaling 100%;

column j of the matrix represents the probability for slot j that each flight i is assigned to that slot, totaling 100%;

when initializing the matrix, assignment is carried out by adopting equal probability, and aiming at the specific problem of time slot allocation, time constraint in a formula (6) exists, and a rule is formulated that for each row i, in M time slots, the time constraint is earlier than S L DT_iIs set to 0, will be later than S L DT_iIs set to 1/z, where z is later than S L DT_iThe number of time slots.

4. The method of claim 3, wherein step 4 comprises:

in the step 4-1, the pareto solution is summarized into two parts: the method comprises a boundary part and a middle part, wherein the boundary part is a single-target optimization part, the middle part is a multi-target optimization part, and a boundary fitness function in a VEGA algorithm is adopted to calculate aiming at the boundary part:

EdgeFitness_f(X)＝obj_value_f(X) (13)

wherein X represents an arbitrary solution, edgeFitness_f(X) represents the fitness function of solution X at target f, obj _ value_f(X) represents the value of the function of solution X on target f;

for the middle part, the following intermediate fitness function is used for calculation:

wherein X represents an arbitrary solution, CentralFitness (X) represents the fitness function of solution X in the middle part, p (X) and q (X) represent the number dominated by X and the number dominated by X, respectively;

the fitness function values of the middle parts of all pareto solutions are more than or equal to 1, and meanwhile, all non-pareto solutions are less than 1 and more than 0;

respectively calculating the boundary fitness and the intermediate fitness according to formulas (13) and (14);

step 4-2, combining the boundary fitness and the intermediate fitness function by adopting the following composite fitness function with exponential weight parameters:

wherein D-Fitness (X) is the comprehensive fitness of X considering the composite middle part fitness and the boundary part fitness, m is the total number of targets, N _ POP is the total number of population, and Ranking_EdgeFitnessi(X) Ranking of solution X based on ith target, Ranking_{CentralFitness}(X) is the ranking of X on the intermediate fitness, ω_iAnd ω_m+1The weight parameters are indexed to distinguish scores of different ranks and simultaneously realize the preference of one target.

5. The method of claim 4, wherein step 5 comprises: establishing a file in a database for storing the trust solution as a trust solution set, wherein the capacity of the file is the same as or less than the N _ POP of the population number in each iteration;

and calculating by using a formula (15) to obtain fitness function values of all solutions of the current generation, selecting α% of solutions from high to low according to the fitness function values, storing the solutions into a trust solution set archive, and selecting 1- α% of trust solutions from the previous generation trust solution set archive according to the fitness function values from high to low, and storing the solutions into the trust solution set archive of the current generation.

6. The method of claim 5, wherein step 6 comprises:

step 6-1, learning the Markov network:

by the formula (16), the variable X is calculated_iAnd Y_jMutual information amount MI (X) of_i；Y_j)：

Wherein, X_iAnd Y_jRepresenting two random variables, wherein in the time slot allocation problem, the variable specifically represents the time slot allocation of one flight; p (x)_iI D) and p (y)_j| D) respectively represent variables X based on a set of trusted solutions D_i＝x_iProbability sum Y of_j＝y_jProbability of p (x)_i,y_jI D) is X_i＝x_iAnd Y_j＝y_jA joint probability of (a);

if the mutual information quantity of the two variables is greater than a given threshold, judging that the two variables have strong correlation, and establishing a connecting line in the Markov network and mutually calling the connecting line as a neighbor; the threshold is set by the user, or dynamically updated according to the known mutual information amount in the calculation process, as shown in formula (17):

where MI represents two variables X_iAnd Y_jNDV represents the total amount of variables, and parameter β is used to control the structural complexity;

step 6-2, updating the probability matrix:

variable X_iIn its neighbor N_iAs shown in equation (18):

wherein P (X)_i|N_i) Represents variable X_iIn all its neighbors N_iConditional probability matrix of { x }_i,1，…，x_i,j，…，x_i,YiRespectively represent variables X_iA value that can be assigned in the slot assignment, i.e. for each slot number to be assigned, x_i,YiRepresents variable X_iThe Y th_iA value that can be assigned, Y_iRepresents change X_iMaximum number of values, x_i,YiIs namely X_iThe maximum value which can be taken as a value is the latest time slot in the time slot allocation problem;

represents variable X_iValue_jIn time, the combined value scene of all the neighbors is Neighbor_kA conditional probability of; z_iRepresenting a set of neighbors N_iThe maximum number of the combined value-taking scene;

for the probabilities in equations (12) and (18), the boundary probability P of the current generation is first calculated_t ^*(X＝x)：

Where N (X ═ X) denotes the number of variables X that take on value X in the trusted solution, and N (d) denotes the total number of trusted solutions; the probabilistic model is updated by equation (20):

P_t(X＝x)＝(1-λ)×P_t-1(X＝x)+λ×P_t ^*(X＝x) (20)

wherein P is_t(X ═ X) and P_t-1(X ═ X) respectively represents the probability that the time variable X of the t-th generation takes the X value and the probability that the time variable X of the t-1 th generation takes the X value in the iterative computation process, wherein lambda is the learning rate from the current generation, and if lambda is 1, the probability model is completely learned from the current generation trust solution set archive;

step 6-3, after updating the probability model by formula (20), performing mutation operation of formula (21) by mutation probability:

P_t(X＝x)＝min(P_t(X＝x)+θ,1-) (21)

where θ represents the variance and is a small positive number.

7. The method of claim 6, wherein step 7 comprises:

step 7-1, randomly generating a feasible solution s ═ s(s)₁,s₂,…,s_n)，s_nIs the last variable X_nValue of (1), s in the slot allocation problem₁,s₂,…,s_nTime slot selection, s, on behalf of each flight_nAnd (3) representing the time slot selection of the nth flight, wherein the feasible solution s is generated by the following method: arrange all flights according to flight schedule, from the first flight x₁Initially, a time slot s is randomly selected from all the remaining time slots₁Is assigned to x₁Requires s₁Not earlier than x₁And time slot s₁Delete slot s without selection by other flights₁Move to the next flight x₂Randomly selecting by the same method until all flights are selected;

step 7-2, converting the solution s to(s)₁,s₂,…,s_n) The sequence is randomly disordered to obtain a new sequence s₁ ^*,s₂ ^*,…,s_n ^*Wherein s is_n ^*A slot selection indicating the last flight in the new flight sequence;

step 7-3, from s₁ ^*At first, set s₁ ^*Corresponding to the flight h, finding out all neighbors of the flight h variable according to the established Markov network, and solving s-s(s) according to the moment₁,s₂,…,s_n) Value of middle neighbor and corresponding conditional probability p (x)_h|N_h) Determining the probability of selecting each time slot for flight h, and generating s by adopting a roulette selection algorithm according to the probability distribution₁ ^*A new value of (d); then move to s₂ ^*Recalculating s by the same method₂ ^*Until all variables are completed, the sampling of a solution is completed;

and 7-4, repeating the methods from the step 7-1 to the step 7-3 to complete the sampling of N _ POP solutions of a group, and finally forming a solution set of the current generation.

Images