基于混合算法的自由曲面網(wǎng)格結(jié)構(gòu)多目標(biāo)優(yōu)化

王志成,曹正罡,趙 林,李展熇,范 峰,孫 瑛

(1.結(jié)構(gòu)工程災(zāi)變與控制教育部重點(diǎn)實(shí)驗(yàn)室(哈爾濱工業(yè)大學(xué)),哈爾濱150090；2.土木工程智能防災(zāi)減災(zāi)工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室(哈爾濱工業(yè)大學(xué)),哈爾濱150090；3.中南大學(xué) 土木工程學(xué)院,長沙410083)

自由曲面結(jié)構(gòu)由于其卓越的建筑表現(xiàn)力得到了廣泛應(yīng)用(圖1)，成為空間結(jié)構(gòu)發(fā)展的主要趨勢之一。自由曲面是無法用單個或幾個解析函數(shù)表達(dá)的曲面，是一種明顯區(qū)別于傳統(tǒng)建筑造型的曲面。自由曲面結(jié)構(gòu)早期是通過物理實(shí)驗(yàn)的方法形成合理形態(tài)[1]。近年來，學(xué)者們陸續(xù)提出了Bezier曲面[2]、B樣條曲面[3]和非均勻有理B樣條曲面(NURBS)[4]等以計(jì)算圖形學(xué)為基礎(chǔ)的形態(tài)創(chuàng)建方法。目前，NURBS曲線與曲面已成為結(jié)構(gòu)設(shè)計(jì)領(lǐng)域的核心技術(shù)[5]。

圖1 阿布扎比Yas酒店

結(jié)構(gòu)優(yōu)化是使結(jié)構(gòu)具有良好受力性能的有效方式，包括形態(tài)優(yōu)化、拓?fù)鋬?yōu)化等。因此，學(xué)者們基于自由曲面數(shù)學(xué)建構(gòu)的方法，通過優(yōu)化來得到合理的自由曲面形態(tài)。Tomas等[6]以曲面控制點(diǎn)高度為優(yōu)化變量，分別以結(jié)構(gòu)的最大位移、單元應(yīng)力和結(jié)構(gòu)整體質(zhì)量為優(yōu)化目標(biāo)進(jìn)行優(yōu)化。Le等[7]提出了基于梯度法來實(shí)現(xiàn)自由曲面的形態(tài)優(yōu)化。崔昌禹等[8]以結(jié)構(gòu)應(yīng)變能為優(yōu)化目標(biāo)，基于應(yīng)變能敏感度調(diào)整曲面節(jié)點(diǎn)高度優(yōu)化曲面，并在工程中得以應(yīng)用。Ohmori等[9]以結(jié)構(gòu)應(yīng)變能及曲面優(yōu)化前后的相似程度為優(yōu)化目標(biāo)，以曲面控制點(diǎn)坐標(biāo)和殼體厚度為優(yōu)化變量進(jìn)行曲面形態(tài)優(yōu)化，但實(shí)際得到的優(yōu)化曲面與初始曲面形狀的差別依然較大。

為使自由曲面結(jié)構(gòu)既有較好的力學(xué)性能又滿足建筑要求，Winslow等[10]基于傳統(tǒng)遺傳算法，通過調(diào)整曲面上兩個主要方向網(wǎng)格線之間的角度進(jìn)行自由曲面結(jié)構(gòu)多目標(biāo)優(yōu)化。馮若強(qiáng)等[11]以結(jié)構(gòu)應(yīng)變能和桿長均方差為優(yōu)化目標(biāo)，曲面控制點(diǎn)高度為優(yōu)化變量，采用加權(quán)法進(jìn)行自由曲面索撐網(wǎng)殼的多目標(biāo)優(yōu)化。目前關(guān)于自由曲面網(wǎng)格結(jié)構(gòu)的多目標(biāo)優(yōu)化均是基于加權(quán)法將多目標(biāo)問題轉(zhuǎn)化為單目標(biāo)問題或采用傳統(tǒng)遺傳算法，且優(yōu)化目標(biāo)并未綜合考慮自由曲面網(wǎng)格的相似性、流暢性以及規(guī)整性。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，元啟發(fā)式算法被廣泛應(yīng)用于工程領(lǐng)域，蘇巖等[12]在自由曲面網(wǎng)格結(jié)構(gòu)形態(tài)創(chuàng)建過程中采用了群智能算法。然而在多目標(biāo)優(yōu)化問題中，現(xiàn)有多目標(biāo)優(yōu)化算法的效率不高，因此需要提高算法的尋優(yōu)效率。

