基于Kriging模型和模擬退火粒子群算法的結(jié)構(gòu)有限元模型修正

康俊濤 柯志涵 胡 佳

(武漢理工大學(xué)土木工程與建筑學(xué)院 武漢 430070)

0 引 言

工程中，實(shí)際結(jié)構(gòu)現(xiàn)場(chǎng)試驗(yàn)響應(yīng)值與建立的結(jié)構(gòu)有限元分析模型計(jì)算值往往存在誤差，需要對(duì)建立的初始有限元模型進(jìn)行修正.模型修正在結(jié)構(gòu)的性能評(píng)價(jià)[1]、健康監(jiān)測(cè)[2]、損傷識(shí)別[3],以及驗(yàn)證結(jié)構(gòu)設(shè)計(jì)等領(lǐng)域中廣泛應(yīng)用.有限元模型修正，是建立的初始有限元分析模型利用結(jié)構(gòu)現(xiàn)場(chǎng)試驗(yàn)響應(yīng)值來(lái)修正模型計(jì)算值，使得修正后模型計(jì)算值與結(jié)構(gòu)試驗(yàn)值誤差最小.直接調(diào)用有限元軟件進(jìn)行模型修正存在計(jì)算效率低，不易于二次開發(fā)等缺點(diǎn)[4-5].采用代理模型結(jié)合智能算法進(jìn)行模型修正，克服了修正過(guò)程中大量調(diào)用有限元模型計(jì)算的不足.代理模型是利用顯式函數(shù)式擬合構(gòu)造出結(jié)構(gòu)有限元模型設(shè)計(jì)參數(shù)與其響應(yīng)之間復(fù)雜的隱式關(guān)系的一種數(shù)學(xué)模型，其避免了復(fù)雜的有限元計(jì)算.任偉新等[6]在模型修正中采用代理模型，提高了修正的效率.

常用的代理模型有響應(yīng)面方法[7](response surface method,RSM)、徑向基函數(shù)[8](radical basis function,RBF)、Kriging模型[9]等.其中Kriging模型由于存在隨機(jī)過(guò)程，更適合于復(fù)雜非線性的結(jié)構(gòu)輸入、輸出模型，具有非常高的精度.

模型修正問(wèn)題歸根結(jié)底可以看成優(yōu)化問(wèn)題，常結(jié)合優(yōu)化算法求解.基本粒子群算法(particle swarm optimization，PSO)的思想源于模擬鳥群的覓食行為，是由Kennedy等[10]提出的一種群體協(xié)作隨機(jī)搜索算法.其實(shí)現(xiàn)容易、概念簡(jiǎn)單、魯棒性強(qiáng)、收斂速度快，在結(jié)構(gòu)模型修正[11]、優(yōu)化設(shè)計(jì)[12]等領(lǐng)域得到廣泛應(yīng)用.但PSO在搜索過(guò)程中容易陷入局部極值解，存在早熟等缺點(diǎn)[13].而模擬退火粒子群算法[14](simulated annealing particle swarm optimization，SimuAPSO)在搜索過(guò)程中受溫度參數(shù)的控制而避免陷入局部極值解，能夠較快的收斂到全局最優(yōu)解,提高計(jì)算效率.

1 Kriging函數(shù)模型及試驗(yàn)設(shè)計(jì)

1.1 Kriging函數(shù)模型

Kriging函數(shù)模型是由回歸模型與隨機(jī)過(guò)程相加而成.其形式為

y(X)=fT(X)β+z(X)

(1)

式中：X為維度為d的設(shè)計(jì)參數(shù)[x1,x2,…,xd]；y(X)為模型響應(yīng)值；fT(X)β為回歸模型，fT(X)=[f1(X),f2(X),…,fm(X)]，β=[β1,β2,…,βm]T為回歸矩陣，常用二次無(wú)交叉項(xiàng)多項(xiàng)式模型為

(2)

