Cooperative interval management conversion program optimization method based on flight demand uncertainty

Abstract

The invention provides a collaborative interval management conversion program optimization method based on flight demand uncertainty. The method comprises the following steps: determining an uncertainty variable of a collaborative interval management conversion program CMCP optimization problem; estimating probability density function key parameters of normal distribution of the uncertainty variables according to historical data; constructing a CMCP opportunity constrained programming model and equivalent transformation thereof; and solving the converted CMCP opportunity constraint programming model based on a particle swarm algorithm. According to the method, the flight demand uncertainty factor, namely the uncertainty of the flight point passing time, is considered. Under the condition that the actual flight passing time is uncertain and interfered, the flight operation based on the optimization result can keep good airspace performance, so that the actual execution effect of the flight plan cannot be seriously fluctuated due to the uncertainty of the flight passing time.

Description

Technical Field

[0001] The present invention belongs to the technical field of flight resource allocation optimization in air traffic flow management, and in particular relates to a collaborative interval management conversion program optimization method based on uncertain flight demand. Background Art

[0002] Miles / Minutes-In-Trail (MIT / MINIT) is an air traffic flow management procedure that requires the distance between two consecutive aircraft to meet specific interval requirements at designated waypoints, sector boundaries, etc. Given the advantages of simple implementation, MIT / MINIT has always been the main tool for cross-border air traffic flow management. However, with the increasing busyness of air transportation, MIT / MINIT has gradually exposed typical problems such as high air delay costs, poor delay predictability, and increased workload for controllers and pilots. As an emerging air traffic flow management procedure, CMCP is the core of airspace resource allocation. However, at present, CMCP does not consider the uncertainty of flight transit time when allocating resources, resulting in the problem of execution deviation in the actual flight resource allocation results. For example, if a certain transit time allocated for a flight does not arrive on time in reality, the actual interval between aircraft will be less than the specified safety interval, resulting in certain safety risks.

[0003] In the technical field of air traffic flow management flight resource allocation optimization, uncertainty factors can be divided into two categories: capacity uncertainty and demand uncertainty. Previous research on the optimization of air traffic flow management tools mainly considered the uncertainty factors of airspace capacity, thereby constructing dynamic models or air speed adjustment optimization schemes to cope with scenarios where airspace capacity changes over time. As for the consideration of demand uncertainty factors, current work has mainly focused on airport surface operation management, such as flight departure time optimization or flight arrival and departure coordinated sequencing optimization considering the uncertainty of flight surface operation time. In the air operation stage, existing literature has conducted statistical analysis on the uncertainty of flight arrival time in each sector, and obtained the distribution characteristic function of uncertainty factors, but no optimization work has been carried out. Summary of the invention

[0004] Purpose of the invention: The technical problem to be solved by the present invention is to provide a method for optimizing the conversion procedure of collaborative interval management based on uncertain flight demand in view of the shortcomings of the prior art, including the following steps:

[0005] Step 1, determine the uncertainty variables of the CMCP optimization problem of the collaborative interval management conversion program;

[0006] Step 2: Estimate the key parameters of the probability density function of the normal distribution of the uncertainty variable based on historical data;

[0007] Step 3: Construct the collaborative interval management transformation program CMCP opportunity constraint programming model and perform equivalent transformation;

[0008] Step 4: Solve the opportunity-constrained programming model of the transformed collaborative interval management transformation program CMCP based on the particle swarm algorithm.

[0009] Step 1 includes: The uncertainty variable in the collaborative interval management conversion program CMCP optimization problem is the time when flight f passes the kth entry point of the sector and the time of exit At the same time, the flight was originally scheduled to pass the kth entry point in the sector and the time of exit is a certain standard time; uncertainty variable and The calculation formula is:

[0010]

[0011] in, and are random variables, representing the uncertainty time deviation of the kth sector entry point and the uncertainty time deviation of the sector exit point.

[0012] Step 2 includes:

[0013] Step 2-1: Randomly select flight trajectory data of flights that passed through the sector entry or exit points on a certain day or multiple consecutive days in history, thereby forming a historical flight data set F sta ;

[0014] Step 2-2, calculate the set F sta The actual time when each flight passes through the sector k and the planned time for the kth entry point through the sector The difference The actual time of the sector exit point and the planned sector exit time The difference

[0015]

[0016] Step 2-3: Filter out the set of flights that fly through the kth entry point in the sector Statistics of corresponding flights Complete the normal distribution fitting of the statistical data and estimate the random variables of all n entry points in the sector The mean μ in the normal distribution probability density function in,k and standard deviation