基于上述研究基礎(chǔ)，本文分別基于NURBS曲線和曲面技術(shù)建立自由曲面索撐網(wǎng)殼、自由曲面空間網(wǎng)格結(jié)構(gòu)，以自由曲面控制點(diǎn)高度為優(yōu)化變量，以結(jié)構(gòu)應(yīng)變能、考慮曲面網(wǎng)格相似性、流暢性以及規(guī)整性的幾何綜合量化指標(biāo)為目標(biāo)函數(shù)。為提高優(yōu)化效率及優(yōu)化結(jié)果精確性，結(jié)合目標(biāo)函數(shù)關(guān)于優(yōu)化變量的敏感度對NSGA-II算法進(jìn)行改進(jìn)，提出敏感度混合進(jìn)化算法(SH-NSGA-II)，采用該算法對自由曲面網(wǎng)格結(jié)構(gòu)進(jìn)行多目標(biāo)優(yōu)化。將優(yōu)化結(jié)果與基于NSGA-II[13]、MOEA/D[14]以及SPEA2[15]算法的優(yōu)化結(jié)果比較，驗(yàn)證敏感度混合進(jìn)化算法的合理性。

1 自由曲面結(jié)構(gòu)幾何建模方法

盡管諸多表示曲面的數(shù)學(xué)模型有各自的優(yōu)缺點(diǎn)，但NURBS在工程應(yīng)用中是最成功的。NURBS能夠比傳統(tǒng)的建模方式更好地控制物體表面的曲線度，從而能夠創(chuàng)建出更逼真、生動的造型。因此，本文基于NURBS技術(shù)生成自由曲面。

自由曲面索撐網(wǎng)殼的網(wǎng)格為平面四邊形，每個單元平面內(nèi)布置對角交叉索[16]。平面四邊形網(wǎng)格由表面平移法形成[17]，即首先基于NURBS曲線構(gòu)造一條準(zhǔn)線和一條母線，并將曲線按等弦長劃分；將準(zhǔn)線沿母線平行移動，再將母線和準(zhǔn)線調(diào)換并再次平移，形成均勻的平面四邊形網(wǎng)格(圖2)。自由曲面空間網(wǎng)格結(jié)構(gòu)采用四邊形網(wǎng)格，基于NURBS曲面技術(shù)構(gòu)造自由曲面，采用等參線分割法[18]將曲面劃分為多個四邊形網(wǎng)格。

圖2 表面平移法幾何原理

2 自由曲面結(jié)構(gòu)優(yōu)化模型

控制點(diǎn)是對曲面形狀影響較大的參數(shù)，當(dāng)同時調(diào)整控制點(diǎn)的三向坐標(biāo)時，易引起結(jié)構(gòu)網(wǎng)格畸變、曲面扭曲。因此選取控制點(diǎn)高度作為優(yōu)化變量，自由曲面結(jié)構(gòu)的優(yōu)化模型表達(dá)為：

(1)

式中：F(x)為優(yōu)化目標(biāo)函數(shù)，f1(x)為自由曲面結(jié)構(gòu)靜力性能優(yōu)化目標(biāo)，f2(x)為幾何優(yōu)化目標(biāo)。

2.1 靜力性能優(yōu)化目標(biāo)

常用的結(jié)構(gòu)靜力性能評價指標(biāo)有位移、內(nèi)力和應(yīng)變能等。位移、內(nèi)力是評價結(jié)構(gòu)局部靜力性能的指標(biāo)，適用于內(nèi)力、位移最大部位明確的結(jié)構(gòu)，但自由曲面結(jié)構(gòu)的形態(tài)復(fù)雜，難以確定內(nèi)力、位移最大的部位。因此位移、內(nèi)力不宜作為優(yōu)化目標(biāo)。應(yīng)變能是能綜合反映結(jié)構(gòu)整體抵抗荷載和變形的指標(biāo)，其定義為

