基于两步优化策略的结构网格负载平衡

2021-10-20段焰辉周鑫高诗婕茜林锦星张怀宝王光学

航空学报 2021年9期

段焰辉，周鑫，高诗婕茜，林锦星，张怀宝，王光学,*

1. 中山大学系统科学与工程学院，广州 510006

2. 中山大学智能工程学院，广州 510006

3. 中山大学航空航天学院，广州 510006

计算流体力学(Computational Fluid Dynamic，CFD)在航空、航天、航海、能源等领域都发挥着重要作用。随着求解问题难度的增加和求解精度要求的提高，CFD的计算量急剧增加，计算周期随之延长。并行计算能够有效地缩短计算周期，其原理是将计算网格分组并分配给高性能计算机的各个进程，进程之间通过传递网格面上的流场信息实现信息交互。分配给各个进程的网格量均匀性以及通信相关的面网格量均匀性是影响并行计算效率的关键，负载平衡是解决网格量和通信相关面网格的分配问题。中国已启动国家数值风洞工程[1]以发展完全自主知识产权的计算流体力学软件，该软件为不同类型网格块的数据交换设计了通用的通信模型及粗粒度MPI/OpenMP(Message Passing Interface/Open Multi-Processing)混合并行模型，因此适用于任意网格类型，且具有良好的扩展性和并行效率[2]。网格负载平衡问题是影响CFD并行计算的另一重要因素，因此也是国家数值风洞工程的重要研究内容。

用于计算流体力学的负载平衡方法大致可分为两类：几何剖分法和图剖分法[3]。几何剖分方法包括Berger和Bokhari提出的递归坐标对分法(RCB)[4]，以解决自适应网格优化策略为目的；还包括RCB的推广，如递归惯性对分法(RIB)[5]和空间填充曲线法[6]等。图剖分法[7]一般用权重图表示计算和通信，包括多层谱对分方法(MSB)[8]、多层剖分方法[9]、超图剖分法[10]以及多约束多目标剖分方法[11]。其中图剖分算法发展迅速，已形成多种应用软件，如ParMETIS[12]和Zoltan[13]。对于多块结构网格，刘宏康等将图剖分用于结构重叠网格的负载平衡[14]；李桂波和杨国伟采用遗传算法对多块结构网格的分配进行研究，但其目标函数仅考虑计算时间，仅按通信时间选择结果[15]；唐波以通信时间为指标、计算时间为约束构造目标函数，采用遗传算法对网格分配进行研究[16]。综上，对于负载平衡问题多采用优化思想求解，但所用优化算法均有局限，如贪婪算法[17-18]效率较高且能在优化过程中剖分大块网格，但不考虑通信时间的平衡；遗传算法能通过设计目标函数的方式同时考虑计算时间和通信时间的平衡，但其结果依赖于目标函数的构造，而且遗传算法本身不能对网格进行剖分。

混合优化策略能结合不同优化算法的优势，并在一定程度上克服不同算法的缺点。混合优化策略已有较多研究与应用[19-20]，其中两步优化策略构造的混合优化方法[21]充分结合两步优化方法的优点，获得了高效的优化算法。本文将贪婪算法和遗传算法相结合，构造两步优化算法：第1步优化采用传统的贪婪算法完成对大块网格的剖分和以进程计算时间为指标的网格块分配；第2步采用遗传算法，目标函数兼顾进程计算时间和通信时间，在第1步优化结果的基础上对网格块在进程上的分配开展二次优化。两步优化策略克服了遗传算法无法对大块网格剖分以及贪婪算法不能考虑通信时间的缺陷。

本文首先介绍了负载平衡算法，包括两步优化策略、用于计算时间和通信时间的模型构建方法和遗传算法，然后采用3个算例对提出的两步优化策略进行验证，最后根据验证结果对两步优化策略的优缺点进行总结。

1 算法介绍

1.1 数学模型及求解思路

