采用花朵授粉算法的等幾何邊界元形狀優化

高翔，王林軍，*，劉洋，陳保家，付君健

(1.三峽大學 水電機械設備設計與維護湖北省重點實驗室，宜昌 443002；2.三峽大學 機械與動力學院，宜昌 443002)

固體力學中常用的離散方法有有限差分法、有限元法、邊界元法、有限體積法、基面力元法[1]，巖土力學中常用的離散方法有離散元法和非連續變形分析，流體力學中常用的離散方法有光滑粒子法、移動粒子半隱式法，斷裂力學中常用的離散方法有擴展有限元法。在固體力學的形狀優化問題中，邊界元法僅需要對邊界進行離散即可，比其他方法更適用于形狀優化問題。

由于CAD 軟件采用樣條曲線建模，而CAE 軟件采用網格建模，樣條曲線在離散過程中會產生較大誤差，導致CAD 模型和CAE 模型的結果存在一定差異。非均勻有理B 樣條（non-uniform relational B-splines, NURBS）曲線靈活多變且精度高，通過等幾何分析將CAD 與CAE 這2 種模型無縫連接后，CAE 的精度會有大幅提升[2]。Simpson 等[3]對二維等幾何邊界元法（isogeometric boundary element method, IGABEM）進行了研究，并與標準邊界元法進行對比。Hao 等[4]通過等幾何有限元法對結構進行分析，并提出了一種基于增強步長調整 (enhanced step length adjustment，ESLA) 法、二階可靠性分析法（second order reliability method, SORM）和加速的序列優化與可靠性評估 (stepped-up sequential optimization and reliability assessment，SSORA)法 的 可 靠性優化設計方法。Lian 等[5]提出了一種基于T 樣條曲線的形狀優化設計方法。Yoon 和Cho[6]通過等幾何邊界元法對結構進行分析，并對線彈性問題的邊界積分方程的敏度進行了推導。Lee 等[7]通過等幾何分析、邊界積分和水平集法提出了一種形狀優化設計方法。Choi 和Cho[8]提出了一種用于等幾何形狀優化的內部控制點更新方法。Lüdeker 等[9]提出了一種逆均勻化結構的等幾何形狀優化設計方法。Kang 和Youn[10]提出了基于修剪樣條曲面的殼結構形狀優化方法。Sun 等[11]通過罰函數法及粒子群算法優化了結構的形狀。由于罰函數法的收斂性依賴于罰因子的初始值，亟需尋找一種收斂性優于罰函數法的計算方法。增廣乘子法的收斂性不依賴于罰因子的初始值，收斂性比罰函數法更好[12]。

梯度優化算法的收斂速度很快，但容易陷入局部最優解，而智能優化算法的全局尋優能力更好。為了提高智能優化算法的全局尋優能力，很多學者提出了不同的改進方法。其中，文化算法能通過雙層進化策略建立的知識庫生成新解[13]，云模型能通過熵及隸屬函數生成服從正態分布的新解[14]，分布估計算法（estimation of distribution algorithm，EDA）根據當前若干最優解的均值及標準差生成新解[15]，基于Alopex 的進化算法（Alopex-based evolutionary algorithm，AEA）通過Alopex 算法確定當前變量步長的正負號[16]，精英反向學習（opposite-based learning ,OBL）策略可生成反向解并逐步縮小搜索范圍[17]，量子算法通過量子比特編碼提高算法的遍歷性[18]，混沌算法[19]、Lévy 飛行[20]、隨機游走[21]這3 種策略均可以提高算法的全局尋優能力。鑒于文獻 [13-21]中算法優良的尋優能力，本文將同時采用精英反向學習策略及EDA 改進優化算法。部分學者采用Copula 理論降低樣本的相關性，提高算法的全局尋優能力[15]，但當相關矩陣非正定時，Copula 理論改進的EDA 不再適用。因此，亟需尋找一種實用性更強的EDA。