(2)

式中：U為曲面結(jié)構(gòu)的位移向量，P為作用在曲面結(jié)構(gòu)的節(jié)點(diǎn)荷載向量。應(yīng)變能越小，結(jié)構(gòu)抵抗荷載和抗變形的能力越高，結(jié)構(gòu)剛度越大，研究人員選擇將應(yīng)變能作為結(jié)構(gòu)靜力性能的評價指標(biāo)[19-20]。因此，本文選取應(yīng)變能作為靜力性能優(yōu)化目標(biāo)。

2.2 幾何優(yōu)化目標(biāo)

除力學(xué)性能外，為符合建筑設(shè)計(jì)要求，自由曲面結(jié)構(gòu)需要有良好的網(wǎng)格質(zhì)量，對于網(wǎng)格質(zhì)量的綜合評價包括曲面相似性、網(wǎng)格流暢性、網(wǎng)格規(guī)整性三方面。相似性指標(biāo)是用于衡量網(wǎng)格劃分前曲面與劃分后曲面的相似程度。在幾何建模過程中，需要對NURBS曲線和等參線進(jìn)行等弦長劃分，因此曲面相似性指標(biāo)D定義為

(3)

式中：L0為曲線原始長度，L為等弦長劃分后曲線長度。因此，D值越小，自由曲面在網(wǎng)格劃分前后的相似程度越好。

網(wǎng)格流暢性是評價自由曲面網(wǎng)格質(zhì)量的一個重要方面。影響曲面網(wǎng)格流暢性的主要因素是網(wǎng)格節(jié)點(diǎn)對邊的夾角以及節(jié)點(diǎn)的對角。若節(jié)點(diǎn)對邊夾角β1,3=180°，且臨邊夾角β1=β3，則節(jié)點(diǎn)具有良好的流暢性。因此，節(jié)點(diǎn)i的流暢性指標(biāo)Fi由式(4)～(6)定義[21]：

(4)

(5)

(6)

網(wǎng)格規(guī)整性主要包括桿件長度和網(wǎng)格形狀質(zhì)量兩方面。為保證自由曲面網(wǎng)格具有良好的規(guī)整性，一方面桿件長度應(yīng)盡量均勻，即桿件長度標(biāo)準(zhǔn)差δL越小越好。另一方面為自由曲面各個網(wǎng)格形狀質(zhì)量，將網(wǎng)格形狀質(zhì)量系數(shù)[17]標(biāo)準(zhǔn)差δG作為網(wǎng)格形狀規(guī)整性評價標(biāo)準(zhǔn)。因此，本文采用桿件長度標(biāo)準(zhǔn)差δL與形狀質(zhì)量系數(shù)標(biāo)準(zhǔn)差δG作為網(wǎng)格規(guī)整性指標(biāo)。

(7)

式中：n為桿件總數(shù)，Li為第i根桿件的長度，l為桿件長度均值。

綜合考慮曲面相似性、網(wǎng)格流暢性以及網(wǎng)格規(guī)整性，提出以綜合量化指標(biāo)T作為幾何優(yōu)化目標(biāo)。其定義為

T=γ1×D+γ2×F+γ3×δL+γ4×δG

(8)

式中γi為各幾何指標(biāo)在綜合量化指標(biāo)中的權(quán)重，可根據(jù)實(shí)際需求設(shè)置，但需滿足∑γi=1。

3 敏感度混合進(jìn)化算法

3.1 目標(biāo)函數(shù)敏感度

優(yōu)化變量的變化對目標(biāo)函數(shù)的影響程度稱為敏感度。敏感度是目標(biāo)函數(shù)增加最快的方向，其數(shù)值為相應(yīng)的變化率，將優(yōu)化變量沿負(fù)敏感度方向移動可使目標(biāo)函數(shù)快速減少。本文目標(biāo)函數(shù)為應(yīng)變能以及幾何綜合量化指標(biāo)，優(yōu)化變量為控制點(diǎn)的高度，目標(biāo)函數(shù)與優(yōu)化變量為非線性函數(shù)關(guān)系，用目標(biāo)函數(shù)對控制點(diǎn)高度的微分形式獲得目標(biāo)函數(shù)的敏感度，可用下式表達(dá)：