结构化网格的负载平衡主要处理两个方面问题：网格块在进程上的分配问题和大网格块的切割问题。在贪婪算法的分配过程中实现对结构化网格负载平衡中的网格块切割问题，因非本文重点，不予赘述，相关方法可参考文献[16]。对于网格块在进程上的分配问题，用如下假设进行处理：设负载平衡问题的进程个数为Np、网格块个数为Nb，满足Nb≥Np。以M表示某种分配方式，第i个进程处理结构网格计算问题的总用时可表示为

(1)

(2)

负载平衡问题的目标函数即是所有进程的最大总用时：

(3)

求解结构化网格的负载平衡问题即是最小化式(3)的最优化问题：

minf(M)

(4)

该优化问题的设计变量是网格块在进程上的分配方式，是一种组合优化问题。

传统求解结构化网格负载平衡问题的方法是贪婪算法，这种方法在优化过程中仅考虑网格单元的数量，没有考虑通讯相关的面网格数量，因此所得结果在大部分情况下都不是最优解。在贪婪算法的基础上，再引入用于组合优化的遗传算法，通过合理组织两步优化策略以克服传统算法未考虑通信时间的缺陷。具体而言，优化策略可分为两步：

1) 采用贪婪算法对初始结构化网格开展初步负载平衡优化。主要完成两方面的内容：对大块网格的分割和以计算时间为指标的网格块分配。

2) 采用遗传算法对第1步的优化结果进行二次优化。将组合优化的遗传操作、引入第1步优化结果的种群初始化方法、综合进程计算时间和通信时间的目标函数计算方法相结合，形成适用于结构网格负载平衡的遗传算法。

1.2 计算时间和通信时间的建模方法

采用统计学的回归思想建立时间计算模型，需要重点解决三方面的问题：如何生成样本、如何建模和如何对模型进行验证。选择恰当的计算网格，通过逐渐增加网格量获取计算时间随网格单元数量以及通信时间随通讯相关面网格数量的样本，同时要充分考虑计算平台硬件和操作系统的影响；建模需要根据样本的分布规律合理选择适当的模型；模型验证通过与真实的模型计算时间对比以检验模型的准确度。

1.2.1 样本生成与预处理

CFD并行计算的计算时间和通信时间都会受计算平台硬件和操作系统的扰动影响。若影响为小量，则通过统计平均消除扰动；若为大量，则作为奇异点处理。通过对多个具有相同硬件配置节点的计算结果对进程数求平均的方式，消除硬件和系统误差；采用先以CFD软件计算多步再对步数求平均的方式，消除操作系统和计算软件自身的影响；重复采样多次，对采样次数做平均以进一步消除硬件、操作系统和计算软件的影响。

j=1,2,…,Nsamp

(5)

式中：Nsamp为每次采样时样本的数量。另需要注意的是，式(5)的计算方式假定每个块的计算量是相同的，但实际情况是物理边界会造成计算量的微小偏差，因为用于计算物理边界的网格量相对于整体网格量很小，因此可以忽略物理边界造成的计算量偏差。

j=1,2,…,Nsamp

(6)

图1 采样网格及边界条件设置示意图Fig.1 Schematic diagram of sampling grid and boundary conditions

图2给出了10次采样的结果。计算时间的结果显示出双线性特点，经过后期多次测试，数据都呈现了相同的趋势。推测应与计算平台有关，但具体原因分析不在研究范畴内，且所有计算均使用相同计算平台，只要建立相应的模型就不会影响后期模型的计算精度。通信时间的结果整体呈线性分布，但部分样本点存在较大偏差，属于奇异点。假设多次采样的数据服从正态分布，采用正态分布数理统计的方法对采样结果进行统计分析。采用1.96σ的置信区间，其中σ为标准差，利用在置信区间内样本点的平均值替换奇异点[22]。处理后的数据点如图3所示。

图2 10组采样数据Fig.2 Ten sets of sampled data

图3 处理奇异点后的采样数据Fig.3 Sampled data without singularities

1.2.2 线性回归建模

鉴于样本数据都呈现线性分布，采用线性回归的思路建立线性模型。对于体网格和计算时间模型分别建立单线性模型和双线性模型，并对结果进行对比。单线性模型的结果如式(7)所示，模型与真值的对比如图4所示。

