基于Kriging模型及代理優(yōu)化算法的老舊橋梁有限元模型修正

楊 鵬

(中國鐵建昆侖投資集團(tuán)有限公司 四川成都 610000)

1 引言

當(dāng)前橋梁有限元模型的建立多以橋梁設(shè)計(jì)圖紙為依據(jù)，在建模時(shí)引入了多種理論假設(shè)，并且存在邊界條件近似以及節(jié)點(diǎn)連接剛度模擬不恰當(dāng)?shù)葐栴}，尤其當(dāng)橋梁服役若干年之后，構(gòu)件開裂等病害使得橋梁剛度下降，結(jié)構(gòu)整體的力學(xué)性能發(fā)生變化。為了更好地模擬舊橋的力學(xué)性能，為日后的監(jiān)測、保養(yǎng)以及加固提供參考依據(jù)，需要依靠模型修正技術(shù)來校正按照圖紙建立的理論有限元模型。

對于模型修正問題而言，可以將數(shù)值分析軟件看成一個(gè)輸入輸出系統(tǒng)，即作為輸出量的目標(biāo)函數(shù)和約束函數(shù)可看成是待優(yōu)化參數(shù)的函數(shù)。通過有限元軟件計(jì)算一次這個(gè)函數(shù)需要很長的時(shí)間，不利于很多啟發(fā)式算法的使用。如果可以通過建立一個(gè)模型繞過有限元計(jì)算，而直接擬合這兩個(gè)函數(shù)，則可以大大提高計(jì)算效率。Kriging模型作為一種無偏估計(jì)模型，具有出色的非線性擬合能力，在特定的隨機(jī)過程中，能夠給出最優(yōu)線性無偏估計(jì)。在現(xiàn)有研究中，Kriging模型大量應(yīng)用在地質(zhì)、水文、氣象、航空航天以及飛行器設(shè)計(jì)領(lǐng)域[1]。

2 Kriging代理模型優(yōu)化算法

2.1 Kriging代理模型

首先在優(yōu)化空間通過抽樣方法，得到n個(gè)樣本點(diǎn):

對這n個(gè)樣本點(diǎn)進(jìn)行有限元分析，得出n個(gè)分析結(jié)果:

對這些結(jié)果進(jìn)行插值，即:

只需要給出加權(quán)系數(shù)ω＝[ω1ω2…ωn]，便可得到任意參數(shù)對應(yīng)的真實(shí)響應(yīng)值。為此，Kriging模型先將預(yù)測值和輸入值之間的關(guān)系假設(shè)為:

式中，μ為常數(shù)項(xiàng)；z(x)是按N(0，σ2)分布的隨機(jī)變量。前者表示預(yù)測值的整體趨勢值，后者反映了預(yù)測值的局部波動(dòng)偏差。

在優(yōu)化空間的不同位置處，隨機(jī)變量存在一定的相關(guān)性，對于隨機(jī)過程z(x)，任意兩點(diǎn)隨著距離的接近，將趨于相同的函數(shù)值。Kriging模型將這種相關(guān)性表示為:

式中，θk為不同維度上相關(guān)性的衰減速率，決定了模型的擬合性能。

基于以上假設(shè)來尋找最優(yōu)的加權(quán)系數(shù)，使得均方誤差MSE最小[2]:

可以得到未知點(diǎn)x處Kriging模型的預(yù)測值為:

采用極大似然估計(jì)的方法尋找θk，將問題轉(zhuǎn)化為尋找θk使式(8)的值最大[3]:

采用全局搜索算法尋找最優(yōu)的θk，θk也隨著樣本庫的更新而更新來獲得更好的性能。

2.2 優(yōu)化加點(diǎn)策略

建立好初始的代理模型之后，下一步就是選擇添加新的樣本點(diǎn)，以較少的樣本點(diǎn)最大程度地提高擬合精度，同時(shí)不斷添加樣本點(diǎn)直到優(yōu)化收斂[4]。

目前國際上已經(jīng)發(fā)展出多種基于代理優(yōu)化方法的加點(diǎn)準(zhǔn)則，有代理模型最優(yōu)預(yù)測準(zhǔn)則(MSP準(zhǔn)則)、改善期望準(zhǔn)測(EI準(zhǔn)則)、均方誤差準(zhǔn)則(MSE準(zhǔn)則)等。