(9)

3.2 敏感度混合進(jìn)化算法流程

NSGA-II算法是基于Pareto 最優(yōu)概念處理多目標(biāo)優(yōu)化問題的優(yōu)化算法，是多目標(biāo)優(yōu)化問題中引用率最高的進(jìn)化算法。該算法基于快速非支配排序和精英策略降低了算法的復(fù)雜度，保持了種群的多樣性。因此基于NSGA-II算法，結(jié)合3.1節(jié)定義的目標(biāo)函數(shù)敏感度，提出了一種敏感度混合進(jìn)化算法，實(shí)現(xiàn)了自由曲面索撐網(wǎng)殼以及自由曲面空間網(wǎng)格結(jié)構(gòu)的多目標(biāo)優(yōu)化。基于敏感度混合進(jìn)化算法的多目標(biāo)優(yōu)化過程如下：

1)設(shè)置種群大小N，交叉概率c，變異因子r和迭代次數(shù)S。本文中參數(shù)設(shè)置為：N=50，c=0.95，r=1/m，m為優(yōu)化變量的個數(shù)。

2)根據(jù)實(shí)際需求，設(shè)置優(yōu)化變量的變化范圍。

3)初始化種群，設(shè)置控制點(diǎn)坐標(biāo)，以控制點(diǎn)高度向量x生成大小為N的初始種群P。

4)建立自由曲面索撐網(wǎng)殼結(jié)構(gòu)或自由曲面空間網(wǎng)格結(jié)構(gòu)模型，生成網(wǎng)格節(jié)點(diǎn)坐標(biāo)。

5)計(jì)算目標(biāo)函數(shù)關(guān)于優(yōu)化變量的敏感度，確定合理步長，按式(10)更新種群，式中：x為控制點(diǎn)高度向量，λ為更新步長，?fn(x)/?x為目標(biāo)函數(shù)的敏感度。

(10)

6)對更新后的種群進(jìn)行選擇、交叉和變異，產(chǎn)生子代種群Q。

7)將父代種群Pi與子代種群Q合并，得到大小為2N的種群R。

8)對種群R進(jìn)行非支配排序(F1表示Pareto序號為1的解集)，將解集F1～Fk放置種群Pi+1，直至|F1∪F2∪…∪Fk+1|≥N。

9)計(jì)算解集Fk+1的擁擠度，對解集Fk+1按擁擠度排序，按擁擠度將Fk+1存入Pi+1，直至Pi+1中解的個數(shù)為N，得到新父代種群Pi+1。

10)若滿足收斂條件則循環(huán)結(jié)束，否則返回第2步，本文以迭代次數(shù)作為收斂條件。

11)輸出解集PS，從中選取合適的解作為最終控制點(diǎn)高度。

4 敏感度混合進(jìn)化算法

本節(jié)提供了兩個基于上述優(yōu)化算法的自由曲面網(wǎng)格結(jié)構(gòu)的多目標(biāo)優(yōu)化算例，將優(yōu)化結(jié)果分別與基于NSGA-II、MOEA/D以及SPEA2算法的優(yōu)化結(jié)果比較，證明提出的算法提高了計(jì)算效率和解的精確性。算法SH-NSGA-II中參數(shù)設(shè)置為：N=50，c=0.95，r=1/m，m為優(yōu)化變量的個數(shù)。另外3種算法中的所有參數(shù)都是根據(jù)文獻(xiàn)[13-15]中的推薦值選取的，且種群大小均為50。

4.1 自由曲面索撐網(wǎng)殼多目標(biāo)優(yōu)化

現(xiàn)有一自由曲面索撐網(wǎng)殼結(jié)構(gòu)，跨度和長度均為35 m。鋼桿件采用100 mm×10 mm的圓鋼管，鋼管的彈性模量為2.06×1011N/m2，索的彈性模量和截面面積分別為1.3×1011N/m2、65 mm2，初始預(yù)應(yīng)力為100 MPa。節(jié)點(diǎn)形式為剛接節(jié)點(diǎn)，支座處剛接，4個角點(diǎn)支承，自由曲面索撐網(wǎng)殼的初始形狀見圖3。荷載形式為均布恒荷載(512 N/m2)和滿跨分布活荷載(500 N/m2)，控制點(diǎn)高度范圍為[-20+x, 20+x]。