图4 单线性模型与样本点对比Fig.4 Comparison between single linear model and sample points

T(xv)=4.342 2×10-6xv-3.518 2×10-2

(7)

式中：T为计算时间；xv为体网格单元数量。

双线性模型的结果如式(8)所示，模型与真值的对比如图5所示。

(8)

由图4可知，单线性模型主要拟合了第2段数据，所以最大误差出现在第1段数据的末端，对第1段数据的拟合程度很差。图5双线性模型的拟合精度明显提高，最大误差位置仅出现在整个数据的前端。表1给出了两种模型的误差对比，双线性模型的均方误差和最大误差均优于单线性模型，因此选择双线性模型。

表1 单线性模型和双线性模型的精度对比Table 1 Accuracy comparison between single linear model and double linear model

图5 双线性模型与样本点对比Fig.5 Comparison between double linear model and sample points

同理建立通信时间的线性模型，如式(9)所示，均方误差为1.373 5×10-12，最大误差为1.366 0×10-6。模型结果与样本的比较如图6所示。

图6 通信时间的模型与样本点对比Fig.6 Communication time comparison between model and sample points

C(xs)=4.485 7×10-9xs-4.875 3×10-7

(9)

式中：C为通信时间；xs为面网格单元数量。

1.2.3 模型验证

采用阻力预测会议的标模DPW2的网格[23]，该网格共有141块，网格量约1千万，如图7所示。在天河二号上，采用40个进程计算100步，通过设置并行节点的使用保证计算节点与模型建立时的节点相同。使用计算模型式(2)、式(8)和式(9)分别计算各个进程上的总时间、计算时间和通信时间，并和准确时间对比。

图7 DPW2网格块Fig.7 Grid block of DPW2

相对误差如表2所示。从平均误差看，进程计算时间的准确率高于通信时间，但总体上两者都保持在10%附近；从最小误差、最大误差和标准差能够看到各个进程的误差分布情况，各个进程计算时间和总时间误差的最大、最小值相差不大，标准差也较小，所以整体分布比较均匀；通信时间误差分布的最大、最小值相差较大，标准差也大，所以分布的差异性较大。计算时间和总时间相近的原因是通信时间的数值很小，对总时间的影响较小；通信时间的差异推测和计算平台的性能相关，应与采样时通信时间存在奇异点的原因一致，具体的误差分布如图8所示，较大误差的进程数较少，但存在误差值很大的进程，类似于奇异点。总体而言，建立的计算时间和通信时间计算模型是可信的，进一步由两个时间结合得到的进程总时间模型也是可信的，可用于遗传算法的目标函数计算。

表2 单个进程的模型计算结果与真实结果误差对比

图8 通信时间误差分布Fig.8 Communication time error distribution

1.3 第1步优化：贪婪算法

贪婪算法是求解结构网格负载平衡问题的常用算法，只给出该方法的思路，算法细节可参考文献[16]。贪婪算法的基本思想是“优先匹配优值资源”，具体到结构化网格的负载平衡就是选择大网格块匹配计算能力强的进程。使用计算平台各个进程的计算能力相同，通常用总网格量与进程数量的比值作为进程计算能力的度量。至此，可给出贪婪算法的大致求解步骤：

1) 计算进程的初始计算能力，即总网格量与进程数量的比值。

2) 从未分配网格块中选择网格量最大的网格块，从未分配进程中选择计算能力最优的进程，判定两者是否匹配(网格量和计算能力数值上近似相等)。

3) 匹配则将网格块分配给相应进程，并将该进程的计算能力标定为0。

4) 不匹配则分两种情况处理：第1种，所选网络块网格量大于进程计算能力，则对网格块进行分割，分出匹配该进程计算能力的网格块按步骤3)处理，剩余网格块作为未分配网格块；第2种，所选网格块网格量小于进程计算能力，则把网格块分配给进程，同时使用进程的计算能力减去网格块网格量的差作为进程新的计算能力。