[0017] Step 2-4: Filter out the flight set that passes through the sector exit point Statistics of corresponding flights Complete the normal distribution fitting of statistical data and estimate the random variables The mean μ in the normal distribution probability density function out and standard deviation σ out .

[0018] Step 3 includes:

[0019] Step 3-1: According to the start and end time of the collaborative interval management conversion program CMCP and the sector object of the restriction constraint, filter out the following related flight sets:

[0020] The flight set F affected by CMCP;

[0021] The set of flights F that are affected by CMCP and fly through the kth entry point in the sector in,k ;

[0022] The set of flights F that are affected by CMCP and fly through the sector exit point out ;

[0023] The set of flights affected by CMCP and subject to ground delays F G ;

[0024] The set of flights F that are affected by CMCP and controlled by interval management M ;

[0025] Set F M The set of flights F that fly through the kth entry point in the sector M,k ;

[0026] Step 3-2: Construct the collaborative interval management conversion program CMCP. The decision variable of the opportunity-constrained programming model is the delay time d assigned to flight f. f ;

[0027] Step 3-3, construct the collaborative interval management conversion program CMCP. The optimization objectives of the opportunity-constrained programming model are efficiency objectives and fairness objectives, where the efficiency objective is to minimize the total flight delay cost, and the fairness objective is to minimize the maximum airline fairness ratio deviation;

[0028] Step 3-4, construct the constraints of the collaborative interval management conversion program CMCP opportunity constraint programming model;

[0029] Step 3-5: Based on the opportunity constraint equivalent transformation theory, the opportunity constraints in the collaborative interval management transformation program CMCP opportunity constraint programming model are equivalently transformed.

[0030] In step 3-3, the total flight delay cost is calculated as follows:

[0031]

[0032] Where C represents the total flight delay cost, and β represents the ratio of the cost per minute of air delay to the cost per minute of ground delay.

[0033] In step 3-3, the maximum airline fair ratio deviation is calculated as follows:

[0034]

[0035] Where E represents the maximum fair proportion deviation of airlines, A represents the set of airlines to which the affected flights belong, l represents the total number of airlines in the set ρ, and ρ a represents the fair proportion of airline a, ρ a The calculation method is:

[0036]

[0037] Among them, C a represents the total flight delay cost of airline a, q a represents the total number of affected flights belonging to airline a, and q represents the total number of affected flights.

[0038] In step 3-4, the constraints of the collaborative interval management conversion program CMCP opportunity constrained programming model are:

[0039]

[0040] Among them, Pr{} represents the probability of satisfying the conditional constraints in the curly brackets {}, and f ′ indicates another flight different from flight f, represents the time interval requirement between two consecutive flights entering from the kth entry point of the sector, α in,k represents flight f and flight f ′ The confidence level of the constraint on the interval between the i-th entry point of the sector is maintained, s out Indicates the time interval requirement between two consecutive flights leaving from the sector exit point, α out represents flight f and flight f ′ The confidence level at which the sector exit point interval constraint holds, and N represents a set of natural numbers.

[0041] In step 3-5, the result of equivalent transformation of the opportunity constraints in the collaborative interval management transformation program CMCP opportunity constraint programming model is:

[0042]

[0043] in, and Δs out The calculation method of is as follows:

[0044]

[0045] in, It represents the additional interval time that needs to be maintained between two consecutive flights entering from the kth entry point of the sector, Δs out Indicates the additional interval time between two consecutive flights departing from the sector exit point. Represents the inverse of the standard normal distribution function.

[0046] Step 4 includes:

[0047] Step 4-1, initialize the first generation particle swarm: set j = 1, set the total number of particles to n p , the maximum number of iterations is j max , randomly generate a particle swarm p j The initial positions and initial velocities of all particles in , where the subscript j represents the round of iterative calculation in the algorithm;

[0048] Step 4-2, allocate flight delays: use the delay allocation operator as the particle swarm p j Each particle in the Corresponding to a set of feasible solutions, where the delay allocation operator refers to the position of particles one by one in order based on the greedy algorithm The flights in the distribution of each point of the earliest available time until the particle position All flights complete the transit time allocation and form a set of resource allocation schemes, that is, a set of feasible solutions;

[0049] Step 4-3, calculate the objective function value: according to the resource allocation plan after the delay allocation operator acts, calculate the particle swarm p j The total flight delay cost and the maximum airline fairness ratio deviation are the two optimization target values ​​for each particle;

