一种基于分布估计算法的航班时隙分配多目标优化方法

摘要

本发明公开了一种基于分布估计算法的航班时隙分配多目标优化方法，本发明针对现有方法对时隙分配多目标优化算法研究不足，同时对于资源约束下航班关联性研究甚少等现状，提出基于分布估计算法的多目标航班时隙分配启发式优化方法。本发明考虑地面等待策略中时隙分配问题，建立包含延误成本和公平性指标在内的资源约束下多目标优化问题，在约束条件中还加入起飞机场容量，防止造成其他繁忙机场容流不平衡。基于分布估计算法，引入马尔科夫网络对问题决策变量进行建模，进行描述、挖掘多变量间的关联性，增加搜索优化效率；同时针对多目标优化过程提出带有指数化权重参数的适应度函数，在低计算复杂度和灵活参数配置下，提高帕累托解的分布性。

说明书

技术领域

本发明属于空中交通流量管理领域，尤其涉及一种基于分布估计算法的航班时隙分配多目标优化方法。

背景技术

在流量管理中，地面等待策略是一种主动实施、非常重要的流量管理方法。实质是将昂贵且高风险的空中等待转化为安全经济的地面等待。流量管理部门预先了解到航空器所经过的航路、扇区或到达的目的地机场即将容量下降导致空中交通流量拥堵，则要求航空器在起飞机场地面等待来调节空中交通网络的流量，实现容需平衡，从而达到优化航班延误成本、提高机场、空域利用率和降低飞行安全风险的目的。

在地面等待策略中，一项核心技术就是航班时隙分配。即机场由于容量受限产生多个时隙资源(可看做是进场时隙)，同时存在多架次进场航班，要求每个时隙只能至多分配给一个航班，并且每个航班必须配置一个比航班计划进场时间晚的时隙。一般来说，延误成本是必须要考虑的目标，除此以外，公平性指标也是需要考量的重要因素。因此，时隙分配问题往往是一种(时隙)资源约束下的多目标优化问题。

目前，时隙分配问题的研究一方面集中在问题建模上，包含确定性、不确定性、动态、多目标等问题；另一方面在优化方法上进行研究，例如整数规划、随机规划、鲁棒优化等。

针对时隙分配问题，现有方法对多目标优化算法研究不足，同时对于资源约束下航班关联性研究甚少。针对上述不足，本发明提出利用马尔科夫网络对航班关联性进行建模，通过基于概率统计的分布估计算法进行搜索优化，同时对多目标优化给出新的适应度函数，以此建立一种针对航班时隙分配的多目标优化方法。

发明内容

发明目的：在确保机场容量约束前提下，充分利用时隙，实现延误成本和公平性指标的多目标优化，并满足计算时间复杂度要求。在优化过程中，给出时隙分配解决方案的同时，还对航班之间关联性进行挖掘，为辅助决策提供数据支持。

技术方案：一种基于分布估计算法的航班时隙分配多目标优化方法，包括如下步骤：

步骤1，基于效率性和公平性建立时隙分配多目标模型；

步骤2，基于时隙分配多目标模型，建立对应马尔科夫网络；

步骤3，建立概率矩阵；

步骤4，采用新的适应度函数，计算并找出当代帕累托解；

步骤5，根据当代帕累托解，更新信任解集存档；

步骤6，学习马尔科夫网络，更新概率矩阵；

步骤7，根据概率矩阵和马尔科夫网络进行采样，生成新解；

步骤8，达到终止条件，获得最终结果。

本发明中，步骤1包括：

步骤1-1，建立如下多目标函数：

其中，N表示航班总数，M表示时隙总数，c_ij表示航班i被分配给时隙j时的延误成本，x_ij表示决策变量，F(x_ij)表示航班i被分配给时隙j时的公平性指标；公式(1)表示最小化延误成本和公平性差异指标的多目标函数。

步骤1-2，设定约束条件：

c_ij＝ω_i×(t_j-SLDT_i)>

t_j≥SLDT_i,当x_ij＝1时>

