CN109544998B - A Multi-objective Optimization Method for Flight Slot Allocation 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

A Multi-objective Optimization Method for Flight Slot Allocation Based on Distribution Estimation Algorithm

technical field

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

Background technique

In traffic management, the ground waiting strategy is an active and very important traffic management method. The essence is to convert expensive and high-risk air waiting into safe and economical ground waiting. The traffic management department knows in advance that the air route, sector or destination airport that the aircraft passes through is about to reduce the capacity of the air traffic flow, so the aircraft is required to wait on the ground at the departure airport to adjust the flow of the air traffic network to achieve a balance between capacity and demand. , so as to optimize flight delay costs, improve airport and airspace utilization and reduce flight safety risks.

In the ground holding strategy, a core technology is flight slot allocation. That is, the airport generates multiple time slot resources (which can be regarded as arrival time slots) due to limited capacity, and there are multiple arrival flights at the same time. It is required that each time slot can only be assigned to at most one flight, and each flight must be configured with A time slot later than the scheduled arrival time of the flight. Generally speaking, the cost of delay is the target that must be considered. In addition, the fairness index is also an important factor to be considered. Therefore, the time slot allocation problem is often a multi-objective optimization problem under (time slot) resource constraints.

At present, the research of time slot allocation problem focuses on problem modeling, including deterministic, uncertain, dynamic, multi-objective and other problems; on the other hand, research on optimization methods, such as integer programming, stochastic programming, robust Stick optimization, etc.

For the time slot allocation problem, the existing methods are insufficient for multi-objective optimization algorithms, and there are few studies on flight correlation under resource constraints. In view of the above deficiencies, the present invention proposes to use the Markov network to model the flight correlation, to perform search optimization through a distribution estimation algorithm based on probability statistics, and to provide a new fitness function for multi-objective optimization, so as to establish a A multi-objective optimization method for flight slot assignment.

SUMMARY OF THE INVENTION

Purpose of the invention: Under the premise of ensuring airport capacity constraints, the time slot is fully utilized to achieve multi-objective optimization of delay cost and fairness indicators, and to meet the computational time complexity requirements. In the optimization process, the time slot allocation solution is given, and the correlation between flights is also mined to provide data support for auxiliary decision-making.

Technical solution: a multi-objective optimization method for flight time slot allocation based on a distribution estimation algorithm, comprising the following steps:

Step 1, establish a multi-objective model of time slot allocation based on efficiency and fairness;

Step 2, based on the time slot allocation multi-objective model, establish a corresponding Markov network;

Step 3, establish a probability matrix;

Step 4, using the new fitness function, calculate and find the contemporary Pareto solution;

Step 5, according to the contemporary Pareto solution, update the trust solution set archive;

Step 6, learn the Markov network and update the probability matrix;

Step 7: Sampling according to the probability matrix and the Markov network to generate a new solution;

Step 8, the termination condition is reached, and the final result is obtained.

In the present invention, step 1 includes:

Step 1-1, establish the following multi-objective function:

where N is the total number of flights, M is the total number of time slots, c _ij is the delay cost when flight i is assigned to time slot j, x _ij is a decision variable, and F(x _ij ) is when flight i is assigned to time slot j Equation (1) represents the multi-objective function that minimizes delay cost and fairness difference index.

Step 1-2, set constraints:

c _ij =ω _i ×(t _j -SLDT _i ) (5)

t _j ≥SLDT _i , when x _ij =1 (6)

Among them, ω _i represents the importance weight of flight i in terms of delay time, ranging from 0.0 to 1.0; Q represents the total number of airlines, N _q represents the total number of flights of airline q, t _j represents the start time of time slot j; x _ij is a decision variable, representing whether flight i is allocated to time slot j, and can only take a value of 0 or 1; SLDT _i represents the scheduled landing time of flight i, c _i ^* represents the optimal delay of flight i, d( x _ij ) represents the difference between when flight i is assigned time slot j and the optimal situation, nd(x _ij ) represents the difference between the regularized d(x _ij ), and f _q (x _ij ) represents the difference of airline q Average regularized difference; formula (2) is the decision variable of the problem, formulas (3) and (4) represent the constraints of flight allocation and time slot allocation, respectively, and formula (5) is used to calculate flight delays Cost, formula (6) is the flight schedule constraint, formulas (7) and (8) calculate the difference between the actual time slot allocation and the optimal allocation and regularize it, and formulas (9) and (10) calculate all participating time slots. The variance of the gap-allocated companies is used as an indicator of airline fairness.

