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

(1.三峽大學(xué) 水電機(jī)械設(shè)備設(shè)計(jì)與維護(hù)湖北省重點(diǎn)實(shí)驗(yàn)室，宜昌 443002；2.三峽大學(xué) 機(jī)械與動(dòng)力學(xué)院，宜昌 443002)

固體力學(xué)中常用的離散方法有有限差分法、有限元法、邊界元法、有限體積法、基面力元法[1]，巖土力學(xué)中常用的離散方法有離散元法和非連續(xù)變形分析，流體力學(xué)中常用的離散方法有光滑粒子法、移動(dòng)粒子半隱式法，斷裂力學(xué)中常用的離散方法有擴(kuò)展有限元法。在固體力學(xué)的形狀優(yōu)化問(wèn)題中，邊界元法僅需要對(duì)邊界進(jìn)行離散即可，比其他方法更適用于形狀優(yōu)化問(wèn)題。

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

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

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

1 等幾何邊界元法

1.1 NURBS 曲線

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

式中：ξ為節(jié)點(diǎn)向量中的一串序列；in為基函數(shù)編號(hào)。

B 樣條曲線上各節(jié)點(diǎn)的坐標(biāo)為

式中：Px(ξ)、Py(ξ)分別為控制點(diǎn)的橫、縱坐標(biāo)。

NURBS 曲線是一種具有非均勻節(jié)點(diǎn)向量 Ξ的有理B 樣條曲線，且有理B 樣條曲線在B 樣條曲線的各節(jié)點(diǎn)上添加了權(quán)重項(xiàng) ωin。NURBS 曲線的節(jié)點(diǎn)向量和節(jié)點(diǎn)坐標(biāo)分別為

階數(shù)p=1,2,3分別對(duì)應(yīng)線性、二次、三次B 樣條曲線，而且B 樣條曲線的p?m階導(dǎo)數(shù)連續(xù)[23]。

1.2 邊界元法

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

根據(jù)式（13）可求得von Mises 應(yīng)力[33]為

2 優(yōu)化算法理論

2.1 增廣乘子法

若數(shù)學(xué)模型的目標(biāo)函數(shù)F(X)有最小值，而且有b個(gè)等式約束h(X)，則該有約束優(yōu)化模型可以通過(guò)增廣乘子法轉(zhuǎn)換為無(wú)約束優(yōu)化模型[12]，如下：

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

2.2 花朵授粉算法

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

規(guī)則①和規(guī)則③中，采用Lévy 飛行的全局授粉表達(dá)式為

Mantegna 算法是一種效率最高的生成服從Lévy分布的偽隨機(jī)數(shù)的算法。該算法生成的步長(zhǎng)為

規(guī)則④中，當(dāng)隨機(jī)數(shù)大于指定概率時(shí)，采用全局授粉策略；當(dāng)隨機(jī)數(shù)小于指定概率時(shí)，采用局部授粉策略。

為了提高算法的全局搜索能力，本文采用了精英反向?qū)W習(xí)策略，如下：

式中：Nv為目標(biāo)函數(shù)下F(X)中自變量X的個(gè)數(shù)；N(0,1)為服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)。

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

圖1 對(duì)數(shù)坐標(biāo)系下的Gauss-Lévy 混合模型Fig.1 Gauss-Lévy mixed model in logarithmic coordinate

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

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

步驟1初始化種群。

步驟2計(jì)算個(gè)體對(duì)應(yīng)的適應(yīng)度，并按升序?qū)€(gè)體進(jìn)行排列。

步驟3當(dāng)前的最優(yōu)解采用精英反向?qū)W習(xí)策略更新。

步驟4最差的25%個(gè)體按照LSEDA 更新。通過(guò)Gauss-Lévy 混合模型生成的隨機(jī)數(shù)，以及當(dāng)前25%的最優(yōu)個(gè)體的均值和標(biāo)準(zhǔn)差生成新解，并替換25%的最差個(gè)體。

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

步驟6達(dá)到迭代上限時(shí)，算法終止。否則，算法轉(zhuǎn)步驟2。

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

設(shè)花朵授粉算法的種群數(shù)量為20 個(gè)，迭代上限為200 步。同時(shí)采用精英反向?qū)W習(xí)策略和LSEDA 這2 種算法的花朵授粉算法與基本花朵授粉算法的迭代過(guò)程如圖3 所示。

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

由圖3 可知，本文算法在第14 步收斂，最優(yōu)解為（?0.294 8, ?0.100 1），最小值為8.881 8×10?16；而基本花朵授粉算法在136 步收斂，最優(yōu)解為（?0.447 0,?0.777 7），最小值為0.001 4。因此，本文算法尋優(yōu)能力更強(qiáng)，且將采用本文算法優(yōu)化等幾何邊界元模型，以實(shí)現(xiàn)結(jié)構(gòu)的形狀優(yōu)化設(shè)計(jì)。

3 基于等幾何邊界元法的形狀優(yōu)化

本文通過(guò)精英反向?qū)W習(xí)策略和LSEDA 這2 種算法改進(jìn)了花朵授粉算法，并將該算法用于優(yōu)化等幾何邊界元模型中的控制點(diǎn)，提出一種新的形狀優(yōu)化設(shè)計(jì)算法。該算法的計(jì)算流程如圖4 所示。

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

3.1 位移最小化的形狀優(yōu)化

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

表1 位移算例的控制點(diǎn)及節(jié)點(diǎn)向量Table 1 Control points and knot vectors for displacement example

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

由式（4）可知，表1 中 Ξ對(duì)應(yīng)的NURBS 曲線階數(shù)p為2 階，[0,1]內(nèi)包含的節(jié)點(diǎn)數(shù)m為10，因此基函數(shù)個(gè)數(shù)為10+2+1=13 個(gè)。式（5）為NURBS 基函數(shù)的計(jì)算式。

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

3.2 應(yīng)力最小化的形狀優(yōu)化

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

表2 應(yīng)力算力的控制點(diǎn)及節(jié)點(diǎn)向量Table 2 Control points and knot vectors for stress example

圖6 應(yīng)力最小化Fig.6 Minimization of stress

4 結(jié) 論

通過(guò)精英反向?qū)W習(xí)策略及LSEDA 對(duì)花朵授粉算法進(jìn)行改進(jìn)，并用于優(yōu)化等幾何邊界元模型的形狀，算例結(jié)果表明：

1）通過(guò)精英反向?qū)W習(xí)策略及LSEDA 這2 種算法改進(jìn)了花朵授粉算法，且Ackley 函數(shù)的測(cè)試結(jié)果表明，本文算法尋優(yōu)能力更強(qiáng)。

2）通過(guò)本文算法求得的控制點(diǎn)坐標(biāo)重構(gòu)了NURBS 曲線，而且該NURBS 曲線非常光滑，能精確地表示結(jié)構(gòu)邊界。

3）結(jié)構(gòu)示意圖中，“單元邊界點(diǎn)”均位于結(jié)構(gòu)的邊界，而結(jié)構(gòu)的其余節(jié)點(diǎn)均通過(guò)NURBS 基函數(shù)及權(quán)值進(jìn)行插值求得，因此“配置點(diǎn)”不一定位于結(jié)構(gòu)的邊界上。