基于兩步優化策略的結構網格負載平衡

段焰輝，周鑫，高詩婕茜，林錦星，張懷寶，王光學,*

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計算時間和通信時間的縮短。

表5 DPW6網格第1步和第2步優化結果的網格塊分配對比

表6 DPW6網格第1步和第2步優化結果的時間對比

分析DPW6算例的結果，編碼長度的增加極大地提高了尋優難度，即使種群增加到8 000個，目標函數才獲得0.000 1 s的收益，這是遺傳算法尋優的不利之處，也是之后需要重點研究的方向之一。

3 結 論

1) 設計的進程計算時間和通信時間的建模方法具有一定的通用性，可以擴展到相關負載平衡問題的研究上；建立的進程計算時間模型和通信時間模型都具有較高的精度，滿足遺傳算法優化設計的需要；通信時間受計算平臺的擾動較大，所以計算精度略低于網格計算時間的模型。

2) 提出的兩步優化策略能夠有效求解負載平衡問題。第1步優化采用貪婪算法，能夠完成大網格塊的剖分和以計算時間為依據的網格塊分配；第2步優化采用遺傳算法，在第1步優化基礎上，以建立的計算時間和通信時間計算模型計算目標函數開展二次優化，能夠得到更優的分配結果。

3) 第2步優化結果顯示，單進程單步的時間收益是較小的。對于計算量巨大的問題，如網格量巨大、迭代步數冗長和計算狀態很多的問題，整體上能夠獲得較大的時間收益；對于計算量較小的問題，可不開展第2步優化。

4) 本文屬于探索性研究，要提升本文兩步優化策略的工程實用能力，還可從以下問題尋求突破：網格塊數和進程數增多造成編碼長度的增加，從而導致設計空間急劇增大，最終造成遺傳算法尋優困難的問題；組合優化遺傳算法交叉和變異操作的改進問題；非組合優化遺傳算法的改進措施在組合優化遺傳算法上的適用性研究。