In the present invention, step 2 includes: establishing a Markov network corresponding to a multi-objective model of time slot allocation, and in the established Markov network, each node represents a decision variable, that is, the time slot allocation selection of flights; node X _i represents the time slot selection of flight i, and its value range D(X _i )={s ₁ , s ₂ ..., s _M }, where M is the total number of available time slots, and s _M represents the M-th available time slot; In the established Markov network, the connection between two nodes indicates that there is strong coupling between the two flights and competition for flight time slots.

In the present invention, step 3 includes: establishing the following probability matrix:

P(t) represents the probability matrix of the t-th generation in the iterative calculation process;

P _NM represents the value of the Nth row and Mth column of the matrix;

Row i of the matrix represents the probability that flight i is assigned to each time slot j, which adds up to 100%;

Column j of the matrix represents, for time slot j, the probability that each flight i is assigned to this time slot, which adds up to 100%;

When initializing the matrix, equal probability is used for assignment. For the specific problem of time slot allocation, there is a time constraint in formula (6), and the following rules are formulated: For each row i, in M time slots, the time slot will be earlier than the time of SLDT _i . The initial value of probability p of slot j is set to 0, and the initial value of probability p of time slot j later than SLDT _i is set to 1/z, where z is the number of time slots later than SLDT _i .

In the present invention, step 4 includes:

Step 4-1, the Pareto solution is divided into two parts: the boundary part and the middle part. The boundary part is the single-objective optimization part, and the middle part is the multi-objective optimization part. For the boundary part, the boundary fitness function in the VEGA algorithm is used to calculate:

EdgeFitness _f (X)=obj_value _f (X) (13)

where X represents an arbitrary solution, EdgeFitness _f (X) represents the fitness function of the solution X on the target f, and obj_value _f (X) represents the function value of the solution X on the target f;

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

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

The fitness function value of the middle part of all Pareto solutions is greater than or equal to 1, and at the same time, all non-Pareto solutions are less than 1 and greater than 0;

According to formulas (13) and (14), the boundary fitness and the intermediate fitness are calculated respectively;

In step 4-2, the following composite fitness function with exponential weight parameters is used to combine the boundary fitness and the intermediate fitness function:

where D-Fitness(X) is the comprehensive fitness of the solution X considering the composite middle part fitness and boundary part fitness, m is the total number of targets, N_POP is the total population, and Ranking _EdgeFitnessi (X) is the solution X based on the ith target The ranking of , Ranking _{CentralFitness} (X) is the ranking of X on the intermediate fitness, ω _i and ω _m+1 are the indexed weight parameters (for example, dual-objective optimization problem, m=2, so there are 3 indexed weight parameters : ω ₁ , ω ₂ , ω ₃ ), which are used to distinguish the scores of different rankings, and at the same time can realize the preference for a target.

In the present invention, step 5 includes: establishing an archive in the database for storing trust solutions, as a trust solution set archive, the capacity is the same as the population number N_POP in each iteration, or less than N_POP; for example, the value is N_POP×1/2 .

Calculate the fitness function value of all solutions of the current generation by formula (15), and select α% solutions from high to low according to the fitness function value and store them in the trust solution set archive. The degree function value from high to low selects 1-α% of the trust solutions and stores them in the current generation's trust solution set archive.

In the present invention, step 6 includes:

Step 6-1, learn Markov network:

By formula (16), the mutual information MI(X _i ; Y _j ) of the variables X _i and Y _j is calculated:

Among them, X _i and Y _j represent two random variables. In the time slot allocation problem, the variables specifically represent the time slot allocation of a flight; p( _xi |D) and p(y _j |D) represent The probability of variables X _i = _xi and Y _j =y _j of the set D of belief solutions, p(x _i ,y _j |D) is the joint probability of X _i = _xi and Y _j =y _j ;