式中：z(X)為均值為0協(xié)方差為σ2的隨機(jī)過(guò)程，提供模擬的局部近似.z(X)中任意兩個(gè)參數(shù)的協(xié)方差表示為

Cov(z(X)z(ω))=σ2R(X,ω)

(3)

式中：R(X,ω)為任意兩個(gè)樣本點(diǎn)X和ω之間的變異函數(shù),其一般形式為

(4)

式中：Rj(dj)為變異函數(shù)的核函數(shù)，通常采用高斯函數(shù),則變異的函數(shù)形式為

(5)

式中：θj為計(jì)算空間相關(guān)函數(shù)中大于0的系數(shù)；n為變異的維數(shù).假設(shè)有N個(gè)樣本點(diǎn)，響應(yīng)陣和多項(xiàng)式陣分別為Y=[y(X1)…y(XN)]，F(xiàn)=[fT(X1)…fT(XN)]T，得到β和σ2值為

(6)

(7)

式中：空間相關(guān)矩陣為

(8)

(9)

式中：|R|為矩陣R的行列式值；θ值可以通過(guò)遺傳算法、模式搜索、粒子群算法等優(yōu)化方法得到.確定了θ1,θ2,…,θd的值，最后利用已得到的模型系數(shù)和變異函數(shù)參數(shù)，得到最好線性無(wú)偏估計(jì)值為

(10)

1.2 試驗(yàn)設(shè)計(jì)

試驗(yàn)設(shè)計(jì)是構(gòu)建Kriging函數(shù)的基礎(chǔ)，對(duì)修正參數(shù)進(jìn)行試驗(yàn)設(shè)計(jì)抽樣得到樣本點(diǎn)作為輸入集，并將樣本點(diǎn)代入有限元計(jì)算得到結(jié)構(gòu)響應(yīng)作為輸出集，根據(jù)輸入和輸出集擬合得到Kriging函數(shù).因此試驗(yàn)設(shè)計(jì)方法直接影響構(gòu)造Kriging模型的精度，同時(shí)也是決定結(jié)構(gòu)有限元計(jì)算量的主要因素.

本文樣本抽樣方法選用的是由Timothy 等[15]提出的一種基于分層抽樣思想的拉丁超立方設(shè)計(jì)(LHD).LHD將變量空間分成n個(gè)子區(qū)間，每個(gè)子區(qū)間概率相同且互不重疊,然后分別獨(dú)立隨機(jī)地在每個(gè)子區(qū)間內(nèi)采樣一次，以保證樣本點(diǎn)的均勻性.拉丁超立方抽樣的樣本為

(11)

式中：n為樣本點(diǎn)數(shù)；i為水平數(shù)；j為維數(shù)；dij為1到n的獨(dú)立隨機(jī)排列數(shù)；eij為[0,1]之間的隨機(jī)數(shù).

2 模擬退火粒子群算法

模擬退火粒子群算法是一種源于固體退火原理的啟發(fā)式搜索算法，其步驟如下.

步驟1隨機(jī)初始化各粒子的位置與速度.

步驟2評(píng)價(jià)每個(gè)微粒的適應(yīng)度，將第i個(gè)微粒所經(jīng)歷的歷史最好點(diǎn)儲(chǔ)存在pi中，將所有pi中的最優(yōu)值儲(chǔ)存在pg中.

(12)

(13)

vi(t+1)=φ{(diào)vi(t)+c1r1[pi-xi(t)]+

c2r2[pg-xi(t)]}

(14)

xi(t+1)=xi(t)+vi(t+1)

(15)

(16)

式中：vi(t)為粒子溫度t時(shí)的速度；xi(t)為粒子溫度t時(shí)的位置；ω為慣性權(quán)重；pi為每個(gè)粒子經(jīng)過(guò)的歷史最優(yōu)點(diǎn)；pg為所有粒子經(jīng)過(guò)的全局最優(yōu)點(diǎn)；c1、c2為學(xué)習(xí)因子；r1,r2為擾動(dòng)因子，r1,r2∈[0,1].