MSP準(zhǔn)則尋找的樣本點(diǎn)是可能的最優(yōu)解，可以快速提高最優(yōu)解附近的精度；MSE準(zhǔn)則尋找的樣本點(diǎn)是誤差期望值最大的點(diǎn)，可以快速提高全局精度。結(jié)合兩種準(zhǔn)則的特點(diǎn)，通過計(jì)算均方誤差最大的點(diǎn)更新預(yù)測誤差最大的區(qū)域，同時(shí)計(jì)算更新代理模型預(yù)測的最優(yōu)解。在迭代的前期，計(jì)算預(yù)測均方誤差最大的點(diǎn)可以快速改善模型的精度；在迭代次數(shù)足夠多時(shí)，可以認(rèn)為模型擬合的精度已經(jīng)達(dá)到要求，這時(shí)不再計(jì)算均方誤差最大的點(diǎn)而只計(jì)算代理模型預(yù)測的最優(yōu)解進(jìn)行加點(diǎn)，加快迭代收斂的速度[5]1716-1718。

2.3 尋找篩選更新點(diǎn)

首先選取好優(yōu)化目標(biāo)，針對每一個(gè)目標(biāo)，分別建立一個(gè)代理模型。根據(jù)兩個(gè)加點(diǎn)準(zhǔn)則，建立兩個(gè)多目標(biāo)優(yōu)化問題[6]。

本文采用NSGA-Ⅱ方法求解多目標(biāo)優(yōu)化問題，應(yīng)用NSGA-Ⅱ算法求解式(9)和式(10)分別得到候選的更新樣本點(diǎn)集H1和H2，為了使Pareto解更加均勻，根據(jù)極大極小距離準(zhǔn)則來選擇正式的更新點(diǎn):

選擇的正式樣本與已有樣本的最小距離，在待更新樣本集中是最大的[7]，即:

式中，x*為選擇的樣本點(diǎn)。

更新樣本集S中，再尋找H1中極大極小距離對應(yīng)的樣本點(diǎn)，這兩個(gè)樣本點(diǎn)就是篩選出的正式更新點(diǎn)。如果沒有求解式(10)，則從求解式(9)得到的樣本集中貫序選擇兩個(gè)點(diǎn)作為正式更新點(diǎn)[5]1720。

2.4 代理優(yōu)化方法

多目標(biāo)優(yōu)化方法具體步驟如下[8]:

(1)首先確定好優(yōu)化計(jì)算的迭代次數(shù)Tm，加點(diǎn)時(shí)計(jì)算預(yù)測方差最大點(diǎn)的次數(shù)Tmse和初始樣本數(shù)n，定義優(yōu)化變量、約束條件、優(yōu)化目標(biāo)。

(2)運(yùn)用拉丁超立方抽樣進(jìn)行實(shí)驗(yàn)設(shè)計(jì)，通過有限元軟件精確求解得到初始樣本集S，從初始樣本集S中篩選出Pareto最優(yōu)解集A。

(3)基于樣本集S構(gòu)建Kriging代理模型。

(4)如果迭代次數(shù)小于等于Tmse，則使用NSGAⅡ算法分別求解本文定義的式(9)和式(10)，分別得到候選樣本集H1和H2；如果加點(diǎn)個(gè)數(shù)大于Tmse，則只求解式(9)，得到候選樣本集H2。

(5)如果有H1，則先篩選H1中的點(diǎn)為樣本點(diǎn)，計(jì)算H1中每一個(gè)點(diǎn)與已有樣本點(diǎn)之間的距離的最小值，每個(gè)候選樣本點(diǎn)都有一個(gè)對應(yīng)的最小值，從這些最小值中選取最大的那個(gè)對應(yīng)的候選點(diǎn)加入樣本集S；再遵循這一準(zhǔn)則從H2中選取一個(gè)點(diǎn)加入樣本集S和最優(yōu)解集A。如果沒有計(jì)算H1，則從H2中選取兩次，得到兩個(gè)點(diǎn)加入樣本集S和最優(yōu)解集A。最后維護(hù)最優(yōu)解集A，只保留非支配個(gè)體。

(6)判斷終止條件，若迭代次數(shù)大于Tm，則輸出A作為最優(yōu)解集合；否則，則轉(zhuǎn)到步驟(3)繼續(xù)迭代。

3 有限元模型

3.1 工程背景

背景工程上構(gòu)為5×16 m預(yù)制鋼筋混凝土簡支T梁，下部結(jié)構(gòu)為柱式墩、樁柱式臺(tái)、樁基礎(chǔ)，如圖1所示。設(shè)計(jì)荷載:汽車—20級(jí)，掛車—100級(jí)(1972標(biāo)準(zhǔn))。