圖3 自由曲面索撐網(wǎng)殼初始形狀

分別采用SH-NSGA-II、NSGA-II、SPEA2以及MOEA/D算法對上述結(jié)構(gòu)進(jìn)行多目標(biāo)優(yōu)化，為了比較4種算法的優(yōu)化效果，圖4給出了迭代步數(shù)為150時4種算法的Pareto解集，圖5給出了4種算法的Pareto最優(yōu)解集。迭代步數(shù)為150時，SH-NSGA-II的Pareto解集的范圍以及均勻性明顯優(yōu)于其他3種算法的Pareto解集。對比分析4種算法在迭代步數(shù)為150和400時各自的Pareto解集，SH-NSGA-II和MOEA/D算法在150步的Pareto解集與400步的Pareto解集幾乎重合，說明當(dāng)?shù)綌?shù)為150時兩種算法收斂；而算法NSGA-II、SPEA2是在迭代步數(shù)為400時收斂。考慮計(jì)算時間，收斂時其他3種算法的計(jì)算時間分別是SH-NSGA-II算法的2.893、2.382、2.783倍。

圖4 各算法的Pareto解集(步數(shù)=150)

圖5 自由曲面索撐網(wǎng)殼的最優(yōu)Pareto解集

將SH-NSGA-II算法得到的Pareto最優(yōu)解集與其他3種算法得到的最優(yōu)解集相比，最優(yōu)解集的均勻性優(yōu)于其他3種算法，多樣性明顯優(yōu)于算法SPEA2、MOEA/D，精確性優(yōu)于MOEA/D。綜上所述，SH-NSGA-II算法既可以獲得精確性較高、多樣性以及均勻性較好的Pareto最優(yōu)解集，又具有較高的計(jì)算效率。

從Pareto最優(yōu)解集中選取合適的解，保證應(yīng)變能和幾何綜合量化指標(biāo)均有所下降。表1給出了優(yōu)化后結(jié)構(gòu)各個指標(biāo)的數(shù)值，圖6為優(yōu)化前后結(jié)構(gòu)的形狀對比。由表1可知，優(yōu)化后結(jié)構(gòu)的應(yīng)變能降低了21.2%，結(jié)構(gòu)的最大位移下降了52.2%，即有效提高了結(jié)構(gòu)的力學(xué)性能；幾何綜合量化指標(biāo)降低了15.4%，且相似性指標(biāo)、索單元長度標(biāo)準(zhǔn)差、網(wǎng)格質(zhì)量系數(shù)標(biāo)準(zhǔn)差以及流暢性指標(biāo)分別降低了24.9%、15.4%、22.2%、7.5%，從而證明了將幾何綜合量化指標(biāo)作為目標(biāo)函數(shù)可有效降低結(jié)構(gòu)的各項(xiàng)幾何指標(biāo)。

表1 自由曲面索撐網(wǎng)殼多目標(biāo)優(yōu)化結(jié)果

圖6 結(jié)構(gòu)優(yōu)化前后形狀比較

4.2 自由曲面空間網(wǎng)格結(jié)構(gòu)多目標(biāo)優(yōu)化

現(xiàn)有一自由曲面空間網(wǎng)格結(jié)構(gòu)，跨度和長度均為25 m。鋼桿件采用100 mm×10 mm的圓鋼管，鋼管的彈性模量為2.06×1011N/m2。節(jié)點(diǎn)形式為剛接節(jié)點(diǎn)，支座處剛接，4個角點(diǎn)支承，自由曲面空間網(wǎng)格結(jié)構(gòu)的初始形狀見圖7。荷載形式為均布恒荷載(512 N/m2)和滿跨分布活荷載(500 N/m2)。非邊界控制點(diǎn)高度變化范圍為[-5+x, 5+x]，邊界控制點(diǎn)高度變化范圍為[-0.5+x, 0.5+x]。

圖7 自由曲面空間網(wǎng)格結(jié)構(gòu)初始形狀