其中，ω_i表示航班i在延误时间方面的重要性权重，取值范围0.0～1.0；Q表示航空公司总数，N_q表示航空公司q的航班总数，t_j表示时隙j的开始时间；x_ij为决策变量，代表航班i是否被分配给时隙j，仅能取值0或者1；SLDT_i表示航班i的时刻表计划着陆时间，c_i^*表示航班i最优情况的延误，d(x_ij)表示航班i被分配时隙j时和最优情况的差值，nd(x_ij)表示正则化后的d(x_ij)的差值，f_q(x_ij)表示航空公司q的平均正则化后的差值；公式(2)为问题的决策变量，公式(3)和(4)分别表示航班分配的约束条件和时隙分配的约束条件，公式(5)用于计算航班延误成本，公式(6)为航班时刻表约束，公式(7)和(8)计算实际时隙分配与最优分配下的差值并进行正则化，公式(9)和(10)计算所有参与时隙分配公司的方差，作为航空公司公平性指标。

本发明中，步骤2包括：建立对应时隙分配多目标模型的马尔科夫网络，在建立的马尔科夫网络中，每个节点代表一个决策变量，即航班的时隙分配选择；节点X_i表示航班i的时隙选择，其取值范围D(X_i)＝{s₁,s₂…,s_M}，其中M为可用时隙总数，s_M表示第M个可用的时隙；在建立的马尔科夫网络中，两个节点之间的连线表示两个航班之间存在较强耦合性和航班时隙争夺。

本发明中，步骤3包括：建立如下概率矩阵：

P(t)表示迭代计算过程中第t代的概率矩阵；

P_NM表示矩阵第N行第M列的值；

矩阵的行i表示航班i分配到每个时隙j上的概率，合计为100％；

矩阵的列j表示针对时隙j，每个航班i被分配到该时隙的概率，合计为100％；

初始化矩阵时，采用等概率进行赋值，针对时隙分配具体问题，存在公式(6)中的时间约束，制定如下规则：对每一行i，在M个时隙中，将早于SLDT_i的时隙j的概率p初始值设置为0，将晚于SLDT_i的时隙j的概率p初始值设置为1/z，其中z为晚于SLDT_i的时隙数量。

本发明中，步骤4包括：

步骤4-1，帕累托解归纳为两部分：边界部分和中间部分，边界部分是单目标优化部分，中间部分是多目标优化部分，针对边界部分，采用VEGA算法中的边界适应度函数来计算：

EdgeFitness_f(X)＝obj_value_f(X)>

其中X表示任意解，EdgeFitness_f(X)表示解X在目标f的适应度函数，obj_value_f(X)表示解X在目标f上函数值；

针对中间部分，采用如下中间适应度函数进行计算：

其中X表示任意解，CentralFitness(X)表示解X在中间部分的适应度函数，p(X)和q(X)分别表示被X支配(dominate)的数量以及支配X的数量；

所有帕累托解的中间部分适应度函数值都是大于等于1，同时，所有非帕累托解都小于1大于0；

根据公式(13)和(14)，分别计算出边界适应度和中间适应度；

步骤4-2，采用如下带有指数化权重参数的复合适应度函数来结合边界适应度和中间适应度函数：

其中D-Fitness(X)是解X考虑复合中间部分适应度和边界部分适应度的综合适应度，m是目标总数，N_POP是人口总数，Ranking_EdgeFitnessi(X)是基于第i个目标的解X的排名，Ranking_{CentralFitness}(X)是针对X在中间适应度上的排名，ω_i和ω_m+1是指数化权重参数(例如双目标优化问题，m＝2，所以指数化权重参数有3个：ω₁，ω₂，ω₃)，，用来区分不同排名的得分，同时能够实现对一个目标的偏好。

本发明中，步骤5包括：在数据库中建立一个存档用于存放信任解，作为信任解集存档，容量和每一次迭代中的人口数量N_POP相同，或者小于N_POP；例如取值N_POP×1/2。

通过公式(15)计算获得本代的所有解的适应度函数值，并根据适应度函数值由高到低选择α％的解存入信任解集存档，从上一代信任解集存档中根据适应度函数值由高到低选择1-α％的信任解存入本代的信任解集存档。

本发明中，步骤6包括：

步骤6-1，学习马尔科夫网络：

通过公式(16)，计算变量X_i和Y_j的互信息量MI(X_i；Y_j)：

其中，X_i和Y_j表示两个随机变量，时隙分配问题中，变量具体代表的是一个航班的时隙分配；p(x_i|D)和p(y_j|D)分别表示基于可信解集合D的变量X_i＝x_i的概率和Y_j＝y_j的概率，p(x_i,y_j|D)是X_i＝x_i和Y_j＝y_j的联合概率；