If the mutual information of the two variables is greater than a given threshold threshold (for example, the value can be 0.5), it is determined that there is a strong correlation between the two variables, and a connection line is established in the Markov network and called neighbors to each other; the threshold threshold The threshold is set by the user, or dynamically updated according to the known mutual information during the calculation process, as shown in formula (17):

Among them, MI represents the mutual information of the two variables X _i and Y _j , NDV represents the total amount of variables, and the parameter β is used to control the structural complexity (usually can be set between 0.5-1.5);

Step 6-2, update the probability matrix:

The probability matrix of variable X _i in its neighbor N _i , as shown in equation (18):

where P(X _i |N _i ) represents the conditional probability matrix of variable X _i in all its neighbors N _i , {x _i,1 ,..., _xi,j ,..., _xi,Yi }respectively represent that variable X _i can The assigned value corresponds to the number of each time slot to be assigned in the time slot allocation, x _{i, Yi} represent the value that can be assigned to the variable X _i , and Y _{i represents the maximum number of the variable X i value} . , x _{i, Yi} is the maximum value that X _i can take, and in the time slot allocation problem, it is the latest time slot; p( _xi,Valuej |n _i,Neighbork ) represents the value of the variable X _i When Value _j , the conditional probability when all its neighbors jointly take the value scene is Neighbor _k ; Z _i represents the maximum number of the joint value scene of the neighbor set _Ni ;

For the probabilities in equations (12) and (18), first calculate the boundary probability P _t ^* (X=x) of the current generation:

Among them, N(X=x) represents the number of values x of the variable X in the trusted solution, and N(D) represents the total number of trusted solutions; the probability model is updated by formula (20):

P _t (X=x)=(1-λ)×P _t-1 (X=x)+λ×P _t ^* (X=x) (20)

Among them, P _t (X=x) and P _t-1 (X=x) represent the probability that the variable X takes the value of x in the t-th generation and the probability that the variable X takes the value of x in the t-1 generation in the iterative calculation process, respectively. , λ is the learning rate from the current generation, if λ=1, it means that the probability model will be completely learned from the current generation trust solution set archive;

Step 6-3, after the probability model is updated by the formula (20), the mutation operation of the formula (21) is performed by the mutation probability (for example, 5%):

P _t (X=x)=min(P _t (X=x)+θ,1-ε) (21)

where θ represents the variation value (eg 0.2) and ε is a small positive number (eg 0.01).

In the present invention, step 7 includes:

Step 7-1, randomly generate a feasible solution s=(s ₁ ,s ₂ ,…,s _n ), where s _n is the value of the last variable X _n , in the time slot allocation problem s ₁ ,s ₂ ,…, s _n represents the time slot selection of each flight, and s _n represents the time slot selection of the nth flight. The feasible solution s is generated as follows: Arrange all flights according to the flight schedule, starting from the first flight x ₁ , starting from all flights In the remaining time slots, randomly select a time slot s ₁ and assign it to x ₁ . It is required that s ₁ is not earlier than the planned departure time of x ₁ , and the time slot s ₁ is not selected by other flights. Delete the time slot s ₁ and move to the next One flight x ₂ , randomly selected in the same way, until all flight selection ends;

Step 7-2, randomly shuffle the order of the solution s=(s ₁ ,s ₂ ,...,s _n ) to obtain a new order s ₁ ^* ,s ₂ ^* ,...,s _n ^* , where _sn ^* represents a new order The time slot selection of the last flight in the flight sequence of ; it should be noted that: the corresponding relationship between the flight and the time slot at this time has not changed, but the order of the arrangement has changed;

Step 7-3, starting from s ₁ ^* , set s ₁ ^* corresponding to flight h, find out all the neighbors of flight h variable according to the established Markov network, and solve s=(s ₁ ,s ₂ , ...,s _n ) and the corresponding conditional probability p(x _h |N _h ) to determine the probability of selecting each time slot for flight h, and then according to the probability distribution, the roulette wheel selection algorithm (Roulette WheelSelection) is used, namely Randomly select according to the probability of various values, the higher the probability of a certain value, the greater the possibility of being selected) to generate a new value of s ₁ ^* ; then move to s ₂ ^* , and use the same method to calculate s ₂ ^* The new value of , until all variables are completed, the sampling of a solution is completed;

