第五章-分布估计算法研究

第五章分布估计算法研究遗传算法是广泛应用的智能优化算法，通过选择、交叉和变异三种基本操作实现全局最优解的搜索。遗传算法获得巨大成功的理论基础是建筑块理论[1]，该理论指出遗传算法具有通过组合个体的建筑块获得最优解的能力。建筑块理论能够有效使用的重要条件是在数据编码中能够明确表示建筑块的位置[2,3]，否则遗传算法中的交叉和变异操作就有可能破坏建筑块，特别是那些较长的或距离较远的块，从而导致算法的收敛速度下降或找不到最优解[4]。分布估计算法（EstimationofDistributionAlgorithm,EDA）源于遗传算法，是为解决遗传算法中存在的因交叉、变异操作而导致的建筑块破坏问题而提出的[5,6]。在分布估计算法中没有交叉和变异操作，取而代之的是对于选择出来的优势群体的概率分布模型进行估计，并根据估计的模型进行采样。由于分布估计算法具有更强的理论基础，因此受到广大学者的重视，成为当前进化计算领域的研究热点[7~10]。5.1引言与遗传算法类似，分布估计算法也是基于群体的一种进化算法，算法由一个初始群体开始，在每代循环中，首先根据每个个体的适应值选择出具有代表性的适应值较好的一些个体组成优势群体。然后根据优势群体估计其中的个体所服从的概率分布模型，比如联合正态分布。第三步便是根据估计的概率分布模型采样产生一些新个体，并用来替换原有群体中的部分个体，从而组成新一代的群体。当满足循环终止条件时，算法结束，当前群体中的适应值最好的个体就是算法优化的结果。与遗传算法相比，分布估计算法中没有了交叉和变异两个操作，而是对选择群体的概率分布模型进行估计和采样。由于建筑块中的变量相关性要比其他变量之间的相关性强。分布估计算法更注重变量的相关结构，对建筑块的分布情况进行整体估计，这样不仅可以避免交叉变异所带来的建筑块破坏问题，反而有助于建筑块的增长。分布估计算法的基本步骤如下：Step1：初始化群体。在搜索空间内按均匀分布随机产生PS个点，组成初始群体，记进化代数g为1。Step2：计算适应值。根据适应值评价函数计算第g代群体中的各点的适应值。Step3：选择优势群体。根据适应值，运用一定的选择策略（如截断选择，轮赌选择等）选出适应值较好的s个个体组成优势群体。Step4：估计优势群体的概率分布模型。一般假设优势群体服从某一类概率分布模型，然后以优势群体为样本，对具体的概率分布F进行估计。Step5：根据估计的概率模型进行采样，产生一些新个体。在搜索空间内按概率分布F随机产生l个点，作为新生代中的部分个体。Step6：更新群体。由第g代群体中m个适应值较好的个体，新生的l个个体，以及随机产生（或者以较好个体为初始点搜索等方式产生）的PS-m-l个个体组成新一代群体，并将进化代数g←g+1。Step7：若满足某种停止准则，则算法结束，g代群体中的最好个体就是优化的结果；否则算法转到Step2继续执行。根据变量相关性的情况，一般认为分布估计算法可分为以下三类：各变量相互独立的分布估计算法，变量两两相关的分布估计算法和多变量相关的分布估计算法。前两类分布估计算法在运算复杂度方面优于第三类算法，但由于其不能充分考虑变量的相关性，因此在复杂问题的优化效果方面要次于第三类算法。下面对这三类算法中的代表性算法进行简单的介绍和分析。5.1.1变量相互独立的分布估计算法PBIL（Population-BasedIncrementalLearning）[11]和UMDA（UnivariateMarginalDistributionAlgorithm）[12]都是典型的变量相互独立的分布估计算法，它们仅考虑离散变量。对于选择出来的优势群体，根据每个变量的取值情况估计其一维边缘分布，由于在这类算法中所考虑的变量是相互独立的，因此反映群体分布情况的联合分布为各边缘分布之积。两种方法的不同之处在于它们更新概率模型的方式不同。PBIL采用线性调整概率模型的方式，在原有概率模型的基础上，根据优势群体的分布情况，各边缘分布按一定的学习速率进行线性调整，即pk(X)=∏ipk(xi)=∏i[α×pk-1(xi)+(1-α)×ni/N]。当α=0时，PBIL等价于UMDA。UMDA则根据优势群体的分布情况直接确定概率模型，即pk(X)=∏ip(xi|Dk-1s)。cGA（compactGeneticAlgorithm）[13]也是变量独立的分布估计算法，它需要较小空间，因为其群体规模仅为2，在每代运行时比较两个个体的适应值，并将较好的个体作为概率模型。在连续域内的变量独立的分布估计算法有PBILc[14,15]和UMDAcG[16],它们分别是PBIL和UMDA在连续域内的扩展，边缘分布使用的是一元正态分布。这类算法在解决优化变量之间独立的问题时表现非常好，并能有效解决遗传算法带来的建筑块破坏问题，但是在优化变量之间存在相关性的问题时则表现的非常差，其原因是算法中没有考虑到优化变量的相关性。5.1.2变量两两相关的分布估计算法最初考虑变量相关性的分布估计算法将变量关系简化为两两相关，典型的算法有MIMIC（MutualInformationMaximizationforInputClustering）[17]、COMIT（CombiningOptimizerwithMutualInformationTree）[18]和BMDA（BivariateMarginalDistributionAlgorithm）[19]。它们分别采用链状结构、树结构和森林来表示两两相关的变量关系。在MIMIC中概率模型表示为pkπ(X)=pk(xi1|xi2)…pk(xin-1|xin)pk(xin)，为了找到变量(x1,x2,…,xn)的最佳排列π=(i1,i2,…,in)，MIMIC算法使用贪婪算法最小化两个概率分布之间的K-L距离。COMIT采用树结构表示变量的两两相关性，通过反复最小化相互信息指标来向树结构中不断加入最大的枝。BMDA使用Pearson的χ2测试决定变量连接到哪棵树中以及树中变量的依赖关系。连续域的此类算法如MIMIC在连续域中的扩展MIMICcG[16]，其中使用了二元正态分布表示边缘分布。为找到最优排列，MIMICcG最小化的条件熵为。两两相关的分布估计算法考虑了优化变量的部分关系并且其结构学习也比较容易，在解决一些变量不独立的优化问题时表现也较好，比如算法COMIT和BMDA能够有效的优化二维spin-glass问题。但对于变量关系更为复杂的优化问题，这类算法则表现得较差。5.1.3多变量相关的分布估计算法多变量相关的分布估计算法能够全面考虑优化问题的结构，但其模型学习的复杂度也相对较高。根据表示变量关系的网络结构可以分为有向网络的分布估计算法、基于无向网络的分布估计算法和基于一般网络结构的分布估计算法。Bayesian网络是有向无环图，典型的基于Bayesian网络的分布估计算法有BOA（BayesianOptimizationAlgorithm）[20,21]、EBNA（EstimationofBayesianNetworkAlgorithm）[22,23]和EGNA（EstimationofGaussianNetworkAlgorithm）[24]。此类算法中使用Bayesian网络表示变量关系，根据选择出的群体估计网络的结构和参数，并根据估计的Bayesian网络模型进行采样生成下一代群体。BOA使用BD量（Bayesian-Dirichletmetric）和贪婪算法学习网络结构。EBNA使用BIC量（BayesianInformationCriteriametric）来衡量网络结构。此类算法在优化复杂的变量相关问题时效果很好，但是Bayesian网络的学习本身是NP难问题，因此算法的时间复杂度较高，多数时间花费在模型的估计上。EGNA是连续域分布估计算法，使用高斯网络表示变量的关系，该网络由贝叶斯得数、惩罚极大似然得数和边排除测试决定，变量的概率密度由高斯分布表示。Markov网络是无向图，基于Markov网络的分布估计算法有FDA（FactorizedDistributionAlgorithm）[25]、DEUM（DistributionEstimationUsingMarkovnetworkalgorithm）[26-28]、MN-FDA（MarkovNetworkFactorizedDistributionAlgorithm）[29]、MN-EDA（MarkovNetworkEstimationofDistributionAlgorithm）[30]、MOA（MarkovianitybasedOptimizationAlgorithm）[31]和MARLEDA（MarkovianLearningEstimationofDistributionAlgorithm）[32]。FDA用固定的概率图模型表示变量关系，但需要有关于变量关系的专家知识来完成概率分布函数的分解。DEUM根据无向图中的集群建立适应函数的模型，并将联合分布分解成Gibbs分布，模型参数通过解群体进行估计，采样时使用了Markov链蒙特卡洛方法（包括Gibbs采样和Metropolis采样）。MN-FDA用无向图近似表示联合分布，并从中构造出连接图，然后根据连接图采样。MN-EDA用联合分布的Kukichi近似来构建概率模型并用Gibbs采样产生新群体。MOA和MARLEDA是依赖于局部Markov性的分布估计算法，它们通过无向网络估计条件概率并由此采样。MARLEDA根据χ2测试估计结构，MOA采用基于交互信息的方式采样并利用退火温度表来控制探索与开采。这类算法同样适用于优化变量强耦合的问题，但其中的FDA对于黑箱问题或形式复杂的优化问题则无能为力。ECGA（ExtendedCompactGeneticAlgorithm）[33]是cGA的扩展，该算法将变量分成若干个相互独立的组，在每一组中变量之间有相关性，而与其它组中的任何变量都不相关。所有变量的联合分布为各组变量的分布函数之积。EMNA（EstimationofMultivariateNormalAlgorithm）[34]用所有变量的联合正态分布表示解分布模型，根据选择的群体估计均值向量和协方差矩阵作为模型参数，并从联合正态分布中采样得到新一代群体。钟伟才等在文献[35]中提出一种建立在一般Gauss网络上的算法，并采用Cholesky分解的方式产生新个体。这类算法全面考虑了变量之间的相关性，并且省去了复杂网络结构的学习，但是估计协方差矩阵的运算量是优化空间维数的平方，其时间复杂度仍然较高。目前已有的分布估计算法中要么不能充分考虑变量之间的相关性，对变量强相关的问题优化效果较差，要么模型估计的时间开销很大。最近的研究表明分布估计算法中90%的运算时间花费在模型估计阶段[36]，为减少算法的时间开销，学者们采用了不同的方法对其进行改进，如并行方式[37,38]、与其他算法融合的方式[39-41]、混合局部搜索的方式[42]，使用迭代的[43]或偶然的模型建立方式[44]以及初始信息优化的方式[45,46]等。如果基础算法中能够用较小的时间开销实现模型估计是最好的。所有变量的联合分布无疑是反映变量关系的准确且简捷的方式。5.2copula分布估计算法基于copula理论的分布估计算法是近几年来提出的一种新型的分布估计算法，据我们所有的中英文资料显示，我们在2009年IEEECEC会议中发表的论文[47]是这方面最早的论文。后来也有文献显示在2007年就有西班牙文的论文将高斯copula函数用于分布估计算法的研究中[48]。美国密苏里州分布估计算法实验室对分布估计算法的综述中特别指出，基于copula理论的分布估计算法是一种利用连续型概率模型对群体结构进行估计的解决连续型优化问题的分布估计算法[49]，并引用了我们的论文作为这方面研究的原始代表性成果。2009年之后，许多国内外学者在基于copula理论的分布估计算法方面进行了研究，这里先简单介绍几个代表性的成果。文献[50]中提出了copula分布估计算法的基本框架，将概率模型的估计操作分成了两部分：copula函数的估计和各边缘分布函数的估计；在采样操作中首先根据copula函数采样得到概率向量（即该向量的每一维都代表一个概率，其值在(0,1)区间内），然后根据概率向量和各边缘分布函数计算得到搜索空间中的点作为新一代中的个体。文献[51~60]在此基础上对copula分布估计算法做了进一步的研究。文献[61]利用copula函数和一种称为D-vine的规则藤结构在分布估计算法中建立概率图模型，并从理论上分析了三维copula函数的熵和条件交互信息之间的关系。文献[62]利用阿基米德copula函数建立变量之间的两两相关的关系，并提出了一种估计copula函数参数的方法。文献[52]和文献[63]分别利用阿基米德copula函数和经验copula函数反映变量之间的相关结构，提出了相应的基于copula理论的分布估计算法。文献[53]在阿基米德copula分布估计算法中用伪极大似然估计的方式估计copula函数的参数，得到比固定参数的copula分布估计算法更好的优化效果。文献[64]提出了选择copula函数的方法，利用依赖树和二维copula函数表示概率图模型，建立了灵活的联合分布函数。5.2.1copula分布估计算法的统一框架在copula理论中，多变量的联合分布可以分解成两部分：（1）各个变量的边缘分布，（2）一个copula函数。copula函数是多元函数，其中每个自变量取值范围为[0,1]。显然对单变量的边缘分布进行估计和采样要比联合分布简单的多。通过copula函数“连接”各边缘分布可以获得所有变量的联合分布。对此联合分布进行采样时，可以首先对copula函数进行采样，产生[0,1]D空间中的点，然后根据各边缘分布函数的反函数求得采样点的位置。根据Sklar定理，随机向量的联合分布可以分成两部分表示，一是各个随机变量的一维边缘分布，另一是反映随机向量相关结构的copula函数。那么在copula分布估计算法中就可以用此定理来完成估计群体分布模型的操作，其分布模型为优化变量的联合分布。事实上，在copula分布估计算法中联合分布不需要显式表示，只需要估计出各个边缘分布Fi(i=1,..,n)并构造出copula函数C即可。在分布模型估计操作完成之后，下一步需要对此模型进行采样。根据copula理论，为产生服从联合分布H(x1,x2,…,xn)=C(F1(x1),F2(x2),…,Fn(xn))的个体，首先需要产生服从联合分布C的[0,1]n空间中的点(u1,u2,…,un)，然后根据边缘分布函数Fi(i=1,..,n)的反函数Fi-1计算得到值xi=Fi-1(ui)，这样得到的向量(x1,x2,…,xn)服从分布H，即(x1,x2,…,xn)为根据估计的概率模型H采样得到的新个体。根据x根据xi(k)=Fi-1(ui(k)),i=1,..,n,k=1,…,l计算得到新的个体当前代群体优势群体选择估计copula函数C(u1(1),u2(1),…,un(1))(u1(2),u2(2),…,un(2))……(u1(l),u2(l),…,un(l))采样F1F2…Fn估计边缘分布是否终止最优个体为优化结果否是图5.1copula分布估计算法的算法原理图综上所述，copula分布估计算法如图7.1所示，首先在搜索空间中按均匀分布产生规模为PS的初始群体，然后循环执行下列4个步骤直至满足终止条件：Step1：选择。在当前群体中按基于适应值的某种选择策略选择s个个体组成优势群体，记为x={xi=(xi1,xi2,…,xin),i=1,2,…,s}。Step2：估计边缘分布。优势群体中的s个体是n维随机向量(X1,X2,…,Xn)的s个样本，那么{xij,i=1,2,…,s}就是随机变量Xj的样本，据此可估计其边缘分布函数Fj(j=1,..,n)。Step3：从copula函数C中采样。产生l个服从联合分布函数C的向量(u1(k),u2(k),…,un(k)),k=1,…,lStep4：产生新一代群体。由以下三部分构成新一代群体：①保留当前代群体中的适应值最好的m个个体至新一代；②通过计算xi(k)=Fi-1(ui(k))(i=1,..,n,k=1,…,l)得到l个新的个体(x1(k),x2(k),…,xn(k)),k=1,…,l；③在搜索空间中按均匀分布随机产生其余的PS-m-l个个体。copula分布估计算法考虑的是所有变量的联合分布，其联合分布用copula函数和各变量的边缘分布表示，采用一般的copula函数可以表示所有变量全部相关的情况。特别地，对于copula函数C=∏，那么变量之间就是相互独立的。采用不同的边缘分布就可以表示变量不相关的一些典型算法。如UMDAc可以认为是copula分布估计算法的一种特殊情况，其copula函数为∏，边缘分布为正态分布。对于变量可分为若干相互独立的组,且每组内变量相关的分布估计算法，则可以用嵌套式的copula分布估计算法表示，其copula函数C=∏(C1,..,Cs)，s表示变量的分组数。copula理论将复杂的联合分布中的边缘分布分离出来，而把注意力集中在了反映变量之间关系的骨架结构上。在copula理论中起重要作用的关键定理是Sklar定理[46]。该定理表明：在连续的情况下，一个联合分布函数可以唯一的表示为一个copula函数和各变量边缘分布函数的复合函数；且对于给定的copula函数及边缘分布函数，由它们确定的联合分布也是唯一的。5.2.2copula分布估计算法的收敛性不失一般性，设优化问题是maxf(x),x∈D，其中x=(x1,x2,…xn)∈Rn，D是Rn中的一个非空有界闭集，而且f(x):D→R是一个正的连续的目标函数。那么存在一个点x*∈D，对任意x∈D，有f(x*)≥f(x)，即x*是优化问题maxf(x)的全局最优解，在此点处目标函数达到最大值G*=f(x*)。设对于任意实数Q<G*，集合B={x|x∈D,f(x)>Q}的Borel测度均为正。设Pop(t)和Pops(t)分别表示第t代的群体和选择群体，它们分别服从联合分布函数H(x,t)和Hs(x,t)。在copula分布估计算法中，选择操作就是从Pop(t)中按一定的策略选出Pops(t)，即：H(x,t)→Hs(x,t)。在估计操作中，对Pops(t)的概率分布模型进行估计，首先估计每一维的边缘分布Fis(xi,t),i=1,2,…,n。由Glivenko-Canteli定理[92]可知，在群体规模无穷大的情况下，第t代选择群体的第i维经验分布函数收敛于Fis(xi,t)，因此Pops(t)的各一维边缘分布函数Fis(xi,t)可以被准确地估计。然后根据Pops(t)估计copula函数Cs(u,t)。根据Sklar定理可知，在D上存在唯一的copula函数Cs(u,t)，使得Hs(x,t)=Cs(F1s(x1,t),F2s(x2,t),…,Fns(xn,t),t)；对于给定的copula函数Cs(u,t)和边缘分布函数Fis(xi,t)，其组合Cs(F1s(x1,t),F2s(x2,t),…,Fns(xn,t),t)一定是分布函数。也就是，在Fis(xi,t),i=1,2,…,n确定的条件下，Hs(x,t)和Cs(u,t)是一一对应的。那么copula分布估计算法的估计操作就可以表示为：根据Pops(t)估计一维边缘分布函数Fis(xi,t),i=1,2,…,n和copula函数Cs(u,t)，设估计值为和，由Glivenko-Canteli定理，。因此，估计操作实际上是根据Pops(t)估计Fis(xi,t),i=1,2,…,n和。在采样操作中，首先根据copula函数采样得到[0,1]n空间中服从分布函数的随机向量(U1,U2,…,Un)的取值，然后根据边缘分布函数的反函数(Fis)-1计算Xi=(Fis)-1(Ui,t)，从而得到新群体(X1,X2,…,Xn)。在copula函数采样无误差的情况下，(X1,X2,…,Xn)服从分布。(X1,X2,…,Xn)构成了第t+1代群体，前面已经说明其分布为H(x,t+1)，其copula函数表示形式为Cs(F1s(x1,t+1),F2s(x2,t+1),…,Fns(xn,t+1),t+1)。由采样操作可知，Fi(xi,t+1)=Fis(xi,t),i=1,2,…,n，根据Sklar定理，。综上所述，copula分布估计算法可以为表示反复执行下面的两个操作：选择：H(x,t)→Hs(x,t)变异：Cs(u,t)→C(u,t+1)设联合分布函数H(x,t)和Hs(x,t)对应的联合密度函数分别为h(x,t)和hs(x,t)，边缘分布函数Fis(xi,t)对应的密度函数为fis(xi,t)，copula函数C(u,t)和Cs(u,t)本身是分布函数，它们对应的密度函数分别是c(u,t)和cs(u,t)。则，其中(5.1)，其中(5.2)定义5.1：设用分布估计算法优化maxf(x)，其最优值为G*，算法中第t代群体服从密度函数h(x,t)，令(5.3)即E(t)表示第t代群体的平均适应值，如果，则称分布估计算法全局收敛。引理5.1：如果选择算子为比例选择，那么一个个体被选中的概率和它的适应值成正比，选择群体的密度函数可表示为(5.4)引理5.2：如果选择算子为截断选择，其选择率为α(t)，那么第t代群体中适应值最好的百分之(100α(t))个体被选中，组成选择群体。从而选择群体的密度函数可表示为(5.5)其中，实数β(t)满足条件(5.6)即第t代群体中适应值不小于β(t)的个体被选中。引理5.3：如果选择算子为二个体锦标赛选择，那么在第t代群体中随机选择两个个体，将其中适应值较大的一个加入到选择群体中。此时选择群体的密度函数可表示为。(7.7)Fatou引理：设{fn(x)}是可测集E上的一系列负可测函数，则(5.8)定理5.1：在群体规模无穷大的情况下，如果copula分布估计算法中的初始群体密度函数h(x,0)为正的连续函数，且copula函数满足C(u,t+1)=Cs(u,t)，那么：（1）采用比例选择算子的copula分布估计算法收敛；（2）如果截断选择算子的选择率α(t)满足：在任意代t，都有α(t)≤θ(0<θ<1)，那么采用此截断选择算子的copula分布估计算法收敛；（3）采用二个体锦标赛选择算子的copula分布估计算法收敛。证明：由于C(u,t+1)=Cs(u,t)，(5.9)故c(u,t+1)=cs(u,t)。(5.10)又在copula分布估计算法中fi(xi,t+1)=fis(xi,t),i=1,2,…,n。(5.11)根据公式(5.1)、(5.2)、(5.10)和(5.11)有(5.12)（1）如果使用比例选择，根据公式(5.4)和(5.12)有(5.13)已知f(x)和h(x,0)在D上是正的连续函数，故h(x,t)也是正的连续函数，且对任意t≥0，有0≤E(t)<G*。因为，(5.14)所以(5.15)从而(5.16)由于h(x,t)在D上是正的连续函数且E(t)≥0，故对任意t>0，有E(t+1)≥E(t)。{E(t)}构成了D上的单调有界序列，故极限存在，记为。下面证明G=G*。用反证法，假设G<G*。由公式(7.13)有(5.17)当f(x)>G时，有(7.18)又因为对任意x∈D，h(x,0)>0，所以对任意满足f(x)>G的x，有。令S={x|x∈D,f(x)>G}，则S的Borel测度非零。从而根据Fatou引理有(5.19)这与h(x,t)是密度函数相矛盾。因此G=G*，copula分布估计算法全局收敛。（2）如果使用截断选择，由公式(5.5)和(5.12)有(5.20)那么对任意f(x)<β(t)有h(x,t+1)=0。从而(5.21)根据公式(5.6)(5.22)比较公式(5.21)和(5.22)可知则极限存在，记为。下面证明β=G*。同样用反证法，假设β<G*。那么对于任意f(x)>β，有(5.23)又因为对任意x∈D，h(x,0)>0，所以对任意满足f(x)>β的x，有。令S={x|x∈D,f(x)>β}，则S的Borel测度非零。从而根据Fatou引理有(5.24)这与h(x,t)是密度函数相矛盾，故。因为E(t)≥β(t)，所以。（3）如果采用二个体锦标赛选择，那么由公式(5.7)和(5.12)有(5.25)对任意ε>0，令，(5.26)，(5.27)，(5.28)由公式(5.25)可得(5.29)显然，对任意x∈Nε，有(5.30)由公式(5.29)和(5.30)可得，(5.31)其中(5.32)由的对称性可知(5.33)那么(5.34)从而(5.35)因为，所以，即。因为，所以，，即。5.3阿基米德copula分布估计算法下面对copulaEDA进行具体分析。设优化问题为。(5.36)记规模为s的优势群体为，则是随机向量的s组观测值。各随机变量的边缘分布可以采用正态分布或者经验分布函数构造。在本节中介绍copula函数为阿基米德copula函数的copulaEDA。根据Sklar定理，由各随机变量的边缘分布函数及copula函数可以确定一个联合分布。此时，解空间的概率模型构建成功。接下来需要对估计的模型进行随机采样产生新种群。在采样中关键的步骤是对copula函数进行采样。5.3.1阿基米德copula函数采样算法设阿基米德copula函数为C，其生成元为φ，(U1,U2,…,Un)是服从联合分布C的随机向量，根据Marshall和Olkin提出的算法[65]，如果存在一个分布函数F，满足F(0)=0，且其拉氏变换与生成元的反函数相等，即，则产生(U1,U2,…,Un)的取值(u1,u2,…,un)可采用如下方法。算法5.1：从阿基米德copula函数中采样的算法Step1：采样,其中是生成元的反拉氏变换；Step2：产生相互独立的值xi~U[0,1],i=1,…,n；Step3：令，则(u1,u2,…,un)为所需要的向量。常见的阿基米德copula函数及其反拉氏变换如表5.1所示，在本节中以Gumbelcopula为例介绍阿基米德copula分布估计算法及其参数估计。表5.1阿基米德copula函数的生成元及对应的的反拉氏变换函数名称生成元反拉氏变换对应的分布函数GumbelClaytonFrank参数为的正整数域内的对数级数5.3.2Gumbelcopula分布估计算法Gumbelcopula函数的生成元是，由表5.1可知，根据Marshall等提出的采样算法，从Gumbelcopula函数和边缘经验分布函数采样新个体的算法如下：算法5.2:(x1(k),x2(k),…,xn(k))=sample_Gumbel(k)Step1：产生服从均匀分布的随机变量取值；Step2：产生与独立的、服从指数分布exp(1)的随机变量取值W；Step3：令,,,,；Step4：通过如下计算得到服从分布的随机变量取值Z：(5.37)Step5：通过如下计算得到服从分布的随机变量取值v：(5.38)Step6：产生相互独立的服从均匀分布U[0,1]的随机变量取值vi,i=1,…,n,令(5.39)Step7：通过如下计算得到新个体(x1(k),x2(k),…,xn(k))(5.40)在上述算法中，(u1,u2,…,un)是服从分布函数为Gumbelcopula函数的In空间中的向量。(x1(k),x2(k),…,xn(k))是搜索空间中服从联合分布C(F1,F2,…,Fn)的向量。在copula分布估计算法框架下，基于Gumbelcopula函数和边缘经验分布函数的copula分布估计算法如下：算法5.3:具有经验边缘分布函数的GumbelcopulaEDAStep1：随机产生规模为PS的初始群体Pg，令当前进化代数g←0；Step2：根据一定的选择策略从当前群体Pg中选择出规模为s的优势群体Sg；Step3：根据Sg估计每一维变量的边缘经验分布函数Fi；Step4：Fork=1tol执行(x1(k),x2(k),…,xn(k))=sample_Gumbel(k)；Step5：由以下三部分构成新群体：Pg中的m个最优个体，l个新产生的个体，在搜索空间中随机产生其余的PS-m-l个个体；令g←g+1；Step6：如果满足终止条件，则算法结束，Pg中适应值最好的个体为优化的结果；否则转到Step2继续执行。实验结果表明了copula分布估计算法的有效性和可行性，并且在优化效果方面有很好的表现。同时在研究过程中发现需要进一步研究的内容：阿基米德copula分布估计算法中关于copula函数的参数对算法的影响，以及如何根据优势群体确定适当的参数。因为copula函数的参数也影响着变量的相关结构，从而影响了所估计的概率分布模型与优势群体实际的概率分布模型之间的吻合度，进一步对copula分布估计算法的性能有所影响。5.3.3基于PMLE的阿基米德copula分布估计算法参数估计法不同的copula函数体现的变量之间的相关结构不同，copula函数的参数不同，所对应的变量之间的相关程度也不同，在copula函数选定的情况下，如何选取参数是一个很重要的问题。不同的参数θ对应的分布函数不同，因此需要根据优势群体确定相应copula函数的参数值，使得该函数能够更准确的模拟优势群体中变量的相关性。在本节中参数估计方法采用PMLE，具体操作如下：Step1：根据样本x={xi=(xi1,xi2,…,xin),i=1,2,…,s}估计边缘经验分布,j=1,2,…,n，即(5.41)其中，x(i)j为第j维向量的s个样本按升序排列后的第i个值。Step2：根据公式(5.41)中的分布函数确定样本x对应的分布函数值u={ui=(ui1,ui2,…,uin),i=1,2,…,s}。Step3：根据u中的值及表7.1中具体的copula函数求得θ的极大似然估计值。Gumbelcopula函数的边缘分布函数仍然是Gumbelcopula函数，其参数θ的值保持不变。因此，当Gumbelcopula函数的维数较高时，可以随机选择两个变量，估计其二维边缘Gumbelcopula函数的参数值θ作为所有变量的联合分布的参数值。对于Gumbelcopula函数，令，,，，则Gumbelcopula函数的密度函数为。在实验中，随机选择两个变量ui和uj，其对数似然函数为，其中。由于θ≥1，可令θ0=1，并取步长为t，θp+1=θp+t。根据样本依次计算，直至，则为极大似然估计的近似值。实验表明这种具有参数估计的copula分布估计算法优于固定参数的copula分布估计算法5.4经验copula分布估计算法经验copula分布估计算法中根据选择群体构造经验copula函数，其原理与一维经验分布函数相同，然后根据其构造原理推导出从经验copula函数中采样的算法。在实现经验copula分布估计算法时，不需要明确构造出经验copula函数，而只需要直接实现从copula函数采样即可。5.4.1多维经验copula函数的构造方式设n维随机向量X，其联合分布为F(x)，copula函数为C(u)，各变量的一维边缘分布为Fi(x)，i=1,…,n。则根据Sklar定理有。产生一组服从分布F(x)的随机数(x1,…,xn)的方法是首先产生服从分布C(u)的一组随机数u1,…,un，然后再由各边缘分布函数的反函数产生（x1,…,xn），其中。设copula函数C(u)相应的密度函数为c(u)，则边缘密度函数(5.42)显然c(u)=cn(u)。为了产生非独立的随机数u1,…,un，需要知道条件分布函数(5.43)设n维随机向量Z，已知其规模为s的样本z=zi=(zi1,zi2,…zin)，i=1,2,..s，首先对随机向量的每一维构造其经验分布函数，然后根据所得的经验分布函数将样本zi对应到其经验分布函数值，这样就得到向量ui，ui的每维取值在区间[0,1]内。从而由随机向量Z的样本得到了随机向量U的样本，其中，向量u∈[0,1]n就是copula函数的自变量。由copula理论，copula函数C(u)实际上是随机向量U的联合分布函数。那么根据概率论的基本理论，产生向量u的方式就是根据U的边缘分布计算出相应的条件分布，然后由条件分布产生向量u。因此首先需要估计U的边缘分布cd(u)，U的样本已知，边缘分布易于估计。具体做法如下：对于随机向量Z的s个样本，分别按维由小到大进行排序，设第d∈{1,2,..n}维排序后的结果是，那么第d维的经验分布函数就是(5.44)由公式(5.44)可以将s个样本映射为s个向量，作为空间[0,1]n中的s个点。将区间[0,1]分成K等份，记为，即。这样，n维单位空间[0,1]n就分解成为Kn个小空间，。统计每个小空间中点的个数，记为。那么密度函数(5.45)令，则。由公式(5.42)，有(5.46)其中，(5.47)则条件分布函数(5.43)就是(5.48)(5.49)其中,。计算公式(5.47)中定义的的方法在[94]中给出，具体如下：算法5.4：计算Step1：根据公式(5.44)将样本点映射为[0,1]n空间上的点；Step2：令数组f(d),d=1,2,…,n中的元素为0；Step3：for(i=1:s)for(d=1:n);endfor(d=2:n) end end产生相互依赖的一组随机数u1,…,un的过程generation(f)如下：算法5.5：generation(f)Step1：在[0,1]区间内按均匀分布产生随机数u1和uStep2：根据公式(7.48)有，其中是满足条件的{0,…,K-1}内的最小整数。Step3：for（d=3：n） Step3.1在[0,1]区间内按均匀分布产生随机数u; Step3.2根据公式(7.49)有，其中是满足条件的{0,…,K-1}内的最小整数。5.4.2经验copulaEDA算法步骤及复杂性分析设优化问题为。基于上述分析，经验copulaEDA的算法步骤如下：Step1:在搜索空间中按均匀分布随机产生PS个个体，确定选择率ss和变异率tt；Step2：根据适应值选择其中的s=ss×PS个主体作为优势群体(z1,…,zs)；Step3：根据(z1,…,zs)计算，即执行算法7.4；Step4：重复下列步骤l=(1-ss-tt)×PS次，产生l个个体加入到群体中；Step4.1:(u1,…un)=generation(f);Step4.2:for（i=1：n） end（x1,…,xn）为产生的新个体。Step5：在搜索空间中按均匀分布随机产生tt×PS个个体加入到群体中；Step6：若满足算法终止条件，则结束，群体中具有最优适应值的个体即优化结果。否则转Step2继续执行。在文献[66]中已证明当s能被K整除时，上述方法构造的函数C(u)是一个copula函数。构造copula函数的时间复杂度是。根据C(u)产生服从该分布的随机向量的时间复杂度是。对于经验copulaEDA，设优化问题的维数为n，群体规模为PS，选择率为ss，则根据选择的群体估计概率模型并采样生成下一代群体的时间复杂度为(5.50)实验表明，经验copulaEDA能够收敛到全局最优解，但算法的局部开采能力较弱，在进化后期收敛的速度较慢。因此可以在经验copulaEDA中引入其他的进化算法以增强其局部开采能力，或者在算法后期采用爬山法等局部搜索法进行改进。5.5基于离散Quasi-Copula的分布估计算法Copula方法是利用随机变量的边缘分布来近似确定其联合分布的数学方法，是在构造多元联合分布以及随机变量间相关结构分析中常用的工具。目前copula理论主要被应用于金融市场的相关性分析、资产定价和风险管理等金融领域。国内对copula理论更深入的研究还不多，且很少利用copula理论研究复杂的组和优化问题。Copula的种类很多，不同的问题应用不同的copula方法，离散Quasi-Copula方法是求解离散变量间相关性的一个有效工具[67-69]。本节根据Quasi-Copula函数在构建反映离散随机变量相关性问题上具有的优势，将其引入分布估计算法中构造反映离散随机数据分布的多变量相关的概率模型，并从中采样产生新个体，从而提高分布估计算法求解复杂的离散优化问题的能力。5.5.1离散Quasi-Copula基本概念Copula是连接、交换的意思,它是把多个随机变量的联合分布与它们的边缘分布相连接的连接函数。Quasi-Copula是copula的一个扩展，它是描述在随机变量的边缘分布函数不可导的条件下，刻画变量间相关性的连接函数。离散Quasi-Copula是描述离散变量条件下的连接函数。n维离散Quasi-Copula可描述如下：设是定义在n维空间In的一个集合，，,…为大于1的自然数。则一个离散Quasi-Copula是一个n元函数C:In→I,I=[0，1]，且满足下列条件：（1），，若中至少有一个为0，则=0.（2）若中除了外全为1，则=.（3），，.5.5.2基于离散Quasi-Copula的概率模型从文献[70,71]中可知，对于给定的n维样本，对应的经验copula为：，其中，…为样本对应的顺序统计量。离散Quasi-Copula的种类很多，很明显，一个经验copula就是一个典型的离散Quasi-Copula。因此，本节以经验copula为例描述个体的相关性，求解联合概率分布，建立概率模型，进而求解问题。确定各变量的边缘分布函数在离散优化问题中，随机获得的样本（个体）的分布函数及其类型都是未知的，因此这里通过顺序统计量获得样本的经验分布函数[72]，再利用经验分布函数估计个体的边缘分布函数。具体过程如下：设优势群体中的个体为Xi=(x1i,x2i,…xDi),i=1,2,…n。按照维数，对所有个体相同维数上的元素按大小递增的顺序重新排列，获得每一维上的顺序统计量，即xd(1),xd(2),…xd(n),d=1,2,…D。针对所有个体的每一维，计算其经验分布函数，获得个体的边缘分布函数，即（5.51）xd.表示第d维上的变量,d=1,2,…D。Fd(xd.)为第d维上的边缘分布函数，是单调非降左连续的跳跃函数（阶梯函数），在每个点上跳跃量都是1/n，0≤Fd(xd.)≤1，Fd(xd.)的可能取值为0,1/n,…n-1/n,1。令ud,i=Fd(xd,i)，即第d维上第i个变量的概率值，ui=(u1i,…uDi),ui为第i个个体Xi的映射点([0,1]D→[0,1])。确定离散Quasi-Copula函数首先对n维空间In进行分割[73]，把空间分成k个子空间，Ω1,Ω2,…,Ωk，Ωj=[((j-1)δ,jδ]，j=1,2,…k,k是一个能被n整除的正整数。则第j个子空间可以表示为，即每个子空间是D维的。分割的细密程度由k决定。令其密度函数为相应点落在对应子空间中的个数的函数，即，，其中，，Nj表示落在第j个子空间上的点的个数。令，，=-1(5.52)则定义各变量的边缘密度为：，d=1,2,…D，(5.53)(5.54)通过边缘密度函数，得到条件分布函数为：(5.55)由上面建立的条件分布函数，获得变量的联合概率分布：(5.56)此即建立的基于离散Quasi-Copula的概率模型，从中抽样产生新个体。5.5.3群体的产生新个体产生分两步[73]，第一步，先从概率模型中抽样产生相关的随机数，然后以此随机数为基础产生新个体，组成新群体。Step1在每一维上，产生相关的随机数（v1,v2,…,vn）首先，在[0,1]之间产生一个随机数，然后按下列公式产生后面的数，(5.57)↓jd是满足式(5.9)的最小整数，↓jd∈{0,1,…k-1}(5.58)Step2产生新个体新一代个体的产生如下：，，，，t为进化代数。5.5.4算法流程算法采用经验分布函数确定各变量的边缘分布函数，然后通过经验Copula方法，通过从样本中估计变量的联合分布函数，获得多变量相关的概率分布模型，并从中抽样产生新群体。具体步骤如下：Step1编码。采用自然数编码方式对问题进行编码.Step2产生初始群体。设群体规模为N。以均匀分布随机产生N个自然数排列作为初始群体。Step3适应度计算和优势群体的选择。Step4建立概率模型。Step5产生新群体。Step6更新概率模型。从新产生的群体中，选择新的优势群体，重新计算个体的联合分布，获得新的概率模型。判断收敛条件，如果收敛或满足终止条件，结束。否则返回Step3。重复上述步骤，直到满足收敛条件。整个过程采用精英保留策略。5.5.5实例仿真为了验证本节提出的基于离散QuasiCopula的分布估计算法（CEDA）的优化性能，针对TSP问题进行仿真实验，测试实例来自TSPLIB数据库。实验结果与算法GA、PBIL和OPBIL[74]的实验结果进行了比较。实验环境为：种群规模等于城市个数，eil51和berlin52问题最大运行2000，eil76最大运行2500代，其余例子最大运行3500代，参数=0.15（PBIL与OPBIL参数设置相同）。优势群体选择比例为0.5。每个例子运行30次，为运行30次所得解的平均值，为均方差。仿真结果如表5.2所示。表5.2CEDA算法与其它算法计算结果的比较例子算法GAμσPBILμσOPBILμσCEDAμσeil51687.9340.87590.1730.6546.0320.94521.1119.91berlin5212294.14859.1810676.56812.609792.80548.739670.42534.88eil76986.9974.75846.9457.71754.2344.49738.8342.90kroA10055285.645290.8245871.623769.3238726.172880.9338160.552531.76kroB10056804.895324.8845775.774352.5838300.592782.2938177.562790.95kroC10055560.714637.7245758.083216.8537351.602653.4937141.022346.81kroD10056120.534833.3644964.253324.7637229.582327.2237219.462077.36kroE10055657.734583.1646057.643396.0938876.372637.2638705.052577.17eil1011323.5277.861112.9760.10960.9959.91940.2756.78ch13016566.871573.2214400.64933.6711638.47922.8911397.02687.25从表5.2中可以看出，算法CEDA得到的均值与方差都大大优于GA与基本PBIL算法的结果。OPBIL算法是一种改进的PBIL算法，改进变异算子提高种群的多样性，从而提高算法的求解能力，而算法CEDA的结果优于OPBIL算法的结果，说明算法CEDA有较好的搜索能力。5.6优良模式连接的分布估计算法5.6.1优良模式连接的思想模式在遗传算法中占有很重要的地位，在基于二进制编码的基本遗传算法中，个体是以二进制字符串形式表示的，遗传算法中串的运算实际上是模式的运算[75]。模式定理可描述如下:引理5.6.1模式定理.设遗传算法采用比例选择策略，其中交叉概率和变异概率分别为pc和pm，且取值较小，模式H的阶为o(H)，长度为δ(H)，第t+1代种群P(t+1)中含有H中元素个数的期望为E[|H∩P(t)|],则以下不等式成立：（5.59）它说明了包含适应度值高、长度短、低阶的建筑块的个体在后代中至少以指数级增长。模式定理在一定程度上解释了基本遗传算法的有效性。说明了阶次低、定义长度短且其适应度值超过种群平均适应度值的模式在种群中的数目的期望值以指数级递增。模式定理是遗传算法的基本理论，保证了较优的模式（遗传算法的较优解）的数目呈指数增长，是解释遗传算法机理的有效工具。模式定理为积木块假设理论提供了有力的理论支持。引理5.6.2积木块假设（buildingblockhypothesis）.遗传算法通过低阶、短定义距以及平均适应度值高于种群平均适应度的模式（积木块），在遗传操作的作用下相互结合，最终接近全局最优解。目前，大量的实践支持积木块假设理论，并在许多领域内取得了成功。模式定理保证了较优的模式（遗传算法的较优解）的样本数呈指数增长，即遗传算法存在找到最优解的可能性，而积木块假设指出，积木块在遗传算子的作用下，能生成低阶、短距和高平均适应度的模式，最终生成全局最优解。在进化领域中，基本遗传算法是基于人工选择和交叉、变异等遗传操作构成的一种优化方法。由于遗传操作的随机性，在对大量的构造块（积木块）进行选择和重组过程中易导致构造块破坏和连锁问题的发生，从而导致算法早熟或陷入局部最优解。因此，研究能更好地识别构造块，从而解决连锁问题的算法成为研究的一个热点问题。在求解问题的过程中，如果只是独立的考虑种群中的各个串，则不能充分考虑变量间的相关性，当通过概率把各个串结合起来考虑，发掘串群体的相似点，我们就会得到大量的串信息来帮助指导搜索。因此，借鉴模式定理与积木块假设的思想与结论，我们提出了一种优良模式连接的思想。主要内容可描述如下，采用自然数编码方式，并令最初的模式只包含两个相邻的确定值，在进化过程中，通过在规模不大的优势群体中考虑相似点的大量信息，对以较高频率出现在后代中的相邻值以概率为基础进行连接，组成包含多个相邻确定值的模式，这些模式中的确定值组成了一个连接块，令其为优良模式连接块，并在进化过程中以块为整体参与进化，使得那些适应度值高于种群平均适应度值的模式在下一代中将会有更多的代表串，而对那些适应度值低于种群平均适应度值的模式，它们在下一代中的代表串将会减少，充分体现了“优胜劣汰”的思想。本节将优良模式连接的思想引入分布估计算法中，改进算法中概率模型的建立方式，为提高算法的优化性能提供了可能。5.6.2模式矩阵的建立模式矩阵是根据个体中两两相邻变量在优势群体中出现频率的大小建立的。为表述方便，以模式中的变量（确定值）表示相应的模式。记初始群体为，初始优势群体为，优势群体规模为n。其中第k个个体表示为，。以为模板建立模式矩阵。由于初始群体以均匀分布产生，因此，令所有变量产生的初始概率为，通过统计中两两相邻变量出现的频率，建立模式矩阵如下：t为进化代数,这里t=1。且，（5.60）（5.61）其中，，为1表示中的第k个个体中存在相邻变量，否则为0。为优势群体中相邻变量出现的频率。公式（5.60）中，越大，上一代对下一代影响越大，反之越小。在建立模式矩阵时，只考虑此相邻变量是否出现在相应的个体结构中，这种方法能增大出现频率高的模式被保留的机会。5.6.3分块优化过程建立分块模型设进化到第s代，对应的模式矩阵为，根据对最优个体结构进行分块，建立分块结构模型如下[76]：最优个体中的模式信息代表优良模式信息。若某相邻变量的模式出现的频率高，则其出现在下一代中的几率大。若存在若干相邻变量，出现的频率均很高，为避免在进化过程中优良连接模式被破坏，同时避免在同一区域的重复搜索，则使其连接成一个块，以块为单位参与进化。这时可设定一个参数，若多个相邻模式出现的频率大于，令其连接成一个块。越大，表示被连接成块的模式出现的频率越高，反之越小。具体分块步骤如下：设最优个体为，且∧考虑相邻变量组成的模式的特征，若对应模式矩阵，的值大于参数，则认为模式与相关，属于同一个块，否则，相互独立。越小，连接条件越易满足，变量越易被连接成子块，反之，变量越难被连接成子块。在优化过程中，每个子块被作为一个整体参与进化。若太小，很容易满足连接条件，多个变量很快被连接成子块，使得算法易陷入局部最优。太大，连接条件不容易满足，算法总在整个模式中搜索，降低搜索速度，不利于搜索到最优解。一般来说选择。产生初始子块群体设最优个体最大分块数为R，R为正整数。随机产生M个1到R的排列作为子块初始群体，根据子块排序计算个体得分，选择最优个体和优势群体，令子块优势群体规模为m。分块概率模型的建立和更新令第t代子块优势群体为，第k个个体表示为，。t=0表示初始子块群体，初始子块优势群体记为。由于各子块间相互独立，且可任意交换位置，因此，在子块优势群体中，针对每一个位置，计算各子块排在该位置的概率，建立分块概率模型如下：其中，为第个子块排在第个位置上的概率，t为进化代数，。由于初始子块群体以均匀分布产生，因此令子块的初始产生概率为。通过计算初始子块优势群体中各子块排列位置的概率，得到为：（5.62）在进化过程中，通过计算新的优势群体中各子块排列位置的概率，对概率模型进行更新，概率的更新如下：（5.63）（5.64），越大，上一代个体，对下一代影响越大，反之越小。与意义相同，取值可相等也可不等，文中取二者相等。产生新子块群体从新的概率模型中重新抽样。抽样方式为以的每一列为基础，采用轮盘法反复抽样，给每个位置选择取值，产生新子块群体。计算子块个体得分，判断收敛性。若收敛则结束；否则，调整子块内部变量间的排序。调整子块内部各变量的排序为避免陷入局部最优，针对每一个子块，随机或以条件概率为基础产生新排列，重新调整子块内部各变量的排列顺序，进一步进行局部搜索。重复执行上述步骤，得到子块的若干新排列。计算得分，若得分有改进，保留新排列，否则原排列不变，执行下一个子块操作。直到所有子块操作完成。判断收敛性，若满足收敛条件，结束。否则，选择新的子块优势群体，返回步骤，更新子块概率模型。重复步骤至，当最优子块个体在若干代内得分不变时，分块优化过程结束。5.6.4算法流程基于优良模式连接的EDA算法的主要内容包括，首先，从种群中选择优势群体，统计所有两两相邻模式出现的频率，构建模式矩阵描述个体结构的特征，出现频率越高的模式，在模式矩阵中所占比例也越大。然后，从模式矩阵中反复抽样，产生新的优势群体，根据相应的模式矩阵对优胜个体结构进行分块，构建优良模式的连接块。令各个块之间相互独立，采用变量无关的分布估计算法PBIL优化连接块的排序，针对每个块内部的模式，有条件的进行局部调整，增强局部搜索能力，最终获得问题的满意解。算法基本流程如下：Step1产生初始群体.Step2计算个体的适应度.Step3整个群体按照适应度值从大到小排序，以一定的比例选择若干适应度好的个体作为优势群体.Step4根据优势群体信息，建立模式矩阵.Step5从模式矩阵中抽样产生新群体。计算个体适应度，判断终止条件。若满足终止条件则结束，否则选择新的优势群体，执行Step6.Step6更新模式矩阵.重复执行步骤Step5与Step6，当最优个体在若干代内得分不变时，执行下一步。Step7分块优化过程.当最优子块个体在若干代内得分不变时，当前代分块优化过程结束。令产生的新的子块优势群体与执行分块优化过程前的优势群体中的所有个体，按照适应度的大小重新排序，选择新的优势群体，返回Step4.重复执行上述过程，直到满足终止条件。为保证算法的收敛性，整个过程采用最优保留策略。5.6.5实例仿真为了验证算法的有效性，采用几个著名的JobShop问题和柔性JobShop问题的测试算例作为测试集进行仿真实验。JobShop问题仿真测试本实验采用几个著名的JobShop问题的测试算例作为测试集进行仿真实验。这些算例来自FisherandThompson设计的三个算例Ft06，Ft10与Ft20和其它几类典型算例中的第一个算例。由于本节提出的EDA算法[77]是一种没有混合其它优化算法的单纯EDA算法，因此仿真结果与单纯PSO算法[78]中的结果进行了比较。仿真结果如表5.3所示。表5.3EDA算法与PSO算法结果比较例子问题规模最优解算法PSOEDAMinMaxAverageMinMaxaverageFT066655555956.1555555.0FT10101093098510841035.6937937937.0FT202051165120813521266.9118411841184.0LA01105666666688668.6666666666.0LA06155926926926926.0926926926.0LA112051222122212221222122212221222.0LA1610109459459561035945946945.5LA2115101046110211471128.4107110731071.6LA2620101218126313511312.6125712611257.8LA3130101784178918971830.4178917891789.0LA3615151268137314361409.2129213151312.7从表5.3中可以看出，EDA算法得到的最小值、均值与最大值均小于单纯PSO算法中的值，说明单纯EDA的优化性能优于单纯PSO算法。同时EDA中得到的最大值、均值与最小值之间的差值非常小，说明提出的EDA算法具有很好的稳定性。柔性JobShop问题仿真测试本实验选择了3个有代表性的柔性JobShop例子进行仿真测试。这3个例子均来自文献[79,80]，文中参数与PSO+TS算法[81]参数设置相同。本节所提EDA算法[82]获得的结果与AL+CGA算法,PSO+SA算法[83],PVNS算法[84]和PSO+TS算法进行了比较，具体结果如下。实验1针对问题45的仿真计算问题45是一个完全柔性问题，每个工件的第一个工序的最早开工时间为0，EDA算法所获得的最优解的Gantt图如图5.1所示，EDA算法的仿真结果与其它几类算法仿真计算结果的比较如表5.4所示。图5.1问题45的最优解Gantt图表5.4问题45的计算结果比较AL+CGAPSO+TSEDAMakespan161111Wk101010Wt343232从表5.4中的结果可以看出，EDA算法获得的结果优于混合算法AL+CGA的结果，与PSO+TS算法的结果相同，但是，这里提出的EDA算法是一种单纯算法，算法中没有混合其它启发式算法，而其它两种算法都是混合算法，EDA算法获得的结果优于或等于这两种混合算法的结果，说明提出的EDA算法具有较好的寻优能力。实验2针对问题88的仿真计算问题88是一个部分柔性问题，只能选择某些加工机器，我们令工件的不可选择机器上的加工时间为一个很大的值，例如9999，远远大于其它机器上的加工时间，因此，当选择加工机器时，一般不可能选中这台机器，从而使得问题由部分柔性问题转变为完全柔性问题。问题88的每个工件的第一道工序的最早开工时间为0，EDA算法所获得的最优解Gantt图如图5.2所示，算法的仿真结果与其它算法仿真结果的比较如表5.5所示。图5.2问题88的最优解Gantt图从表5.5中的结果可以看出，EDA算法获得的结果优于混合算法AL+CGA和PSO+SA的结果，与PSO+TS算法和PVNS算法的结果相同，EDA算法是一种单纯算法，算法中没有混合其它启发式算法，仿真结果说明单纯EDA算法有较好的寻优能力。表5.5问题88的计算结果比较算法参数AL+CGAPSO+SAPVNSPSO+TSEDAsMakespan1516151614141514Wk1313121313121212Wt7975757377777579实验3针对问题1010的仿真计算问题1010是一个完全柔性问题，每个工件的第一个工序的最早开工时间为0，EDA算法所获得的最优解Gantt图如图5.3所示，EDA算法的仿真结果与其它算法仿真结果的比较如表5.6所示。图5.3问题1010的最优解Gantt图表5.6问题1010的计算结果比较算法参数AL+CGAPSO+SAPVNSPSO+TSEDAsMakespan77777Wk56-66Wt4544-4343从表5.6中可以看出，这几类算法获得的最优完工时间相同，其它两个参数的结果基本相似，但是本文提出的EDA算法是一种单纯算法，而其它几种算法都是混合算法，从表中的结果看以看出，EDA算法具有较好的寻优能力。5.7基于Bayesian统计推理的分布估计算法Bayesian统计推断理论是数理统计中的一个重要分支。Bayesian统计推断方法与贝叶斯网络优化方法不同，它不需要优化复杂的网络结构，而是直接利用样本提供的信息，通过建立一个后验概率分布对样本进行推断，由于这个后验分布中包含了先验分布信息和样本信息，因而对大样本和小样本均有较好的统计推断效果。因此，本节针对离散优化问题，把Bayesian统计推断方法引入分布估计算法中建立概率模型，提出了一种基于Bayesian统计推理的分布估计算法。5.7.1Bayesian统计推断理论Bayesian统计推断理论是数理统计中的一个重要分支，它有鲜明的特点和独到的处理问题的方法[85]。Bayesian统计推断的主要特点是使用先验分布，在得到样本观测值之后，由x与先验分布提供的信息，得到后验分布。后验分布包含了先验信息与样本提供的信息，是Bayesian统计推断的基础。由于利用了先验知识，Bayesian统计推断对大样本和小样本均有较好的统计推断效果。Bayesian公式为：（5.65）其中，事件构成互不相容的事件完备组，，为先验分布概率，构成先验分布。由于事件B的发生可以对事件发生的概率重新估计，以后验分布体现出来，因此Bayesian公式综合了先验信息与实验提供的样本信息，获得后验信息，并实现了先验分布向后验分布的转化，是Bayesian统计推断的数学模型。5.7.2概率模型的建立很多离散优化问题，尤其是排序优化问题，个体中各个变量的排列位置直接影响个体适应度的大小，因此，建立能反映各变量排列位置信息的概率模型，是分布估计算法求解这类问题的关键。Bayesian统计推断理论为我们提供了一个利用个体中的排列位置信息建立概率模型的有力工具。Bayesian统计推断方法主要是通过建立一个后验分布对样本进行推断，而这个后验分布是利用了先验分布信息与样本提供的条件概率信息，并通过Bayesian公式的转化获得的。因此，在建立后验分布概率模型以前，首先需要建立相应的先验分布概率模型和条件分布概率模型。先验分布概率模型的建立算法采用自然数的排列作为编码方式对离散问题进行编码。充分利用个体中变量的排列位置信息建立先验分布概率模型。设个体长度为L，x=(π(1),π(2),…π(L))为1,2,…L的一个排序，表示一个个体，π(j)为排在第j个位置上的变量。个体中每个变量的排列位置直接影响个体适应度的大小，因此为描述变量排在各个位置上的信息，令个体中的L个变量分别对应L个位置。针对每一个位置上的变量建立一个先验分布概率向量，则所有变量对应的先验概率向量组成一个先验分布概率模型，第j个位置上对应的个体取值与先验分布概率向量的对应关系如下：，概率向量pj=(p1j,p2j,…pLj)T描述了所有变量出现在第j个位置上的先验概率，j=1,2,…L。则所有位置上的概率向量组成一个概率矩阵，即P=(p1,p2,…,pL),称P为所建立的先验概率分布模型。在Bayesian统计推断理论中，先验分布是Bayesian统计模型中的重要组成部分。因此，先验分布的选取是一个很重要的问题。选取先验分布的方法很多，由于本算法是针对复杂优化问题的求解，其解的概率分布都是未知的，需要在求解过程中不断地搜索、更新、重新估计，从而获得一个近似的分布概率模型，从中抽样产生新群体，求得问题的满意解。因此，在此算法中初始的先验分布概率采用均匀分布，这符合先验分布选取中的最大熵原则。熵是信息论中的一个基本概念，是随机变量不确定性的度量，不确定性越大，则熵越大。在“无信息”的情况下，应取熵最大的分布为先验分布，而均匀分布的熵最大，数据随机性最大，包含的信息量也最大。因此，算法中选取初始的先验分布为均匀分布，然后在进化过程中，通过提取优势群体中的优良信息，不断地更新先验分布，从而获得越来越接近最优解的先验分布概率模型。条件分布概率模型的建立若个体中的所有变量之间无先后约束，则每个变量可在所有位置的任意一个位置排列。若对应第j个位置上出现符号（变量）π(j)，而每个变量都有可能出现在该位置上，所以可用一个条件概率向量表示变量π(j)出现在对应位置上的信息，该条件概率向量可表示为qj=(q(π(j)／1),q(π(j)／2),…,q(π(j)／L))T。这个概率向量的值来自优势群体中的样本信息，它描述了当其它变量可能出现在该位置的条件下，变量π(j)出现的条件概率。这时，L个位置对应L个条件概率向量，组成一个条件概率矩阵。设在第t代，第k个个体表示为xkt=(πkt(1),πkt(2),…πkt(L))，第t代优胜个体为xt*=(πt*(1),πt*(2),…πt*(L))。根据优胜个体与优势群体，建立第t代条件矩阵为：其中，qt(πt*(j)／i)表示在符号i出现的条件下，符号πt*(j)出现的条件概率。根据条件概率公式，qt(πt*(j)／i)=qt(πt*(j)·i)／qt(i)成立，其中，q