5) 判断是否所有进程的计算能力都为0,是则计算完成，否则转步骤2)。

1.4 第2步优化：遗传算法

采用与文献[15]中一致的遗传算法，但在处理细节上略有不同。用于结构网格负载平衡的遗传算法处理的是将网格块在各个进程上分配的问题，本质上是一个组合优化问题。遗传算法需要针对结构网格分配的特点，解决编码、交叉、变异和种群初始化的问题。以10个网格块分配到4个进程上的问题为例，对上述方法依次展开介绍。

1.4.1 编码

采用两段编码。第1段是所有网格块序号按某个顺序的排列，第2段是分配到各个进程的网格块的总数。如图9所示，其中第1个进程分配1块网格，网格编号为7；第2个进程分配2块网格，网格编号为8和1；第3个进程分配5块网格，网格编号为10、3、5、4和9；第4个进程分配2块网格，网格编号为2和6。采用这种两段编码方式需要满足如下约束条件：第1段编码中的网格序号不允许有重复；第2段编码中的数字之和等于总的网格块数，也即第1段编码的长度。

图9 遗传算法编码示意图Fig.9 Schematic diagram of genetic algorithm coding

1.4.2 交叉操作

交叉操作仅对第1段编码操作，图示不画出第2段编码。采用Goldberg[24]提出的部分匹配交叉方法，该方法能够保证交叉后的编码满足不重复的约束条件。选择两个父代编码如图10所示，图中三角形标识的位置为随机确定的两个交叉位置5和7，交叉操作可分为3步进行。

1) 交换两个父代交叉位置之间的编码，如图11(a) 所示。根据交换的编码确定映射关系，如图11(b)所示。

2) 在子代1的剩余编码中，保留不在父代1交叉位置之间的编码，对子代2则保留不在父代2交叉位置之间的编码，如图11(c)所示。

3) 利用图12的映射关系对剩余编码进行处理，如子代1第1个“*”可通过将父代相同位置的8换成3得到，由此得到最后的交叉结果，如图11(d) 所示。

图10 编码交叉位置示意图Fig.10 Schematic diagram of coding cross position

图11 交叉操作Fig.11 Cross operation

在交叉操作中特别需要注意的是递归问题，仍然以图10为例，将交叉位置改为4和7，如图12(a) 所示。按照上述3步操作得到图12(b)所示结果，此时两个子代都出现了重复编码10(红色部分)。这是因为在用交叉位置之间的编码建立的映射关系中存在递归关系10和5、10和2。

图12 交叉操作递归问题Fig.12 Recursion problem in cross operation

以图12(a)中父代的第9个编码2为例，将其通过映射换为10之后，需要检查10是否在父代2的映射编码之内，如果有就再做一次映射，得到5，即是所求。综上，对于递归问题，需要做两次或以上的映射，直至映射结果不存在递归关系。对图12(a)按递归问题处理得到的结果如图12(c)所示，仍满足编码的约束条件。

1.4.3 变异操作

采用交换算子进行变异操作，对两段编码分别进行变异操作。以第1段编码为例，假设父代编码如图13(a)所示，随机选取的变异位置如图13(a) 中三角形标记所示。通过变异操作，可得子代如图13(b)所示。第2段的变异操作与第1段相同，不再赘述。显然，换位变异操作能够保证编码的第1段和第2段满足约束条件。

图13 变异操作Fig.13 Mutation operation

1.4.4 种群初始化

发展的遗传算法作为二次优化方法需要在第1步优化的基础上开展优化，同时还要考虑第2次优化的设计空间，因此采用随机初始化和第1步优化结果混合的初始化种群。

对第1段和第2段编码分别进行随机初始化。对于第1段编码，对网格块编号采用随机不放回抽样即可得到随机编码。对于第2段编码，采用多次随机拆分的思路：先将网格块总数随机拆分为两个整数，选择其中最大的再次拆分，此时共得到3个整数，再对其中最大数做随机拆分，直到拆分得到的整数个数与负载平衡的进程数相等。