Step 7-4, repeat the method from step 7-1 to step 7-3, complete the sampling of N_POP solutions of a population, and finally form the solution set of the current generation.

Beneficial effect: Compared with the existing method, the advantages of the present invention are:

Using distribution estimation algorithm, compared with heuristic algorithms such as genetic algorithm, the global search ability and convergence speed are improved;

The Markov network structure is adopted to better describe and explore the relationship of variables and provide more auxiliary decision-making capabilities;

A fitness function based on exponential weight parameters is adopted to achieve low computational complexity and search flexibility.

Under the premise of ensuring airport capacity constraints, the method of the invention makes full use of time slots, realizes multi-objective optimization of delay cost and fairness index, and meets the requirements of computational time complexity. In the optimization process, the time slot allocation solution is given, and the correlation between flights is also mined to provide data support for auxiliary decision-making.

Description of drawings

The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments, and the advantages of the above or other aspects of the present invention will become clearer.

Figure 1 is the flow of the multi-objective distribution estimation algorithm based on Markov network.

Figure 2 is the flow of the distribution estimation algorithm.

Figure 3 is an example of a Markov network structure.

Figure 4 is a Pareto solution containing a boundary part and an intermediate part.

Figure 5 is the fitness function for the middle part.

Figure 6a is an example of a minimization problem.

Figure 6b is a function with an exponential weight parameter.

Figure 7 is an update archive candidate solution set.

FIG. 8 is the experimental verification result (part).

FIG. 9 is a flowchart of a multi-objective optimization method for flight slot allocation based on a distribution estimation algorithm.

Detailed ways

The present invention will be further described below with reference to the accompanying drawings and embodiments.

As shown in Figure 1, it is a schematic flowchart of the multi-objective distribution estimation algorithm based on Markov network. The algorithm is based on the Population-to-Population theory. Each generation generates new candidate solutions, and combines the trust solutions selected by the previous generation to perform multi-objective Pareto selection through the fitness function to update the trust solutions. Then, based on the newly obtained trust solution, the next generation of candidate solutions is obtained through probability model estimation and sampling in the distribution estimation algorithm, and the cycle continues until the termination condition occurs. It should be noted that the introduction of the Markov network is combined with the 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 adopted. The focus is on sampling to generate new solutions and estimating probabilistic models.

As shown in Figure 9, the concrete process of the present invention is as follows:

Step 1: Establish a multi-objective model of slot allocation based on efficiency and fairness.

(1) Parameters

Parameter name Parameter description

N total number of flights

M total number of slots

Q Total number of airlines

i flight i, i=1,...,N

j time slot j, j=1,...,M

q airline q, q=1,...,Q

Total number of flights from N _q airline q

c _ij delay cost when flight i is assigned to slot j

Is x _ij flight i assigned to slot j?

F(x _ij ) fairness index when flight i is assigned to slot j

t _j start time of slot j

Schedule of SLDT _i flight i planned landing time

c _i ^* flight i optimal-case delay

d(x _ij ) when flight i is assigned time slot j, the difference from the optimal case

nd(x _ij ) normalized d(x _ij ) difference

f _q (x _ij ) average regularized difference of airline q

ω _i The importance weight of flight i in terms of delay time

(2) Multi-objective function

(3) Constraints

c _ij =ω _i ×(t _j -SLDT _i ) (5)

t _j ≥SLDT _i , when x _ij =1 (6)

Among them, formula (1) represents the multi-objective function that minimizes the delay cost and fairness difference index. In formula (2), x _ij is a decision variable, representing whether flight i is allocated to time slot j, and can only take a value of 0 or 1. Equations (3) and (4) ensure that for any flight to be allocated, there is one and only one time slot allocated to it, and for any time slot to be allocated, at most one flight can be allocated. Equation (5) indicates that the delay cost is calculated according to the difference between the start time of the allocated time slot and the time of the schedule. Equation (6) guarantees that for any flight, the allocated time slot must be equal to or later than the schedule time. Equations (7) and (8) calculate and regularize the difference between the delay cost under the allocation scheme and the delay cost in the optimal case. Equations (9) and (10) calculate the variance of all participating companies in slot allocation as an airline fairness indicator.