步驟6計(jì)算各微粒新的目標(biāo)值,更新各微粒的pi值及群體的pg值.

步驟7進(jìn)行退溫操作，tk+1=λtk，其中0<λ<1.

步驟8若滿足運(yùn)算精度或迭代次數(shù)等停止條件時(shí)，輸出當(dāng)前最優(yōu)值，搜索結(jié)束，否則返回步驟4繼續(xù)搜索.

3 桁架數(shù)值模擬算例

通過(guò)一空間桁架進(jìn)行數(shù)值模擬，對(duì)設(shè)計(jì)參數(shù)人為給定攝動(dòng)量，構(gòu)造Kriging函數(shù)模型代替空間桁架Midas結(jié)構(gòu)，結(jié)合SimuAPSO進(jìn)行尋優(yōu)，完成對(duì)桁架結(jié)構(gòu)的修正.

具體修正流程圖見圖1.

圖1 模型修正流程圖

3.1 桁架模型

利用一個(gè)空間桁架結(jié)構(gòu)作為算例驗(yàn)證構(gòu)建Kriging法用于模型修正的可行性，見圖2.桁架豎直和水平方向每根桿單元長(zhǎng)度均為0.5 m，總長(zhǎng)度為2 m、寬度為0.5 m、高度0.5 m，彈性模量為2.06×105MPa，采用管形截面,外徑為16 mm，管壁厚2.5 mm，截面面積為7.85×10-4m2.桁架材料的初始密度(用ρ表示)為7 852 kg/m3,用Midas建立空間桁架模型見圖3.

圖2 空間桁架結(jié)構(gòu)

圖3 Midas建立的空間桁架模型

桁架由34個(gè)桿件單元和16個(gè)節(jié)點(diǎn)組成，節(jié)點(diǎn)1約束住三個(gè)平動(dòng)自由度，節(jié)點(diǎn)5，6約束住z方向平動(dòng)自由度，節(jié)點(diǎn)10約束住y，z方向平動(dòng)自由度，且該四個(gè)節(jié)點(diǎn)均約束住繞x，z軸方向的轉(zhuǎn)動(dòng)自由度.對(duì)34根桿件的密度參數(shù)進(jìn)行擾動(dòng)，代入有限元模型計(jì)算其響應(yīng)，經(jīng)過(guò)敏感性分析，選取ρ1，ρ5，ρ9，ρ20四個(gè)敏感參數(shù)，人為地?cái)_動(dòng)這四個(gè)參數(shù)數(shù)值，ρ1增加20%至9 420 kg/m3，ρ5增加15%至9 028 kg/m3，ρ9增加30%至10 205 kg/m3，ρ20增加10%至8 635 kg/m3.未損傷模型代表初始模型，擾動(dòng)后的損傷模型代表實(shí)際模型，經(jīng)過(guò)Midas空間桁架模型計(jì)算得到實(shí)際模型與初始模型的前6階頻率見表1.

表1 初始模型和真實(shí)模型的前6階頻率

3.2 試驗(yàn)設(shè)計(jì)

采用LHD抽樣方法，對(duì)選取的ρ1，ρ5，ρ9，ρ20這四個(gè)參數(shù)進(jìn)行抽樣40次得到試驗(yàn)點(diǎn)并代入有限元模型中計(jì)算得到對(duì)應(yīng)的頻率數(shù)據(jù).設(shè)計(jì)參數(shù)的變化范圍為±50%，即抽樣點(diǎn)密度下界為3 925 kg/m3，上界為11 775 kg/m3.將ρ1，ρ9抽樣的樣本散點(diǎn)圖繪于圖4，試驗(yàn)點(diǎn)均勻地散布于輸入?yún)?shù)空間中.

圖4 利用LHD得到的ρ1，ρ9樣本散點(diǎn)圖

3.3 構(gòu)建Kriging模型及精度檢驗(yàn)