[0050] Step 4-4, update the elite set: according to the objective function value of each particle, calculate the dominance level and crowding degree of each particle, and then based on the elite retention strategy, update the elite set from the particle group p j Select n e =ξ·n p Excellent particles as elite set Among them, ξ is the preset ratio threshold parameter;

[0051] Step 4-5: Implement game operation: According to the elite retention strategy, from the particle swarm p jSelect half of the particles as the parent particle group Then, from the parent particle swarm Take out the i-th particle to be updated From the elite collection Randomly select two candidate elite particles from the , calculate the vector angles θ1 and θ2 between the two candidate elite particles and the current i-th particle to be updated, and use the candidate elite particle with the smaller angle as the guide particle of the current i-th particle to be updated. The position of the guide particle is recorded as Among them, the i-th particle to be updated and the z-th particle The vector angle θ of the candidate elite particle z The calculation formula is:

[0052]

[0053] in, and They represent the efficiency target value and fairness target value of the i-th particle to be updated, respectively. and They represent the efficiency target value and fairness target value of the zth candidate elite particle respectively;

[0054] Step 4-6, Update particle information: parent particle group The speed and position of each particle in the Update to obtain the descendant particle swarm

[0055] Step 4-7, merge the parent particle swarm: and Merge to generate a new generation of particle swarm p j+1 , when the number of iterations meets the set maximum number of iterations j nax Then, go to step 4-8, otherwise, let j=j+1 and go to step 4-2;

[0056] Step 4-8: Output Pareto front solution set: for particle swarm p j+1 The delay allocation, objective function calculation and dominance ranking are redone to determine the Pareto front solution set of the current latest generation of particle swarm.

[0057] In steps 4-6, the parent particle swarm The speed and position of each particle in the To update, the formula is:

[0058]

[0059] in, represents the particle swarm p j The position of the i-th particle in represents the particle swarm p j+1 The position of the i-th particle in i j represents the particle swarm p j The velocity of the ith particle in i j+1 represents the particle swarm p j+1 is the velocity of the ith particle in , r1 and r2 are non-negative constants, η1 and η2 are random numbers in [0,1], The operation is used to subtract the positions of two particles, and the result is the flight speed of the generated particle; The operation is used to add the particle position and particle velocity, and the result is the new position of the added particle; The operation is used to merge two particle speeds, and the result is to generate a new particle speed; The operation is used to multiply the constant by the particle velocity. The result is the generation of a new particle velocity. The operation form is based on the constant λ C There are three different cases for different values ​​of :

[0060] When the constant λ C = 0, the new particle velocity after multiplication is an empty set;

[0061] When 0<λ c <1, the new particle speed after multiplication is the original particle speed The remainder of the commutators, where h represents the length of the multiplied velocity, Represents h and λ C The integer rounded up after multiplication;

[0062] When C ≥1, C Split into two new constants λ C,1 and λ C,2 , satisfying λ C =λ C,1 +λ C,2 , and λ C,1 is a non-negative natural number, 0≤λ C,2 <1, so the new particle speed after multiplication is λ C,1 The speed of the original particles is connected end to end and then merged into λ C,2 The result after multiplying by the original particle velocity is

[0063] In summary, the present invention aims at the uncertainty factor of the flight's transit time in the air (which belongs to the concept of flight demand uncertainty), and adopts the chance constraint programming method to optimize the flight resource allocation in the CMCP program, so as to achieve a more feasible resource allocation result.

[0064] Beneficial effects: For the CMCP optimization problem, the present invention takes into account the uncertainty factor of flight demand, that is, the uncertainty of flight passing time. In the case of uncertainty interference of the actual flight passing time, compared with the results of other optimization methods that do not consider uncertainty, the flight operation based on the optimization results of the present invention can maintain good airspace performance, so that the actual execution effect of the flight plan will not fluctuate seriously due to the uncertainty of the flight passing time. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments, and the above and / or other advantages of the present invention will become more clear.

[0066] Figure 1 It is a schematic diagram of the overall steps of CMCP optimization of the present invention.

[0067] Figure 2 It is a description diagram of the CMCP optimization problem of the present invention.

[0068] Figure 3 It is a schematic diagram of uncertainty variables of the present invention.

[0069] Figure 4 1 is a diagram for analyzing the flight transit time error at the sector entry point U1 according to an embodiment of the present invention.

[0070] Figure 5 1 is a diagram for analyzing the flight transit time error at the sector entry point U2 according to an embodiment of the present invention.