Step 2: Based on the time slot allocation problem, establish the corresponding Markov network.

A Markov network is a meshed undirected graph (or bidirectional graph), which consists of a structure G and a parameter Ψ. Structure G contains node Node and connecting edge E, so a Markov network can be expressed as: MN={G,Ψ}={Node,E,Ψ}. As shown in Figure 3, the Markov network contains 5 variables (x ₁ -x ₅ ), and the connection between the variables indicates that there is a strong correlation between these two variables, and they are called neighbors. For example, variable X ₁ contains 2 neighbors X ₂ and X ₃ ; X ₂ contains 3 neighbors, X ₁ , X ₃ and X ₄ .

The value of variable X _i (D(X _i )={x _i ¹ ,...,x _i ⁿⁱ }), (x _i ⁿⁱ represents the maximum value that X _i can take) is calculated from the conditional probability of all its neighbors , as shown in formula (11).

p(x _i |x-{x _i })=p(x _i |N _i ) (11)

Among them, _{Ni is the set of neighbors of variable Xi, p(x i |} x-{x _i }) _represents the conditional probability of variable Xi under the joint value of all other variables, p( _xi |N _{i )} represents Conditional probability of variable X _i under the joint value of all its neighbors.

For the optimization problem of time slot allocation, this step is responsible for establishing the 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 the choice of time slot allocation for flights. Node X _i represents the time slot selection of flight i, and its value range is D(X _i )={s ₁ , s ₂ . . . , s _M }, that is, all time slots to be allocated. The connection between the two nodes indicates that there is a strong coupling and competition for resources (flight time slots) between the two flights. On the one hand, the function is to provide neighbor information for the evolution process and calculate the probability distribution; Gap exchange and other operations provide a certain basis.

Step 3: Build the probability matrix.

At this time, the probability matrix is based on the Markov network in the previous step. The probability model is shown in formula (12), where P(t) represents the name of the probability matrix of the t-th generation in the iterative calculation process, which is an N×M matrix. Row i represents the probability of flight i being assigned to each time slot j, which adds up to 100%. Viewed from column j, it represents the probability that for time slot j, each flight i is assigned to this time slot, and the total is also 100%.

During initialization, the assignment is generally performed with equal probability. For the specific problem of time slot allocation, there is a time constraint in formula (6). At the same time, in order to speed up the efficiency of sampling and searching, the following rules are formulated: For each row _i , in M time slots, the The initial value of probability p of time slot j is set to 0, and the initial value of probability p of time slot j later than SLDT _i is set to 1/Q, where Q is the number of time slots later than SLDT _i .

Step 4: Using the new fitness function, compute and find the contemporary Pareto solution.

The mathematical model contains two objectives that need to be optimized at the same time, and it is very important to design a set of efficient fitness functions to complete the solution search. In different multi-objective optimization algorithms, the fitness function plays a considerable role in the search efficiency and effect. The invention proposes a fitness function with an exponential weight parameter, which can improve the search efficiency and bring about higher flexibility.

As shown in Figure 4, taking the minimization multi-objective problem as an example, the Pareto solution is divided into two parts: the boundary part (single-objective optimization part) and the middle part (multi-objective optimization part). The boundary fitness function to calculate:

EdgeFitness _f (X)=obj_value _f (X) (13)

Where X represents an arbitrary solution, EdgeFitness _f (X) represents the fitness function of the solution X on the target f, and obj_value _f (X) represents the function value of the solution X on the target f.

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

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

It can be seen from Figure 5 that the fitness function given by formula (14) can well distinguish whether it belongs to a Pareto set. The CentralFitness values of all dominant solutions (Pareto solutions) are greater than or equal to 1, and the more dominant other solutions, the higher the value. Meanwhile, all non-dominated solutions are less than 1 and greater than 0. Therefore, a value of 1 can clearly distinguish whether a solution is the dominant solution.

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