LHD抽樣得到樣本點(diǎn)作輸入集，利用Midas計(jì)算得到的前六階頻率作輸出集，構(gòu)建Kriging模型.選取桿件單元1和9的密度參數(shù)作為自變量，桿件單元5和20的參數(shù)值為初始值不變，做Kriging函數(shù)擬合圖，見圖5.由圖5可知,構(gòu)造的Kriging函數(shù)響應(yīng)面平滑，圖6是構(gòu)建Kriging模型的均方誤差圖，其MSE值數(shù)量級(jí)在10-4，表明擬合構(gòu)建的Kriging函數(shù)設(shè)計(jì)參數(shù)估計(jì)值與真值誤差較小，精確度較高，能代替有限元模型.

圖5 頻率f1和密度ρ1，ρ9之間的函數(shù)擬合圖

圖6 Kriging函數(shù)的均方誤差圖

3.4 構(gòu)建目標(biāo)函數(shù)

以構(gòu)建的Kriging代理模型得到的前六階頻率與實(shí)際結(jié)構(gòu)的前六階頻率相對(duì)誤差的平方和作為目標(biāo)函數(shù).

(17)

式中：fkrigi為構(gòu)建的Kriging代理模型得到的第i階頻率響應(yīng)；fri為實(shí)際損傷結(jié)構(gòu)的第i階頻率響應(yīng).

3.5 SimuAPSO優(yōu)化

SimuAPSO算法參數(shù)設(shè)置如下：初始粒子總數(shù)為100，迭代次數(shù)200次，粒子維數(shù)四維，個(gè)體認(rèn)知部分權(quán)重與社會(huì)認(rèn)知部分權(quán)重均取2.05.優(yōu)化結(jié)果及誤差見圖7.由圖7可知，四個(gè)參數(shù)的初始值為7 850 kg/m3,經(jīng)模擬退火粒子群優(yōu)化得到的修正值分別為9 807，8 521,10 147,8 901 kg/m3,損傷結(jié)構(gòu)的真實(shí)值分別為9 420,9 028，10 205,8 635 kg/m3.圖7中條形圖上數(shù)字表示修正值與真實(shí)值的相對(duì)誤差，可見最大相對(duì)誤差為5.61%，最小相對(duì)誤差為0.57%，滿足誤差要求.

圖7 模擬退火粒子群算法參數(shù)修正圖

圖 8為模擬退火粒子群算法和經(jīng)典粒子群算法對(duì)目標(biāo)函數(shù)的最優(yōu)解進(jìn)化的曲線圖.

圖8 兩種算法對(duì)目標(biāo)函數(shù)的最優(yōu)解進(jìn)化曲線圖

由圖8可知，兩種算法修正后的目標(biāo)函數(shù)值基本一致.但SimuAPSO優(yōu)化效率較高，可以避免陷入局部極值解，更快的收斂到全局最優(yōu)解.標(biāo)準(zhǔn)粒子群算法得到最優(yōu)解的迭代次數(shù)需要20次以上，而本文采取的模擬退火粒子群將迭代次數(shù)縮小到10次左右，提高了優(yōu)化效率.

4 結(jié) 論

1) 基于Kriging代理模型的結(jié)構(gòu)有限元模型修正避免了重復(fù)調(diào)用有限元進(jìn)行計(jì)算，縮短了運(yùn)算時(shí)間.

2) 基于MATLAB編程的Kriging模型修正過(guò)程易于與優(yōu)化算法相結(jié)合，采用模擬退火粒子群算法，較快收斂到全局最優(yōu)解，顯著提高優(yōu)化效率.

3) 通過(guò)空間桁架算例表明，基于Kriging代理模型，結(jié)合模擬退火粒子群算法實(shí)現(xiàn)有限元模型修正是一個(gè)行之有效的方法.在結(jié)構(gòu)優(yōu)化設(shè)計(jì)、模型修正與確認(rèn)、損傷識(shí)別等工程領(lǐng)域有很好的應(yīng)用價(jià)值.