如果两个变量的互信息量大于给定阈值threshold(例如可以取值为0.5)，则判定两变量之间存在较强关联，在马尔科夫网络中建立连接线并互称为邻居；阈值threshold由用户设定，或者在计算过程中，根据已知的互信息量来动态更新阈值，如公式(17)所示：

其中，MI表示两个变量X_i和Y_j的互信息量，NDV表示变量的总量，参数β用于控制结构复杂度(通常可设置为0.5-1.5之间)；

步骤6-2，更新概率矩阵：

变量X_i在其邻居N_i的概率矩阵，如公式(18)所示：

其中P(X_i|N_i)表示变量X_i在所有其邻居N_i的条件概率矩阵，{x_i,1，…，x_i,j，…，x_i,Yi}分别表示变量X_i能够被赋予的值，在时隙分配中即对应每一个待分配的时隙编号，x_i,Yi表示变量X_i第Y_i个能够被赋予的值，Y_i表示变X_i取值的最大编号，x_i,Yi即为X_i能够被取值的最大值，时隙分配问题中，即为最晚的一个时隙；p(x_i,Valuej|n_i,Neighbork)表示变量X_i取值Value_j时，在其所有邻居联合取值场景为Neighbor_k下的条件概率；Z_i表示邻居集合N_i的联合取值场景的最大编号；

对于公式(12)和(18)中的概率，首先计算当前代的边界概率P_t^*(X＝x)：

其中，N(X＝x)表示在可信任解中变量X取值x的数量，N(D)表示可信任解的总量；通过公式(20)更新概率模型：

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

其中P_t(X＝x)和P_t-1(X＝x)分别表示在迭代计算过程中第t代时变量X取x值的概率和第t-1代时变量X取x值的概率，λ是从当前代的学习率，如果λ＝1，则表示概率模型将完全从当前代信任解集存档中学习；

步骤6-3，在通过公式(20)更新概率模型之后，通过变异概率(例如5％)执行公式(21)的变异操作：

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

其中θ表示变异值(例如取0.2)，ε是一个很小的正数(例如取0.01)。

本发明中，步骤7包括：

步骤7-1，随机生成一个可行解s＝(s₁,s₂,…,s_n)，s_n为最后一个变量X_n的取值，时隙分配问题中s₁,s₂,…,s_n代表每个航班的时隙选择，s_n表示第n架航班的时隙选择，可行解s生成方法如下：将所有航班根据航班时刻表排列，从第一个航班x₁开始，从所有剩余时隙中，随机选择一个时隙s₁分配给x₁，要求s₁不早于x₁的计划起飞时间，并且时隙s₁没有被其他航班选择，删除时隙s₁，移动到下一个航班x₂，用同样的方法随机选择，直至所有航班选择结束；

步骤7-2，将解s＝(s₁,s₂,…,s_n)顺序随机打乱，得到新的顺序s₁^*,s₂^*,…,s_n^*，其中s_n^*表示新的航班排序中最后一个航班的时隙选择；需要说明的是：此时的航班以及时隙对应关系并没有改变，只是排列的次序发生变化；