In the next step, in order to obtain the final solution, it is necessary to comprehensively consider the above two fitness functions. There are two schemes, one is to use multi-island optimization, that is, some populations are optimized according to the boundary fitness function, and the other part of the population is optimized according to the intermediate fitness. This method has high execution efficiency, but the problem is that it is difficult to finally form a unified solution and the optimality of the Pareto solution is not high.

In the algorithm, a function with an exponential weight parameter is used to combine the boundary fitness and the intermediate fitness function, as shown in formula (15):

where D-Fitness(X) is the comprehensive fitness of the solution X considering the fitness of the composite middle part and the fitness of the boundary part, m is the total number of targets (dual target problem, m=2), N_POP is the total population, Ranking _EdgeFitnessi (X) is the ranking of X based on the ith target (the higher the EdgeFitness, the higher the ranking), the Ranking _{CentralFitness} (X) is the ranking of X in the middle fitness (the higher the CentralFitness, the higher the ranking), ω _i and ω _m+1 is an indexed weight parameter (for example, a dual-objective optimization problem, m=2, so there are 3 indexed weight parameters: ω ₁ , ω ₂ , ω ₃ ), which is used to distinguish the scores of different rankings. To achieve a preference for a certain goal.

The exponentiation weight parameter is illustrated by taking Fig. 6a and Fig. 6b as an example. The problem is assumed to be a minimization problem. In Figure 6a, the black solution and the white solution are considered. The black solution is a Pareto solution, while the white solution is not. The characteristic of the black solution is that it is very good in one dimension (boundary, f ₂ ). Although the white solution is slightly worse than the black solution in the dimension of the boundary f ₂ , it is better than the black solution in other dimensions. A Pareto solution only needs to be "excellent" in one dimension, not "good" in every dimension. Because in order to clearly distinguish and quickly find such a Pareto solution, Figure 6b is a function with an exponential weight parameter, that is, compared with the common linear relationship, the higher the ranking of a solution in a certain dimension, its The score (fitness function value) is much higher than the other solutions.

At the same time, the use of the exponential weight parameter also has a very important role: to achieve the flexibility of the search. Similar to the concept of weight, if the number of a certain dimension in the Pareto solution is insufficient, it is necessary to increase the exploration of this dimension. By increasing the corresponding exponential weight parameter, the dimension is focused on the calculation of the final fitness function value, "encouraging" this dimension to appear more excellent solutions, and improving the distribution of Pareto solutions.

Step 5: Update the trust solution set archive according to the contemporary Pareto solution.

As shown in Figure 7, based on the candidate solutions of the current generation (Candidate Solutions) and the trust solutions of the previous generation (Promising Data), through the composite fitness function in the previous step, some of the solutions in both parts become Pareto solutions. The trust solution archive (Archive) is used to store the updated trust solution. The fixed-size method is adopted in this algorithm, and its advantage is that there are more solutions for the subsequent probability estimation, which prevents the problem of premature maturity and falling into local optimum. For which candidate solutions are selected as non-Pareto solutions to enter the archive, in order to ensure the flexibility of search, the diversity of solutions and the stability of archives, the fitness function value of all solutions of the current generation is calculated by formula (15), And according to the fitness function value from high to low, select α% solution and save it into the trust solution set archive, and select (1-α%) trust solution from the previous generation trust solution set archive according to the fitness function value from high to low. Enter the trust resolution archive of this generation.

Step 6: Learn the Markov network and update the probability matrix.

The difference between this algorithm and the traditional distribution estimation algorithm is that the Markov network is introduced. Therefore, in the distribution estimation, in addition to the probability matrix needs to be updated, the structure of the Markov network also needs to be re-learned. The following is an explanation from two aspects:

(1) Learning Markov Networks

The connection between two nodes in the Markov network indicates that there is a strong correlation (non-causal relationship) between the two variables, and the assignment of the variable is also based on the "neighbors" associated with it, so Markov The learning process of the network is an important part of this algorithm.