鑒于增廣乘子法和大規模分布估計算法（large scale EDA，LSEDA）這2 種算法對于優化算法的收斂速度有顯著的改進效果，本文提出一種基于花朵授粉算法和等幾何邊界元的形狀優化算法。該算法有如下優點：①采用邊界元法對結構分析，并通過NURBS 曲線對結構內部的參數進行插值，避免對整個結構進行離散，僅僅離散結構的邊界。②等幾何邊界元法采用NURBS 函數表示物體的邊界，使物體的邊界具有精確的幾何表示，計算精度更高。③采用增廣乘子法對優化模型進行轉換，收斂性優于罰函數法。④梯度優化算法的收斂速度較快，但智能優化算法的全局尋優能力更強，可得到性能更好的結構。⑤LSEDA 是一種融合Gauss 分布和Lévy 分布的混合模型[22]，不僅可以保證樣本集中于最優解附近，而且樣本的分布范圍更廣，全局尋優能力更強。⑥采用精英反向學習策略和LSEDA 改進花朵授粉算法，提高算法的收斂速度和魯棒性。Ackley 函數的測試結果表明，本文改進的花朵授粉算法尋優能力更強，收斂更快。2 個形狀優化算例表明，NURBS 曲線重構的結構邊界靈活多樣，且僅需要對結構邊界離散即可。

1 等幾何邊界元法

1.1 NURBS 曲線

階數相同的情況下，B 樣條曲線的基函數數量多于貝塞爾（Bézier）樣條曲線，計算精度更高[23]。若采用p階的B 樣條曲線，則B 樣條曲線的基函數為[24-28]

式中：ξ為節點向量中的一串序列；in為基函數編號。

B 樣條曲線上各節點的坐標為

式中：Px(ξ)、Py(ξ)分別為控制點的橫、縱坐標。

NURBS 曲線是一種具有非均勻節點向量 Ξ的有理B 樣條曲線，且有理B 樣條曲線在B 樣條曲線的各節點上添加了權重項 ωin。NURBS 曲線的節點向量和節點坐標分別為

階數p=1,2,3分別對應線性、二次、三次B 樣條曲線，而且B 樣條曲線的p?m階導數連續[23]。

1.2 邊界元法

當 源 點x承 受 沿方向i的 載荷ei(x)時 ，場 點x′沿方向j的位移為uj(x′)。根據邊界元法的理論，載荷ei(x)與 位移uj(x′)之間的關系式為[29-30]

根據式（13）可求得von Mises 應力[33]為

2 優化算法理論

2.1 增廣乘子法

若數學模型的目標函數F(X)有最小值，而且有b個等式約束h(X)，則該有約束優化模型可以通過增廣乘子法轉換為無約束優化模型[12]，如下：

式中：r0為 罰因子；λ′為拉格朗日乘子；a為b個等式約束條件h(X)的序號。

2.2 花朵授粉算法

Yang[34]提出了花朵授粉算法。該算法中包括4 種花朵授粉規則：①全局授粉（生物授粉和異花授粉）；②局部授粉（非生物授粉和自花授粉）；③昆蟲授粉；④局部授粉與全局授粉。

規則①和規則③中，采用Lévy 飛行的全局授粉表達式為

Mantegna 算法是一種效率最高的生成服從Lévy分布的偽隨機數的算法。該算法生成的步長為

規則④中，當隨機數大于指定概率時，采用全局授粉策略；當隨機數小于指定概率時，采用局部授粉策略。

為了提高算法的全局搜索能力，本文采用了精英反向學習策略，如下：

式中：Nv為目標函數下F(X)中自變量X的個數；N(0,1)為服從標準正態分布的隨機數。

由圖1 中的概率密度曲線可知，由于普通EDA通過Gauss 分布生成新解，且新解完全分布于最優解附近，分布范圍太小，導致算法的效率較低。然而，LSEDA 采用了Gauss-Lévy 混合模型，新解主要分布于最優解附近，而且分布范圍更廣，因此LSEDA的效率高于EDA。其中，Cauchy 分布又稱柯西分布，是一種特殊的Lévy 分布。

圖1 對數坐標系下的Gauss-Lévy 混合模型Fig.1 Gauss-Lévy mixed model in logarithmic coordinate

本文改進后的花朵授粉算法的計算步驟如圖2所示，具體步驟如下：

圖2 改進的花朵授粉算法流程Fig.2 Flow chart of improved flower pollination algorithm

步驟1初始化種群。

步驟2計算個體對應的適應度，并按升序對個體進行排列。

步驟3當前的最優解采用精英反向學習策略更新。

步驟4最差的25%個體按照LSEDA 更新。通過Gauss-Lévy 混合模型生成的隨機數，以及當前25%的最優個體的均值和標準差生成新解，并替換25%的最差個體。

步驟5其他未更新的個體采用基本花朵授粉算法更新。隨機數大于0.5 時，全局搜索；隨機數小于0.5 時，局部搜索。