步骤7-3，从s₁^*开始，设定s₁^*对应航班h，根据建立的马尔科夫网络找出航班h变量的所有邻居，并根据此时解s＝(s₁,s₂,…,s_n)中邻居的取值以及对应的条件概率p(x_h|N_h)来决定航班h选择各个时隙的概率，再按照概率分布，采用轮盘赌选择算法(Roulette>1^*的新值；然后移动到s₂^*，采用同样方法再计算s₂^*的新值，直到完成所有变量，则完成了一个解的采样；

步骤7-4，重复步骤7-1至步骤7-3的方法，完成一个群体N_POP个解的采样，最终形成当前代的解集。

有益效果：与现有方法相比，本发明的优点在于：

采用分布估计算法，相比遗传算法等启发式算法，全局搜索能力和收敛速度有所提高；

采用马尔科夫网络结构，更好描述和发掘变量关系，提供更多辅助决策能力；

采用基于指数化权重参数适应度函数，实现低计算复杂度和搜索灵活性。

本发明方法在确保机场容量约束前提下，充分利用时隙，实现延误成本和公平性指标的多目标优化，并满足计算时间复杂度要求。在优化过程中，给出时隙分配解决方案的同时，还对航班之间关联性进行挖掘，为辅助决策提供数据支持。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述或其他方面的优点将会变得更加清楚。

图1是基于马尔科夫网络的多目标分布估计算法流程。

图2是分布估计算法流程。

图3是马尔科夫网络结构示例。

图4是包含边界部分和中间部分的帕累托解。

图5是针对中间部分的适应度函数。

图6a是最小化问题的示例。

图6b是带有指数化权重参数的函数。

图7是更新存档候选解集。

图8是实验验证结果(部分)。

图9是基于分布估计算法的航班时隙分配多目标优化方法的流程图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

如图1所示，是基于马尔科夫网络的多目标分布估计算法流程示意图。算法是基于种群进化理论(Population-to-Population)，每一代产生新的候选解，并联合上一代选择出的信任解，通过适应度函数进行多目标帕累托选择，以此更新信任解。然后基于新得到的信任解，通过分布估计算法中概率模型估计和采样，得到下一代的候选解，一直循环直到终止条件出现。其中需要说明的是，马尔科夫网络引入是和分布估计算法做结合，增强概率模型的估计和采样过程。

图2是采用的分布估计算法流程示意图。其中重点是采样生成新解，和估算概率模型两个方面。

如图9所示，本发明具体过程如下：

步骤1：基于效率性和公平性建立时隙分配多目标模型。

(1)参数

参数名称参数说明

N 航班总数

M 时隙总数

Q 航空公司总数

i 航班i，i＝1,…,N

j 时隙j，j＝1,…,M

q 航空公司q，q＝1,…,Q

N_q航空公司q的航班总数

c_ij航班i被分配给时隙j时的延误成本

x_ij航班i是否被分配给时隙j

F(x_ij)>

t_j>

SLDT_i航班i的时刻表计划着陆时间

c_i^*航班i最优情况的延误

d(x_ij)>

nd(x_ij)正则化后的d(x_ij)差值

f_q(x_ij)航空公司q的平均正则化后的差值

ω_i航班i在延误时间方面的重要性权重

(2)多目标函数

(3)约束条件

c_ij＝ω_i×(t_j-SLDT_i)>

t_j≥SLDT_i,当x_ij＝1时>

其中，公式(1)表示最小化延误成本和公平性差异指标的多目标函数。公式(2)中x_ij为决策变量，代表是否航班i被分配给时隙j，仅能取值0或者1。公式(3)和(4)确保对于任何一个待分配的航班，有且仅有一个时隙分配给它，而对于任何一个待分配的时隙，最多仅可以被分配给一个航班。公式(5)表示根据分配的时隙的开始时间和时刻表时间的差值，计算延误成本。公式(6)保证对于任何一个航班，分配的时隙必须要等于或者晚于时刻表时间。公式(7)和(8)计算分配方案下延误成本与最优情况的延误成本的差值并正则化。公式(9)和(10)计算所有参与时隙分配公司的方差，作为航空公司公平性指标。

步骤2：基于时隙分配问题，建立对应马尔科夫网络。

马尔科夫网络为一网状无向图(或称为双向图)，由结构G和参数Ψ组成。结构G中包含节点Node和连接边E，因此一个马尔科夫网络可以表示为：MN＝{G,Ψ}＝{Node,E,Ψ}。如图3所示，马尔科夫网络中含有5个变量(x₁-x₅)，变量间的连线代表这两个变量间存在某种强关联，互相之间被称为邻居。例如变量X₁含有2个邻居X₂和X₃；X₂含有3个邻居，X₁，X₃和X₄。

变量X_i(D(X_i)＝{x_i¹,…,x_iⁿⁱ})的取值，(x_iⁿⁱ表示X_i能取的最大值)是由其所有邻居的条件概率计算得来的，如公式(11)所示。

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

其中，N_i是变量X_i的邻居的集合，p(x_i|x-{x_i})表示变量X_i在其他所有变量联合取值下的条件概率，p(x_i|N_i)表示变量X_i在它所有邻居联合取值下的条件概率。

针对时隙分配优化问题，此步骤负责建立对应马尔科夫网络，包含网络中各节点的定义和取值范围。

在马尔科夫网络中，每个节点代表一个决策变量，具体来说就是航班的时隙分配选择。节点X_i表示航班i的时隙选择，其取值范围D(X_i)＝{s₁,s₂…,s_M}，即所有的待分配时隙。两个节点之间的连线表示两个航班之间存在较强耦合性和资源(航班时隙)争夺，作用是一方面为进化过程提供邻居信息，计算概率分布；另一方面为后续例如时隙交换等操作提供一定依据。

步骤3：建立概率矩阵。

此时概率矩阵是基于上步骤中的马尔科夫网络的。概率模型如公式(12)所示，P(t)表示迭代计算过程中第t代的概率矩阵名称，为一N×M的矩阵。行i表示航班i分配到每个时隙j上的概率，合计为100％。以列j来看，表示针对时隙j，每个航班i被分配到该时隙的概率，合计也为100％。

初始化时，一般采用等概率进行赋值。针对时隙分配具体问题，存在公式(6)中的时间约束，同时为了加快采样和搜索的效率，制定如下规则：对每一行i来说，在M个时隙中，将早于SLDT_i的时隙j的概率p初始值设置为0，将晚于SLDT_i的时隙j的概率p初始值设置为1/Q，其中Q为晚于SLDT_i的时隙数量。

步骤4：采用新的适应度函数，计算并找出当代帕累托解。

数学模型中包含两个目标需要同时进行优化，设计出一套高效率的适应度函数来完成解的搜索是十分重要的。在不同的多目标优化算法中，适应度函数对于搜索效率和效果起到相当大的作用。本发明提出一种带有指数化权重参数的适应度函数，能够提高搜索效率的同时，带来更高的灵活性。

如图4所示，以最小化多目标问题为例，帕累托解归纳为两部分：边界部分(单目标优化部分)和中间部分(多目标优化部分)，针对边界部分，采用VEGA算法中的边界适应度函数来计算：

EdgeFitness_f(X)＝obj_value_f(X)>

其中X表示任意解，EdgeFitness_f(X)表示解X在目标f的适应度函数，obj_value_f(X)表示解X在目标f上函数值。

针对中间部分，采用如下中间适应度函数进行计算：

其中X表示任意解，CentralFitness(X)表示解X在中间部分的适应度函数，p(X)和q(X)分别表示被X支配(dominate)的数量以及支配X的数量；

从图5中可以看出，公式(14)给出的适应度函数能够很好的区分出是否属于帕累托集合。所有支配解(帕累托解)的CentralFitness值都是大于等于1，而且支配其他解越多，数值越高。同时，所有非支配解都小于1大于0。因此，数值1可以清楚地分辨一个解是否为支配解。

根据公式(13)和(14)，可以分别计算出边界适应度和中间适应度，而且计算复杂度非常低。

下一步为了能够获得最终解，需要综合考虑上述两个适应度函数。有两种方案，一种是采用多岛(Multi-island)优化，即部分种群依照边界适应度函数优化，另外一部分种群依照中间适应度优化。这种方法执行效率高，但带来的问题是很难最终形成统一解以及帕累托解的最优性不高。

算法中采用带有指数化权重参数的函数来结合边界适应度和中间适应度函数，如公式(15)所示：

其中D-Fitness(X)是解X考虑复合中间部分适应度和边界部分适应度的综合适应度，m是目标总数(双目标问题，m＝2)，N_POP是人口总数，Ranking_EdgeFitnessi(X)是基于第i个目标的X的排名(EdgeFitness越高，排名越靠前)，Ranking_{CentralFitness}(X)是针对X在中间适应度上的排名(CentralFitness越高，排名越靠前)，ω_i和ω_m+1是指数化权重参数(例如双目标优化问题，m＝2，所以指数化权重参数有3个：ω₁，ω₂，ω₃)，，用来区分不同排名的得分，同时可以实现对某一个目标的偏好。

以图6a和图6b为例说明指数化权重参数。设定问题为最小化问题，在图6a中，考量黑色解和白色解，黑色解为帕累托解，而白色解不是。黑色解的特征是在某一个维度上非常优秀(边界，f₂)，白色解虽然在边界f₂这个维度上比黑色解略微差一些，但是其他维度上都要比黑色解好。帕累托解只需在某一个维度中非常“杰出”的，而不是在每个维度都“不错”的。因为为了能够明显区分并且快速找寻这样的帕累托解，图6b是带有指数化权重参数的函数，即和常见的线性关系相比，在某一维度一个解的排名越靠前，它的得分(适应度函数值)比其他解高得越多。

同时，采用指数化权重参数还有一个非常重要的作用：实现搜索的灵活性。类似于权重的概念，如果某一个维度出现在帕累托解中的数量不足，则需要加大对该维度的发掘。通过增大对应的指数化权重参数，对该维度从最终适应度函数值的计算上侧重，“鼓励”该维度能够出现更多的优秀解出现，也提高了帕累托解的分布性。

步骤5：根据当代帕累托解，更新信任解集存档。

如图7所示，基于本代的候选解(Candidate Solutions)和上代的信任解(Promising Data)，通过上步骤中的复合适应度函数，两部分中都有一部分解成为帕累托解。信任解集存档(Archive)用于存放更新后的信任解。本算法中采用固定大小的方式，其的优势是有更多的解用于之后的概率估算，防止出现过早成熟、陷入局部最优的问题。对于选取哪些候选解作为非帕累托解进入存档中，为了保证搜索的灵活性、解的多样性和存档的稳定性，通过公式(15)计算获得本代的所有解的适应度函数值，并根据适应度函数值由高到低选择α％的解存入信任解集存档，从上一代信任解集存档中根据适应度函数值由高到低选择(1-α％)的信任解存入本代的信任解集存档。

步骤6：学习马尔科夫网络，更新概率矩阵。

本算法和传统的分布估计算法不同的是引入了马尔科夫网络，因此在分布估计时，除了概率矩阵需要更新，马尔科夫网络的结构也需要重新学习。下面从两个方面分别说明：

(1)学习马尔科夫网络

马尔科夫网络中两个节点之间的联系表示两个变量之间具有较强的关联性(非因果关系)，并且变量的赋值也是基于和它存在关联的“邻居”的，因此马尔科夫网络的学习过程对于本算法是很重要的一个环节。

一种建立网络结构的方法是由问题专家负责，但是对于大部分问题很难找到专家能够对变量之间的关联性给出可信意见。本算法采用一种条件独立测试(ConditionalIndependence Test)——互信息(Mutual Information，MI)来估算马尔科夫网络结构。互信息的优势在于使用简单并且计算快速。

通过公式(16)，计算变量X_i和Y_j的互信息量MI(X_i；Y_j)：

其中，X_i和Y_j表示两个随机变量(时隙分配问题中，变量具体代表的是某个航班的时隙分配)，p(x_i|D)和p(y_j|D)是基于可信解集合D的变量X_i＝x_i的概率和Y_j＝y_j的概率，p(x_i,y_j|D)是X_i＝x_i和Y_j＝y_j的联合概率；

如果两个变量的互信息量大于给定阈值threshold(例如取值0.5)，则判定两变量之间存在较强关联，在马尔科夫网络中建立连接线并互称为邻居；阈值threshold由用户设定，或者在计算过程中，根据已知的互信息量来动态更新阈值，如公式(17)所示：

其中，X_i和Y_j表示两个随机变量，MI表示两个变量X_i和Y_j的互信息量，NDV表示变量的总量，参数β用于控制结构复杂度(通常可设置为0.5-1.5之间)。如果β被设置为一较高值，则阈值较高，意味着马尔科夫网络中拥有较少的连线和较短的计算时间，但会损失一部分关联信息。如果β被设置为一较低值，则阈值变低，意味着马尔科夫网络中拥有较多的连线和更多的关联信息被保留，有助于提高搜索的质量，但是需要较长的计算时间。因此，参数β可以用来控制搜索质量和效率。

结构学习往往需要很长时间，因此为了搜索效率，并不需要每一代都进行学习。例如可以设置为每10代进行结构学习一次(概率矩阵每一代都需要更新)。

(2)更新概率矩阵

在公式(12)中，概率矩阵表示为每个变量的独立概率。同时由于在马尔科夫网络中存在关联性，还需要建立和更新基于邻居的条件概率矩阵。变量X_i在其邻居N_i的概率矩阵，如公式(18)所示，

其中P(X_i|N_i)表示变量X_i在所有其邻居N_i的条件概率矩阵，{x_i,1，…，x_i,j，…，x_i,Yi}分别表示变量X_i能够被赋予的值，在时隙分配中即对应每一个待分配的时隙编号，Y_i表示变X_i取值的最大编号，x_i,Yi表示X_i能够被取值的最大值(时隙分配问题中，即为最晚的一个时隙)；表示变量X_i的某取值Value_j时，在其所有邻居联合取值场景为Neighbor_k下的条件概率，Z_i表示邻居集合N_i的联合取值场景的最大编号。

对于公式(12)和(18)中的概率，都是根据计算获得的可信任解或者可信任解集存档中的内容进行更新。由于是迭代搜索的，首先计算当前代的边界概率：

其中，N(X＝x)表示在可信任解中变量X取值x的数量，N(D)表示可信任解的总量，1/|X|的设置是为了防止出现0的可能，对概率计算不产生实质影响。

计算完当前代的概率，需要更新概率模型，为了保证模型的稳定性，通过公式(20)进行计算：

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

其中P_t(X＝x)和P_t-1(X＝x)分别表示在迭代计算过程中第t代时变量X取x值的概率和第t-1代时变量X取x值的概率，λ是从当前代的学习率，如果λ＝1，则表示概率模型将完全从当前代得到的可信解中学习。考虑到算法中已经设计了归档机制，可以将λ的值设置大于0.5。

为了能够进一步增加解的多样性，避免陷入局部最优，类似GA的变异操作也被加入到概率更新的计算过程中。在公式(20)更新概率之后，通过变异概率(例如5％)执行公式(21)的变异操作：

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

其中θ表示变异值(例如取0.2)，ε是一个很小的正数(例如取0.01)，为了保证所有的概率都要小于100％。

步骤7：根据概率矩阵和马尔科夫网络进行采样，生成新解。

通过上个步骤已经获得概率模型和网络模型，下一步需要生成下一代的解。马尔科夫网络和贝叶斯网络(Bayesian network，BN)有很大的不同，并不存在因果关系，因此在采样生成新解过程中，并没有明确的顺序。为了得到新解，算法采用一种基于马尔科夫链的蒙特卡洛方法——吉布斯采样(Gibbs Sampler)。具体过程如下：

(1)随机生成一个可行解s＝(s₁,s₂,…,s_n)，s_n为最后一个变量X_n的取值(时隙分配问题中代表每个航班的时隙选择，s_n表示第n架航班的时隙选择)。生成方法如下：将所有航班根据航班时刻表排列，从第一个航班x₁开始，从所有剩余时隙中，随机选择一个时隙s₁(要求s₁不早于x₁的计划起飞时间，并且时隙s₁没有被其他航班选择)分配给x₁，删除时隙s₁，移动到下一个航班x₂，用同样的方法随机选择，直至所有航班选择结束。

(2)将解s＝(s₁,s₂,…,s_n)顺序随机打乱，得到新的顺序s₁^*,s₂^*,…,s_n^*，其中s_n^*表示新的航班排序中最后一个航班的时隙选择。(需要说明的是：此时的航班以及时隙对应关系并没有改变，只是排列的次序发生变化)。

(3)从s₁^*(假设对应航班h)开始，根据建立的马尔科夫网络找出航班h变量的所有邻居，并根据此时解s＝(s₁,s₂,…,s_n)中邻居的取值以及对应的条件概率p(x_h|N_h)来决定航班h选择各个时隙的概率，再按照概率分布，采用轮盘赌选择算法(Roulette>1^*的新值；然后移动到s₂^*，采用同样方法再计算s₂^*的新值，直到完成所有变量，则完成了一个解的采样。

重复上述步骤(1)至(3)，完成一个群体N_POP个解的采样，最终形成当前代的解集。

步骤8：达到终止条件，获得最终结果。

重复以上步骤，直至满足终止条件。终止条件可以设置为迭代次数，例如1000代。当达到终止条件时，此时获得的帕累托解即为最终的结果。

为了验证本发明方法中模型和算法的正确性，基于华北地区流量管理环境，选取了18：00-24：00共计6个小时的模拟航班数据用作验证，并假设22：00-23：00之间北京首都机场(ZBAA)出现容量下降和超容情况。

航班时隙分配的结果如图8所示(部分结果)。图中结果给出了受影响航班的预计起飞时间和计算起飞时间，并给出了地面延误时间的建议。从结果可以比较直观地看出，一方面有较多航班的公司的部分航班是不需要延误的，说明符合多目标优化中的效率性(延误时间)指标，另一方面也可以发现平均延误时间比较类似，并没有出现某一家航空公司的航班出现过多延误，也证明了公平性指标得到保证。通过实验仿真验证，说明本方法可以有效地解决时隙分配问题，并按照问题建模的多目标优化计算获得满意的结果。

本发明提供了一种基于分布估计算法的航班时隙分配多目标优化方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。