One way to build the network structure is for problem experts to be responsible, but for most problems it is difficult to find experts who can give credible opinions on the correlations between variables. This algorithm uses a conditional independence test (ConditionalIndependence Test)-Mutual Information (MI) to estimate the structure of Markov network. The advantage of mutual information is that it is simple to use and fast to compute.

By formula (16), the mutual information MI(X _i ; Y _j ) of the variables X _i and Y _j is calculated:

Among them, X _i and Y _j represent two random variables (in the time slot allocation problem, the variables specifically represent the time slot allocation of a flight), and p(x _i |D) and p(y _j |D) are based on The probability of variables X _i = _xi and Y _j =y _j of the trusted solution set D, p(x _i ,y _j |D) is the joint probability of X _i = _xi and Y _j =y _j ;

If the mutual information of the two variables is greater than the given threshold threshold (for example, the value is 0.5), it is determined that there is a strong correlation between the two variables, and a connection line is established in the Markov network and called each other as neighbors; the threshold threshold is determined by the user Set, or dynamically update the threshold value according to the known mutual information during the calculation process, as shown in formula (17):

Among them, X _i and Y _j represent two random variables, MI represents the mutual information of the two variables X _i and Y _j , NDV represents the total amount of variables, and the parameter β is used to control the structural complexity (usually can be set to 0.5- 1.5). If β is set to a higher value, the threshold is higher, which means that there are fewer connections and shorter computation time in the Markov network, but some associated information is lost. If β is set to a lower value, the threshold becomes lower, which means that there are more connections in the Markov network and more related information is preserved, which helps to improve the quality of the search, but requires a longer time calculating time. Therefore, the parameter β can be used to control the search quality and efficiency.

Structural learning tends to take a long time, so for search efficiency, learning does not need to be done every generation. For example, it can be set to perform structure learning every 10 generations (probability matrix needs to be updated every generation).

(2) Update the probability matrix

In Equation (12), the probability matrix is expressed as the independent probability of each variable. At the same time, due to the existence of correlation in the Markov network, it is also necessary to establish and update the conditional probability matrix based on neighbors. The probability matrix of variable X _i in its neighbor N _i , as shown in formula (18),

表示变量X_i的某取值Value_j时，在其所有邻居联合取值场景为Neighbor_k下的条件概率，Z_i表示邻居集合N_i的联合取值场景的最大编号。where P(X _i |N _i ) represents the conditional probability matrix of variable X _i in all its neighbors N _i , {x _i,1 ,..., _xi,j ,..., _xi,Yi }respectively represent that variable X _i can The _{assigned value corresponds} to the number of each to-be-allocated time _slot in the time slot allocation. In the allocation problem, it is the latest time slot);

Represents the conditional probability of a value _j of the variable X _i when the joint value scene of all its neighbors is Neighbor _k , and Z _i represents the maximum number of the joint value scene of the neighbor set _Ni .

For the probabilities in formulas (12) and (18), both are updated according to the trusted solution obtained by calculation or the content in the archive of the trusted solution set. Since it is an iterative search, first calculate the boundary probability of the current generation:

Among them, N(X=x) represents the number of values x of the variable X in the trusted solution, N(D) represents the total number of trusted solutions, and the setting of 1/|X| is to prevent the possibility of 0 from appearing. Probability calculations have no real impact.

After calculating the probability of the current generation, the probability model needs to be updated. In order to ensure the stability of the model, the calculation is performed by formula (20):

P _t (X=x)=(1-λ)×P _t-1 (X=x)+λ×P _t ^* (X=x) (20)

Among them, P _t (X=x) and P _t-1 (X=x) represent the probability that the variable X takes the value of x in the t-th generation and the probability that the variable X takes the value of x in the t-1 generation in the iterative calculation process, respectively. , λ is the learning rate from the current generation, if λ = 1, it means that the probability model will be completely learned from the trusted solutions obtained from the current generation. Considering that the filing mechanism has been designed in the algorithm, the value of λ can be set larger than 0.5.

In order to further increase the diversity of solutions and avoid falling into local optimum, a mutation operation similar to GA is also added to the calculation process of probability update. After Equation (20) updates the probability, perform the mutation operation of Equation (21) with a mutation probability (eg, 5%):