步驟6達到迭代上限時，算法終止。否則，算法轉步驟2。

現采用Ackley 函數測試本文算法的改進效果。該函數中，變量的取值范圍為(?32.768,32.768)，且函數的最小值0 所在位置為(0, 0, ···, 0)。通常，Ackley 函數的各參數為：a0=20，b0=0.2，c0=2π，d0為 變量X的維度，且Ackley 函數的表達式為

設花朵授粉算法的種群數量為20 個，迭代上限為200 步。同時采用精英反向學習策略和LSEDA 這2 種算法的花朵授粉算法與基本花朵授粉算法的迭代過程如圖3 所示。

圖3 兩種花朵授粉算法的對比Fig.3 Comparation of two kinds of flower pollination algorithms

由圖3 可知，本文算法在第14 步收斂，最優解為（?0.294 8, ?0.100 1），最小值為8.881 8×10?16；而基本花朵授粉算法在136 步收斂，最優解為（?0.447 0,?0.777 7），最小值為0.001 4。因此，本文算法尋優能力更強，且將采用本文算法優化等幾何邊界元模型，以實現結構的形狀優化設計。

3 基于等幾何邊界元法的形狀優化

本文通過精英反向學習策略和LSEDA 這2 種算法改進了花朵授粉算法，并將該算法用于優化等幾何邊界元模型中的控制點，提出一種新的形狀優化設計算法。該算法的計算流程如圖4 所示。

圖4 本文算法流程Fig.4 Flowchart of the proposed algorithm

3.1 位移最小化的形狀優化

結構左端面上方3 個控制點固定，右端面下方3 個控制點承受載荷，橫向載荷50 N，縱向載荷200 N。彈性模量為 2×105MPa，泊松比為0.3。初始結構的參數設置如表1 所示，結構示意圖如圖5(a)所示，初始形狀如圖5(b)所示，NURBS 基函數如圖5(c)所示。

表1 位移算例的控制點及節點向量Table 1 Control points and knot vectors for displacement example

圖5 位移最小化Fig.5 Minimization of displacement

由式（4）可知，表1 中 Ξ對應的NURBS 曲線階數p為2 階，[0,1]內包含的節點數m為10，因此基函數個數為10+2+1=13 個。式（5）為NURBS 基函數的計算式。

現在需要將結構的面積減少至60%，并且最小化每個控制點的最大位移。設罰因子為 1×102，拉格朗日乘子的初始值為0，且每迭代一次，拉格朗日乘子增加 1×10?6，花朵授粉算法所得結果如圖5(d)所示，迭代過程如圖5(e)所示。由圖5(e)可知，花朵授粉算法在10 步以內收斂，收斂性非常好。70 次迭代后，由圖5(c)中的NURBS 基函數構建的結構如圖5(d)所示，且結構中的各控制點的最大位移下降至0.297 2 mm，體積分數為0.602 2。

3.2 應力最小化的形狀優化

結構上端面承受縱向載荷2 000 N，下端面兩端固定，彈性模量為 2×105MPa，泊松比為0.3。初始結構的參數設置見表2，結構示意圖見圖6(a)，初始結構見圖6(b)，NURBS 基函數如圖6(c)所示。現需要將結構的面積減少，同時最小化每個控制點的最大應力。設罰因子為 1×104，拉格朗日乘子的初始值為0，且每迭代一次，拉格朗日乘子增加 1×10?5。由圖6(e)可知，花朵授粉算法在20 次迭代以內，曲線趨于平緩，收斂性較好。由圖6(c)構建的結果見圖6(d)。

表2 應力算力的控制點及節點向量Table 2 Control points and knot vectors for stress example

圖6 應力最小化Fig.6 Minimization of stress

4 結 論

通過精英反向學習策略及LSEDA 對花朵授粉算法進行改進，并用于優化等幾何邊界元模型的形狀，算例結果表明：

1）通過精英反向學習策略及LSEDA 這2 種算法改進了花朵授粉算法，且Ackley 函數的測試結果表明，本文算法尋優能力更強。

2）通過本文算法求得的控制點坐標重構了NURBS 曲線，而且該NURBS 曲線非常光滑，能精確地表示結構邊界。

3）結構示意圖中，“單元邊界點”均位于結構的邊界，而結構的其余節點均通過NURBS 基函數及權值進行插值求得，因此“配置點”不一定位于結構的邊界上。