[0071] Figure 6 4 is a diagram of the flight transit time error analysis of the fan-out point D according to an embodiment of the present invention.

[0072] Figure 7 Schematic diagram of the Pareto front solution set of the calculation result of the embodiment of the present invention. DETAILED DESCRIPTION

[0073] In the embodiment of the present invention, the optimization problem of the collaborative interval management conversion program (CMCP) is first described (eg Figure 2 shown):

[0074] For sector B, due to the limited airspace capacity of its downstream traffic management unit, sector B receives the interval restriction s issued by the downstream unit to its exit point D. out To control the air delay time in sector B, sector B needs to k (k∈{1,2,…,n}) Set appropriate interval constraints In addition, since sector B has control agreements with some upstream airports (these airports can be called cooperative airports), sector B can directly control the take-off time of flights departing from the cooperative airports. Therefore, in the coordinated interval management conversion procedure CMCP, there are two types of flights, namely flights subject to ground delay control (referred to as GDP flights, such as Figure 2 f G ) and flights subject to interval management control (referred to as MINIT flights, such as Figure 2 f M GDP flights consume the delay by waiting on the ground, so these flights have the right of exemption when entering sector B, that is, they are not subject to the entry point U k The interval limit of the flight only needs to meet the interval limit of the exit point D, while the MINIT flight needs to meet the interval limit of the entry point U at the same time. k And the interval limit of the fan-out point D.

[0075] The embodiment of the present invention discloses a method for optimizing a collaborative interval management conversion program based on uncertain flight demand, which specifically includes the following steps (eg Figure 1 shown):

[0076] Step 1, determine the uncertainty variables of the collaborative interval management conversion program (CMCP) optimization problem;

[0077] Step 2: Estimate the key parameters of the probability density function of the normal distribution of the uncertainty variable based on historical data;

[0078] Step 3: Construct the CMCP opportunity-constrained programming model and its equivalent transformation;

[0079] Step 4: Solve the transformed CMCP opportunity constrained programming model based on the particle swarm algorithm.

[0080] Step 1 includes: the uncertainty variable in the CMCP optimization problem is the time when flight f passes the kth entry point of the sector and the time of exit At the same time, the flight was originally scheduled to pass the kth entry point in the sector and the time of exit The uncertainty variable is the determined standard time and The calculation method of is shown in formula (1) and formula (2):

[0081]

[0082] in, and are random variables, representing the uncertainty time deviation of the kth sector entry point and sector exit point, such as Figure 3 shown.

[0083] In step 2, the specific steps of estimating key parameters of the probability density function of the normal distribution of the uncertainty variable based on historical data include:

[0084] Step 2-1: Randomly select flight trajectory data of flights that passed through the sector entry or exit points on a certain day or multiple consecutive days in history, thereby forming a historical flight data set F sta ;

[0085] Step 2-2, calculate the set F sta The actual time when each flight passes through the sector k and the planned time for the kth entry point through the sector The difference The actual time of the sector exit point and the planned sector exit time The difference and The calculations are shown in formula (3) and formula (4) respectively;

[0086]

[0087]

[0088] Step 2-3: Filter out the set of flights that fly through the kth entry point in the sector Statistics of corresponding flights The normal distribution fitting of statistical data is completed through the normfit function in MATLAB, and the random variables of all n entry points in the sector are estimated. The mean μ in the normal distribution probability density function in,k and standard deviation σ in,k ;

[0089] Step 2-4: Filter out the flight set that passes through the sector exit point Statistics of corresponding flights The normfit function in MATLAB is used to fit the normal distribution of statistical data and estimate the uncertainty variables. The mean μ in the normal distribution probability density function out and standard deviation σ out ;

[0090] In step 3, the specific steps of constructing the CMCP opportunity constrained programming model and its equivalent transformation are:

[0091] Step 3-1: According to the start and end time of the collaborative interval management conversion program CMCP and the restricted sector object, filter and select the following related flight sets: (1) the flight set F affected by CMCP, (2) the flight set F affected by CMCP and flying through the kth entry point of the sector in,k , (3) The set of flights F that are affected by CMCP and fly through the sector exit point out , (4) The set of flights F affected by CMCP and controlled by ground delays G , (5) The set of flights F that are affected by CMCP and controlled by interval management M , (6) Set F M The set of flights F that fly through the kth entry point in the sector M,k ;

[0092] Step 3-2: Construct the CMCP opportunity-constrained programming model with the decision variable being the delay time d assigned to flight f f ;

[0093] Step 3-3: The optimization objectives of the CMCP opportunity-constrained programming model are efficiency objectives and fairness objectives. The efficiency objective is to minimize the total flight delay cost (the total flight delay cost is calculated as shown in formula (5)), and the fairness objective is to minimize the maximum airline fairness ratio deviation (the maximum airline fairness ratio deviation is calculated as shown in formula (6)):

[0094]

[0095] In formula (5), C represents the total flight delay cost, β represents the ratio of the cost of each minute of air delay to the cost of each minute of ground delay, and in formula (6), E represents the maximum airline fair proportion deviation, A represents the set of airlines to which the affected flights belong, and ρ a represents the fair proportion of airline a (calculated as shown in formula (7)), l represents the total number of airlines in set A, and in formula (7), C a represents the total flight delay cost of airline a, q a represents the total number of affected flights belonging to airline a, and q represents the total number of affected flights;

[0096] Step 3-4: Construct the constraint conditions of the CMCP opportunity constraint programming model, as shown in formula (8) to formula (10), where formula (8) and formula (9) are opportunity constraints:

[0097]

[0098]

[0099] In formula (8) and formula (9), Pr means probability, which means the probability of satisfying the conditional constraints in the curly brackets behind it, and f ′ indicates another flight different from flight f, represents the time interval requirement between two consecutive flights entering from the kth entry point of the sector, α in,k represents flight f and flight f ′ The confidence level of the constraint on the interval between the i-th entry point of the sector is maintained, s out Indicates the time interval requirement between two consecutive flights leaving from the sector exit point, α out represents flight f and flight f ′ The confidence level of maintaining the sector exit point interval constraint is established. In formula (10), N represents a natural number set;

[0100] Step 3-5: Based on the opportunity constraint equivalent transformation theory, the opportunity constraints in the CMCP opportunity constraint programming model are equivalently transformed. Formula (8) and formula (9) are transformed into formula (11) and formula (12):

[0101]

[0102] In formula (11), represents the additional interval time that needs to be maintained between two consecutive flights entering from the kth entry point of the sector. The calculation is shown in formula (13). In formula (12), Δs out It represents the additional interval time between two consecutive flights leaving from the sector exit point. The calculation is shown in formula (14). In formula (13) and formula (14), Represents the inverse function of the standard normal distribution function;

[0103] Step 4: The specific steps of solving the transformed CMCP opportunity constraint programming model based on the particle swarm algorithm are as follows:

[0104] Explain the relevant symbols appearing in step 4. Suppose the particle swarm is The particle position Expressed as a flight sequence: The particle velocity V i j Expressed as a basic commutative sequence: V i j =(s1,s2,…,s v ,…,s h ). Among them, n p represents the total number of particles in the particle swarm, q represents the total number of affected flights, h represents the length of the basic exchange sequence, and f x Indicates the flight number that ranks xth in the flight sequence, sv represents the commutator ranked in the vth position in the basic commutation sequence. Furthermore, the commutator s v Note: s v =s(x,y) means swapping the x-th flight with the y-th flight in the flight sequence;

[0105] Step 4-1, initialize the first generation particle swarm: set j = 1, set the total number of particles to n p , the maximum number of iterations is j max , randomly generate a particle swarm p j The initial positions and initial velocities of all particles in , where j represents the number of rounds of iterative calculations in the algorithm;

[0106] Step 4-2, allocate flight delays: use the delay allocation operator as the particle swarm p j Each particle in the Corresponding to a set of feasible solutions, where the delay allocation operator refers to the position of particles one by one in order based on the greedy algorithm The flights in the distribution of each point of the earliest available time until the particle position All flights complete the transit time allocation and form a set of resource allocation schemes, that is, a set of feasible solutions;

[0107] Step 4-3, calculate the objective function value: according to the resource allocation plan after the delay allocation operator acts, calculate the particle swarm p j The total flight delay cost and the maximum airline fairness ratio deviation are the two optimization target values ​​for each particle;

[0108] Step 4-4, update the elite set: according to the objective function value of each particle, calculate the dominance level and crowding degree of each particle, and then based on the elite retention strategy, update the elite set from the particle group p j Select n e =ξ·n p Excellent particles as elite set Among them, ξ is the preset ratio threshold parameter. The calculation method of dominance level, crowding degree and elite retention strategy can refer to the paper "A fast and elitist multi objective genetic algorithm: NSGA-II" by K.DEB et al.;

[0109] Step 4-5: Implement game operation: According to the elite retention strategy, from p j Select half of the particles as the parent particle group Then, from the parent particle swarm Take out the i-th particle to be updated From the elite collection Randomly select two candidate elite particles from the , calculate the vector angles θ1 and θ2 between the two candidate elite particles and the current i-th particle to be updated, and use the candidate elite particle with the smaller angle as the guide particle of the current i-th particle to be updated. The position of the guide particle is recorded as Among them, The vector angle θ between the i-th particle to be updated and the z-th candidate elite particle z The calculation of is shown in formula (15):

[0110]

[0111] In formula (15), and They represent the efficiency and fairness target values ​​of the i-th particle to be updated, respectively. and They represent the efficiency and fairness target values ​​of the zth candidate elite particle respectively;

[0112] Step 4-6, update particle information: The speed and position of each particle in the Update to obtain the descendant particle swarm The update is shown in formula (16) and formula (17):

[0113]

[0114] In formula (16), represents the particle swarm p j The position of the i-th particle in represents the particle swarm p j+1 The position of the i-th particle in i j represents the particle swarm p j The velocity of the ith particle in i j+1 represents the particle swarm p j+1 The velocity of the ith particle in the equation, r1 and r2 are non-negative constants, η1 and η2 are random numbers in the range [0,1]. Further, in equations (16) and (17), we define a new operation: and as follows:

[0115] The operation is used to "subtract" the positions of two particles, and the result is the flight speed of the generated particles;

[0116] The operation is used to "add" the particle position and velocity, and the result is the new position of the particle.

[0117] The operation is used to merge two particle velocities. The result is to generate a new particle velocity. In particular, when one of the two merged particle velocities is an empty set, the result is to keep the other particle velocity unchanged.

[0118] The operation is used to "multiply" the constant by the particle velocity. The result is to generate a new particle velocity. The operation is divided into three cases according to the value of the constant (assuming the constant is λ C ):

[0119] (1) When λ C = 0, the new velocity of the particle after "multiplication" is an empty set;

[0120] (2) When 0<λ C <1, the new velocity of the particle after "multiplication" is equal to the original velocity of the particle. The remainder of the commutators, where h represents the length of the multiplied velocity (basic commutative sequence), Represents h and λ C The integer rounded up after multiplication;

[0121] (3) When λ C ≥1, C Split into two new constants λ C,1 and λ C,2 , satisfying λ C =λ C,1 +λ C,2 , and λ C,1 is a non-negative natural number, 0≤λ C,2 <1, so the new particle speed after multiplication is λ C,1 The speed of the original particles is connected end to end and then merged into λ C,2 The result after multiplying by the original particle velocity is

[0122] Step 4-7, merge the parent particle swarm: and Merge to generate a new generation of particle swarm p j+1 , when the number of iterations meets the set maximum number of iterations j max Then, go to step 4-8, otherwise, let j=j+1 and go to step 4-2;

[0123] Step 4-8: Output Pareto front solution set: for particle swarm p j+1 The delay allocation, objective function calculation and dominance ranking are redone to determine the Pareto front solution set of the current latest generation of particle swarm.

[0124] In this embodiment, a simulation example is designed to illustrate the feasibility of the method of this embodiment:

[0125] The simulation environment adopts Sanya Flight Information Region. The example sets two upstream entry points: U1 and U2, and one downstream exit point. The number of MINIT flights passing through U1 is 13, the number of MINIT flights passing through U2 is 15, the number of GDP flights is 24, and the total number of flights is 52.

[0126] In order to obtain the probability distribution function of the uncertainty variable, the flight data of the flights passing through the entry points U1 and U2 and the exit point D in the Sanya Flight Information Region on March 1, 2019 were analyzed. The frequency histogram of the flight transit time error is shown in Figure 2. Figure 4 , Figure 5 and Figure 6 The results of normal distribution fitting all pass the Jarque-Bera test (such as Figure 4 , Figure 5 and Figure 6 ).

[0127] According to the data analysis results, σ in,1 =352.4,σ in,2 =350.4,σ out =384.5, take the confidence level α in,1 =α in,2 =α out =0.8, then the model parameters are: Δs out =7, and set the remaining model parameters as follows: s out =8, β=2, and set the algorithm parameters as follows: n p =100,j max =500, ξ=0.2, r1=1, r2=1. According to the steps of the method of this embodiment, the Pareto front solution set of the example is calculated, such as Figure 7 shown.

[0128] The present invention provides a method for optimizing the conversion program of collaborative interval management based on uncertain flight demand. There are many methods and ways to implement the technical solution. The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. 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

1. A collaborative interval management conversion program optimization method based on flight demand uncertainty, characterized in that: The following steps are involved: Step 1, determine the uncertainty variables of the CMCP optimization problem of the collaborative interval management conversion program; Step 2: Estimate the key parameters of the probability density function of the normal distribution of the uncertainty variable based on historical data; Step 3: Construct the collaborative interval management transformation program CMCP opportunity constraint programming model and perform equivalent transformation; Step 4: Solve the opportunity-constrained programming model of the transformed collaborative interval management transformation program CMCP based on the particle swarm algorithm.

2. The method according to claim 1, characterized in that Step 1 includes: The uncertainty variable in the collaborative interval management conversion program CMCP optimization problem is the time when flight f passes the kth entry point of the sector and the time of exit At the same time, the flight was originally scheduled to pass the kth entry point in the sector and the time of exit is a certain standard time; uncertainty variable and The calculation formula is: in, and are random variables, representing the uncertainty time deviation of the kth sector entry point and the uncertainty time deviation of the sector exit point.

3. The method according to claim 2, characterized in that Step 2 includes: Step 2-1: Randomly select flight trajectory data of flights that pass through the sector entry or exit points, thereby forming a historical flight data set F sta ; Step 2-2, calculate the set F sta The actual time when each flight passes through the sector k and the planned time for the kth entry point through the sector The difference The actual time of the sector exit point and the planned sector exit time The difference Step 2-3: Filter out the set of flights that fly through the kth entry point in the sector Statistics of corresponding flights Complete the normal distribution fitting of the statistical data and estimate the random variables of all n entry points in the sector The mean μ in the normal distribution probability density function in,k and standard deviation Step 2-4: Filter out the flight set that passes through the sector exit point Statistics of corresponding flights Complete the normal distribution fitting of statistical data and estimate the random variables The mean μ in the normal distribution probability density function out and standard deviation σ out .

4. The method according to claim 3, characterized in that: Step 3 includes: Step 3-1: According to the start and end time of the collaborative interval management conversion program CMCP and the sector object of the restriction constraint, filter out the following related flight sets: The flight set F affected by CMCP; The set of flights F that are affected by CMCP and fly through the kth entry point in the sector in,k ; The set of flights F that are affected by CMCP and fly through the sector exit point out ; The set of flights affected by CMCP and subject to ground delays F G ; The set of flights F that are affected by CMCP and controlled by interval management M ; Set F M The set of flights F that fly through the kth entry point in the sector M,k ; Step 3-2: Construct the collaborative interval management conversion program CMCP. The decision variable of the opportunity-constrained programming model is the delay time d assigned to flight f. f ; Step 3-3, construct the collaborative interval management conversion program CMCP. The optimization objectives of the opportunity-constrained programming model are efficiency objectives and fairness objectives, where the efficiency objective is to minimize the total flight delay cost, and the fairness objective is to minimize the maximum airline fairness ratio deviation; Step 3-4, construct the constraints of the collaborative interval management conversion program CMCP opportunity constraint programming model; Step 3-5: Based on the opportunity constraint equivalent transformation theory, the opportunity constraints in the collaborative interval management transformation program CMCP opportunity constraint programming model are equivalently transformed.

5. The method according to claim 4, characterized in that In step 3-3, the total flight delay cost is calculated as follows: Where C represents the total flight delay cost, and β represents the ratio of the cost per minute of air delay to the cost per minute of ground delay.

6. The method according to claim 5, characterized in that In step 3-3, the maximum airline fair ratio deviation is calculated as follows: Where E represents the maximum fair proportion deviation of airlines, A represents the set of airlines to which the affected flights belong, l represents the total number of airlines in set A, and ρ a represents the fair proportion of airline a, ρ a The calculation method is: Among them, C a represents the total flight delay cost of airline a, q a represents the total number of affected flights belonging to airline a, and q represents the total number of affected flights.

7. The method according to claim 6, characterized in that In step 3-4, the constraints of the collaborative interval management conversion program CMCP opportunity constrained programming model are: Among them, Pr{} represents the probability of satisfying the conditional constraints in the curly brackets {}, and f ′ indicates another flight different from flight f, represents the time interval requirement between two consecutive flights entering from the kth entry point of the sector, α in,k represents flight f and flight f ′ The confidence level of the constraint on the interval between the i-th entry point of the sector is maintained, s out Indicates the time interval requirement between two consecutive flights leaving from the sector exit point, α out represents flight f and flight f ′ The confidence level at which the sector exit point interval constraint holds, and N represents a set of natural numbers.

8. The method according to claim 7, characterized in that In step 3-5, the result of equivalent transformation of the opportunity constraints in the collaborative interval management transformation program CMCP opportunity constraint programming model is: in, and Δs out The calculation method of is as follows: in, It represents the additional interval time that needs to be maintained between two consecutive flights entering from the kth entry point of the sector, Δs out Indicates the additional interval time between two consecutive flights departing from the sector exit point. Represents the inverse of the standard normal distribution function.

9. The method according to claim 8, characterized in that Step 4 includes: Step 4-1, initialize the first generation particle swarm: set j = 1, set the total number of particles to n p , the maximum number of iterations is j max , randomly generate a particle swarm p j The initial positions and initial velocities of all particles in , where the subscript j represents the round of iterative calculation in the algorithm; Step 4-2, allocate flight delays: use the delay allocation operator as the particle swarm p j Each particle in the Corresponding to a set of feasible solutions, where the delay allocation operator refers to the position of particles one by one in order based on the greedy algorithm The earliest available time for each point in the flight allocation until the particle position All flights complete the transit time allocation and form a set of resource allocation schemes, that is, a set of feasible solutions; Step 4-3, calculate the objective function value: according to the resource allocation plan after the delay allocation operator acts, calculate the particle swarm p j The total flight delay cost and the maximum airline fairness ratio deviation are the two optimization target values ​​for each particle; Step 4-4, update the elite set: according to the objective function value of each particle, calculate the dominance level and crowding degree of each particle, and then based on the elite retention strategy, update the elite set from the particle group p j Select n e =ξ·n p Excellent particles as elite set Among them, ξ is the preset ratio threshold parameter; Step 4-5: Implement game operation: According to the elite retention strategy, from the particle swarm p j Select half of the particles as the parent particle group Then, from the parent particle swarm Take out the i-th particle to be updated from the elite set Randomly select two candidate elite particles from the , calculate the vector angles θ1 and θ2 between the two candidate elite particles and the current i-th particle to be updated, and use the candidate elite particle with the smaller angle as the guide particle of the current i-th particle to be updated. The position of the guide particle is recorded as Among them, the vector angle θ between the i-th particle to be updated and the z-th candidate elite particle is z The calculation formula is: in and They represent the efficiency target value and fairness target value of the i-th particle to be updated, respectively. and They represent the efficiency target value and fairness target value of the zth candidate elite particle respectively; Step 4-6, Update particle information: parent particle group The speed and position of each particle in the Update to obtain the descendant particle swarm Step 4-7, merge the parent particle swarm: and Merge to generate a new generation of particle swarm p j+1 , when the number of iterations meets the set maximum number of iterations j max Then, go to step 4-8, otherwise, let j=j+1 and go to step 4-2; Step 4-8: Output Pareto front solution set: for particle swarm p j+1 The delay allocation, objective function calculation and dominance ranking are redone to determine the Pareto front solution set of the current latest generation of particle swarm.

10. The method according to claim 9, characterized in that In steps 4-6, the parent particle swarm The speed and position of each particle in the To update, the formula is: in, represents the particle swarm p j The position of the i-th particle in represents the particle swarm p j+1 The position of the i-th particle in represents the particle swarm p j The velocity of the i-th particle in represents the particle swarm p j+1 is the velocity of the ith particle in , r1 and r2 are non-negative constants, η1 and η2 are random numbers in [0,1], The operation is used to subtract the positions of two particles, and the result is the flight speed of the generated particle; The operation is used to add the particle position and particle velocity, and the result is the new position of the added particle; The operation is used to merge two particle speeds, and the result is to generate a new particle speed; The operation is used to multiply the constant by the particle velocity. The result is the generation of a new particle velocity. The operation form is based on the constant λ C There are three different cases for different values ​​of : When the constant λ C = 0, the new particle velocity after multiplication is an empty set; When 0<λ C <1, the new particle speed after multiplication is the original particle speed The remainder of the commutators, where h represents the length of the multiplied velocity, Represents h and λ C The integer rounded up after multiplication; When C ≥1, C Split into two new constants λ C,1 and λ C,2 , satisfying λ C =λ C,1 +λ C,2 , and λ C,1 is a non-negative natural number, 0≤λ C,2 <1, so the new particle speed after multiplication is λ C,1 The speed of the original particles is connected end to end and then merged into λ C,2 The result after multiplying by the original particle velocity is

Images