P _t (X=x)=min(P _t (X=x)+θ,1-ε) (21)

Where θ represents the variation value (eg 0.2), ε is a small positive number (eg 0.01), in order to ensure that all probabilities are less than 100%.

Step 7: Sampling according to the probability matrix and Markov network to generate a new solution.

The probabilistic model and network model have been obtained in the previous step, and the next step is to generate the next generation of solutions. Markov network is very different from Bayesian network (BN), and there is no causal relationship, so there is no clear order in the process of sampling to generate new solutions. In order to obtain the new solution, the algorithm adopts a Monte Carlo method based on Markov chain - Gibbs Sampler. The specific process is as follows:

(1) Randomly generate a feasible solution s=(s ₁ , s ₂ ,...,s _n ), where s _n is the value of the last variable X _n (in the time slot assignment problem, it represents the time slot selection of each flight, s _n represents the time slot selection of the nth flight). The generation method is as follows: Arrange all flights according to the flight schedule, starting from the first flight x ₁ , from all remaining time slots, randomly select a time slot s ₁ (requires that s ₁ is not earlier than the planned departure time of x ₁ , and time slot s ₁ is not selected by other flights) assigned to x ₁ , delete time slot s ₁ , move to the next flight x ₂ , select randomly in the same way, until all flight selection ends.

(2) Randomly shuffle the order of the solution s=(s ₁ , s ₂ ,...,s _n ) to obtain a new order s ₁ ^* ,s ₂ ^* ,...,s _n ^* , where s _n ^* represents a new flight Slot selection for the last flight in the sequence. (It should be noted that the correspondence between flights and time slots at this time has not changed, but the order of arrangement has changed).

(3) Starting from s ₁ ^* (assuming the corresponding flight h), find out all the neighbors of the flight h variable according to the established Markov network, and solve s=(s ₁ ,s ₂ ,...,s _n ) according to this time. The value of the neighbors and the corresponding conditional probability p(x _h |N _h ) determine the probability of selecting each time slot for flight h, and then according to the probability distribution, the roulette wheel selection algorithm (Roulette WheelSelection) is used. The probability of s is randomly selected, the higher the probability of a value, the greater the possibility of being selected) to generate a new value of s ₁ ^* ; then move to s ₂ ^* , and use the same method to calculate the new value of s ₂ ^* until Completion of all variables completes the sampling of a solution.

Repeat the above steps (1) to (3) to complete the sampling of N_POP solutions in a population, and finally form the solution set of the current generation.

Step 8: The termination condition is reached and the final result is obtained.

Repeat the above steps until the termination condition is met. The termination condition can be set to the number of iterations, such as 1000 generations. When the termination condition is reached, the Pareto solution obtained at this time is the final result.

In order to verify the correctness of the model and algorithm in the method of the present invention, based on the traffic management environment in North China, the simulated flight data for a total of 6 hours from 18:00-24:00 is selected for verification, and it is assumed that 22:00-23:00 The Beijing Capital Airport (ZBAA) experienced a capacity drop and overcapacity.

The results of flight slot allocation are shown in Figure 8 (partial results). The results in the graph give the estimated and calculated departure times of the affected flights, as well as suggestions for ground delay times. It can be seen intuitively from the results that, on the one hand, some flights of companies with more flights do not need to be delayed, indicating that they meet the efficiency (delay time) index in multi-objective optimization. On the other hand, the average delay time can also be found. Similarly, there is no excessive delay of a certain airline's flight, which also proves that the fairness index is guaranteed. The experimental simulation and verification show that the method can effectively solve the time slot allocation problem, and obtain satisfactory results according to the multi-objective optimization calculation of problem modeling.

The present invention provides a multi-objective optimization method for flight time slot allocation based on a distribution estimation algorithm. There are many specific methods and approaches for realizing the technical solution. The above are only the preferred embodiments of the present invention. For those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and modifications can also be made, and these improvements and modifications should also be regarded as the protection scope of the present invention. All components not specified in this embodiment can be implemented by existing technologies.

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