采用SH-NSGA-II、NSGA-II、SPEA2以及MOEA/D算法對自由曲面空間網(wǎng)格結(jié)構(gòu)進(jìn)行多目標(biāo)優(yōu)化，圖8給出了4種算法得到的Pareto解集隨迭代步數(shù)的變化。SH-NSGA-II和NSGA-II算法在300步的Pareto解集與400步的Pareto解集相同，說明當(dāng)?shù)綌?shù)為300時算法收斂；SPEA2、MOEA/D算法在迭代步數(shù)為750、800、850時的Pareto解集幾乎重合，說明在迭代步數(shù)為750時收斂。考慮計(jì)算時間(表2)，收斂時其他3種算法的計(jì)算時間分別是SH-NSGA-II的0.953、3.066、6.091倍。

圖8 各算法Pareto解集隨迭代步數(shù)的變化

圖9為4種算法得到的Pareto最優(yōu)解集，從4種算法的Pareto最優(yōu)解集中選取合適的解，表2給出了優(yōu)化后結(jié)構(gòu)各個指標(biāo)的數(shù)值。SH-NSGA-II算法得到Pareto最優(yōu)解集的精確性優(yōu)于其他3種算法，其結(jié)構(gòu)應(yīng)變能的下降幅度大于其他3種算法，幾何指標(biāo)的下降幅度明顯大于NSGA-II、MOEA/D算法，但略小于SPEA2。因此，SH-NSGA-II算法明顯提高了計(jì)算效率和解集的精確性。

圖9 自由曲面空間網(wǎng)格結(jié)構(gòu)的Pareto最優(yōu)解集

圖10為基于SH-NSGA-II算法得到的自由曲面空間網(wǎng)格結(jié)構(gòu)優(yōu)化后的形狀，且由表2可知，基于SH-NSGA-II算法得到的結(jié)構(gòu)最大位移小于其他3種算法得到的結(jié)構(gòu)最大位移。同時各項(xiàng)幾何指標(biāo)分別降低了63.5%、18.6%、19.3%、23%，驗(yàn)證了將幾何綜合量化指標(biāo)作為目標(biāo)函數(shù)的合理性。

表2 自由曲面網(wǎng)格結(jié)構(gòu)多目標(biāo)優(yōu)化結(jié)果

圖10 結(jié)構(gòu)優(yōu)化后形狀

5 結(jié) 論

為了生成既有良好的力學(xué)性能，又滿足建筑要求的自由曲面網(wǎng)格結(jié)構(gòu)，對自由曲面網(wǎng)格結(jié)構(gòu)的多目標(biāo)優(yōu)化方法進(jìn)行改進(jìn)，提高優(yōu)化方法的效率和優(yōu)化后結(jié)構(gòu)的性能。

基于NURBS技術(shù)、并利用表面平移法和等參線法分別建立自由曲面索撐網(wǎng)殼和自由曲面空間網(wǎng)格結(jié)構(gòu)模型。選取曲面控制點(diǎn)的高度作為優(yōu)化變量，同時為充分考慮曲面相似性、網(wǎng)格流暢性以及網(wǎng)格規(guī)整性，提出以幾何綜合量化指標(biāo)T作為幾何優(yōu)化目標(biāo)，選取結(jié)構(gòu)應(yīng)變能為靜力性能優(yōu)化目標(biāo)，從而建立優(yōu)化模型。

提出敏感度混合進(jìn)化算法，基于該算法及上述優(yōu)化模型進(jìn)行多目標(biāo)優(yōu)化。通過與NSGA-II、SPEA2以及MOEA/D算法對比驗(yàn)證了所提出的算法不僅可以獲得精確性、均勻性更好的Pareto最優(yōu)解集，而且具有較高的計(jì)算效率。

兩種結(jié)構(gòu)的優(yōu)化結(jié)果表明：通過優(yōu)化結(jié)構(gòu)的形狀使結(jié)構(gòu)應(yīng)變能下降，結(jié)構(gòu)自身的力學(xué)性能有所提高。以綜合量化指標(biāo)為目標(biāo)函數(shù)可有效提高自由曲面的相似性、流暢性以及網(wǎng)格規(guī)整性，具有較高的工程實(shí)用價值。