随机初始化能够得到随机种群，但无法利用第1步贪婪算法的优化结果。将遗传算法的初始种群分为两部分，一部分使用随机初始化种群，另一部分使用贪婪算法的结果。

2 算例验证

2.1 特殊结构网格算例

专门设计特殊结构网格算例以说明贪婪算法的问题以及两步优化的必要性。该算例初始网格由10块网格组成，网格块的编号如图14所示。大块网格的网格量为50.0万，小块网格的网格量为12.5万。将上述网格分配到4个进程上。

图14 算例1的初始网格及编号Fig.14 Initial grid and number of Case 1

第1步优化采用贪婪算法，得到的结果如图15 所示，4种颜色分别代表分到4个进程的网格块。可知，第1、2两块分配合理，第3～10块分配仅考虑计算量平衡，未考虑通信时间平衡。由式(2)计算的目标函数值为2.141 77 s。

图15 算例1的第1步优化结果(贪婪算法)Fig.15 1st step optimization result of Case 1 (greedy algorithm)

第2步优化采用遗传算法，种群选取100个，交叉概率0.6，变异概率0.1，初始化种群中贪婪算法优化结果占比10%，优化100代，所得结果如图16所示。可知，优化结果是最优的，此时的目标函数值为2.141 65 s。虽然计算时间缩短很小，但是对于网格量大、迭代步数多、计算状态复杂等计算量巨大的情况效果更为明显。

另一方面，贪婪算法对网格块的排列敏感。如把图15中的网格块编号改为如图17所示，则分配结果与图16相同，也可以得到最优结果。但划分网格时很难兼顾网格块的排列顺序，且当网格块数量巨大时，人工排列网格块顺序几乎难以实现。所以对于计算量巨大的问题，第2步基于遗传算法的优化是必要的。

图16 算例1的第2步优化结果(遗传算法)Fig.16 2nd step optimization result of Case 1 (genetic algorithm)

图17 改变网格块排列顺序的初始网格Fig.17 Initial grid changing grid blocks order

2.2 DPW2算例

DPW2的网格与1.2.3节的网格相同，分配到12个进程上，网格如图7所示。第1步优化采用贪婪算法，得到的目标函数值为3.598 49 s。第2步优化采用遗传算法，种群选取2 000个，交叉概率0.6，变异概率0.1，初始化种群中贪婪算法优化结果占比10%，优化200代。优化结果为3.596 0 s。

第1步和第2步优化结果的网格块分配变化体现在第8、11和12个进程上，具体的变化如表3 所示，3个进程上分配的网格块进行了交叉调整，如表3中“*”标示的网格块编号。调整后各个进程的计算时间如表4所示，计算时间最长的进程由进程8变为进程2，发生变化的主要原因来自进程8计算时间的缩短。

表3 DPW2网格第1步和第2步优化结果的网格块分配对比

表4 DPW2网格第1步和第2步优化结果的时间对比

DPW2的结果说明即使经过贪婪算法的优化，网格块的分配也不能保证计算时间指标达到最优，而遗传算法能够对第1步优化结果进行有效改善。

2.3 DPW6算例

为进一步检验本文方法对大规模网格的处理能力，采用网格量接近1亿的DPW6网格。该网格包含203个网格块，共有近1亿的网格单元，分配到60个进程上，网格如图18所示。第1步优化采用贪婪算法，得到的目标函数值为6.817 21 s。第2步优化采用遗传算法，种群选取8 000个，交叉概率0.6，变异概率0.1，初始化种群中贪婪算法优化结果占比10%，优化200代。优化结果为6.817 3 s。

图18 DPW6网格块Fig.18 Grid block of DPW6

第1步和第2步优化结果的网格块分配变化体现在第30、33、36、37和42个进程上，具体的变化如表5所示，5个进程上分配的网格块进行了交叉调整，如表5中“*”标示的网格编号所示。调整后各个进程的计算时间如表6所示，计算时间最长的进程由进程36变为进程46，发生变化的主要原因来自进程36计算时间和通信时间的缩短。