圖1 橋型示意

3.2 建立模型

利用有限元軟件對橋梁第一跨的上部結(jié)構(gòu)建模，主梁、橋面鋪裝以及路緣石等使用實(shí)體結(jié)構(gòu)單元模擬，鋼筋采用三維桿系單元模擬[9]，如圖2所示。

圖2 有限元模型(第一跨)

3.3 施加荷載

根據(jù)橋梁檢測報(bào)告，橋梁靜載試驗(yàn)加載車采用290 kN雙后軸車進(jìn)行加載，如圖3所示。

圖3 軸重及軸距示意

荷載布置如圖4所示。

圖4 荷載布置示意

4 模型修正

4.1 參數(shù)及目標(biāo)選取

根據(jù)橋梁檢測報(bào)告，已知橋梁結(jié)構(gòu)中存在損傷，實(shí)際結(jié)構(gòu)與圖紙的理論模型的區(qū)別在于理論模型無法模擬出實(shí)際結(jié)構(gòu)由于長期使用而產(chǎn)生的剛度變化，因此選擇4片梁的跨中以及1/4跨撓度(見圖5)和結(jié)構(gòu)整體的一階自振頻率共9個(gè)量為修正目標(biāo)[10]。根據(jù)檢測報(bào)告的損傷情況，選取3、4、5號(hào)梁(梁編號(hào)如圖4所示)跨中混凝土彈性模量e1以及密度ρ1，端部混凝土彈性模量e2以及密度ρ2，1、2、6、7號(hào)梁跨中混凝土彈性模量e3以及密度ρ3，端部混凝土彈性模量e4以及密度ρ4，橫隔板混凝土彈性模量e5，共計(jì)9個(gè)變量的系數(shù)作為修正參數(shù)，其變化區(qū)間設(shè)[0.2，2]，據(jù)此可以建立以下最優(yōu)化問題:

圖5 跨中及端部混凝土位置示意

式中，yn(X)為數(shù)值模型計(jì)算的撓度以及頻率與檢測結(jié)果的殘差；X為修正系數(shù)。

考慮到混凝土中箍筋、架立鋼筋等的布置以及混凝土后期強(qiáng)度的增長，并結(jié)合現(xiàn)場觀察橫隔板等結(jié)構(gòu)的損壞情況，各參數(shù)初始取值定為:

對有限元模型進(jìn)行靜力分析以及特征值分析，得到中梁、邊梁跨中、1/4跨撓度以及結(jié)構(gòu)一階自振頻率[11]，見表 1。

表1 有限元模型修正前結(jié)果與檢測結(jié)果對比

可以看出，修正前的有限元模型與橋梁的實(shí)際靜動(dòng)力特性存在較大差異。

4.2 修正結(jié)果

在進(jìn)行了76次迭代后，最終得到一個(gè)包括15個(gè)最優(yōu)解的解集。為評(píng)估修正后模型的準(zhǔn)確性，選取偏載工況下的位移檢測值與有限元模型計(jì)算值進(jìn)行對比[12]，選取其中一組結(jié)果如表2、表3所示。

表2 修正后結(jié)果與檢測結(jié)果對比

表3 修正后偏載結(jié)果與檢測結(jié)果對比

結(jié)合有限元分析結(jié)果以及橋梁檢測報(bào)告，第一組修正參數(shù)相較于其他幾組更加符合橋梁實(shí)際狀況。

5 結(jié)論

(1)修正后的模型計(jì)算的最大相對誤差為3%，修正后的模型能夠較為精確地反映該橋的現(xiàn)有狀態(tài)，為后期該橋的加固計(jì)算提供了一個(gè)基準(zhǔn)模型。

(2)改良后的Kriging代理優(yōu)化方法可以通過較少的有限元計(jì)算次數(shù)得到較好的修正效果，且能夠忽略檢測數(shù)據(jù)中的奇異值，具有良好的魯棒性。

(3)相較于將多目標(biāo)問題解決方法轉(zhuǎn)化為單目標(biāo)問題進(jìn)行求解，直接進(jìn)行多目標(biāo)優(yōu)化能夠獲得多個(gè)可行的修正結(jié)果，結(jié)果之間具有均勻性、多樣性。能結(jié)合工程實(shí)際和進(jìn)一步的檢測結(jié)果選擇較合理的參數(shù)，從而保證修正之后的結(jié)果不失真。