最優(yōu)階次多項式響應面法及其在模型修正中的應用

王涵，史治宇

(南京航空航天大學 機械結(jié)構(gòu)力學及控制國家重點實驗室，江蘇 南京 210016)

0 引言

有限元方法是一種運用離散思想對連續(xù)性問題進行近似求解的數(shù)值方法，在工程領(lǐng)域中應用廣泛[1]。實際工程問題中由于邊界條件、連接等簡化建模，模擬參數(shù)常常存在誤差，有限元模型往往不能準確反映真實模型情況。因此有限元模型(finite element model,FEM)修正技術(shù)應運而生，旨在通過試驗數(shù)據(jù)修正有限元模型，從而獲得更符合真實情況的有限元模型。經(jīng)典參數(shù)型模型修正理論中的迭代計算需反復調(diào)用有限元軟件，尤以靈敏度矩陣的計算耗費大量時間[2]，因此在處理實際工程問題時存在計算效率低的問題。為了減少計算耗時同時保證計算精度，響應面法(response surface method,RSM)開始應用于模型修正領(lǐng)域。響應面法是統(tǒng)計學中的綜合試驗技術(shù)，通過擬合初始模型中輸入與輸出之間的隱式關(guān)系，獲得具有顯式函數(shù)關(guān)系的響應面模型[3]。采用響應面法將初始模型簡化為代理模型進行的修正計算，避免了靈敏度矩陣模塊的調(diào)用，簡化了輸入、輸出間的軟件計算，對動力學分析和優(yōu)化具有重要意義。

響應面發(fā)展出多種近似方法，例如多項式法、神經(jīng)網(wǎng)絡(luò)法、徑向基法等。針對不同問題背景選取合適的響應面建立方法，有益于提高代理模型精度。相較于其他響應面近似方法，多項式法因其計算更簡便、可處理大多數(shù)非線性問題的優(yōu)勢而得到了廣泛應用[4]。常規(guī)的多項式響應面法通常采用二階完全多項式進行構(gòu)建，階次和變量形式固定，存在精度較低、無法準確反映高階模型的缺點[2]。面對非線性程度較高的輸入、輸出，可以通過適當提升多項式階次來達到更精確的擬合。然而隨著階次增高，參數(shù)組成的變量將呈指數(shù)爆炸增長；各交叉項[5]的顯著度存在差異，取舍不當會影響擬合效果；最高階次選擇不當，易產(chǎn)生過擬合的問題，且傳統(tǒng)多項式中所有參數(shù)擬合階次相同，忽略了不同參數(shù)對輸出映射之間的差異。因此盲目采用高階完全多項式，可能出現(xiàn)計算時間增長、模型精度下降等情況，違背了使用代理模型的初衷。同時，代理模型的精度較低時，會引發(fā)模型修正誤差的增大，因此建立高精度的響應面模型是其應用于模型修正的前提。綜上，可以探尋一種區(qū)分不同參數(shù)對輸出的最優(yōu)階次，使得參數(shù)間的擬合階次不再固定的多項式響應面構(gòu)造方法，并將其應用于模型修正領(lǐng)域。

本文在此提出一種最優(yōu)階次多項式響應面的構(gòu)造方法，在傳統(tǒng)2階多項式基礎(chǔ)上提升并區(qū)分了參數(shù)對響應的擬合階次，以突破階次固定統(tǒng)一的限制，提升代理模型的精度。同時將改進的響應面法應用于模型修正領(lǐng)域，以解決經(jīng)典修正過程中耗時長的問題，并保證修正結(jié)果的精度。本文首先闡述了最優(yōu)階次響應面的構(gòu)造方法，包括最優(yōu)階次判定、多項式構(gòu)造以及精度檢驗，接著建立了運用最優(yōu)階次響應面進行模型修正的基本理論。以功能函數(shù)為例，采用改進后的響應面有效實現(xiàn)了函數(shù)模型的還原；以GARTEUR飛機模型為例,以真實試驗數(shù)據(jù)為目標值,利用基于最優(yōu)階次多項式響應面法的有限元模型修正方法對彈簧單元剛度等3個設(shè)計參數(shù)進行了有效修正。

1 基本理論

獲得高精度響應面模型的基礎(chǔ)是試驗設(shè)計，即選擇合適的結(jié)構(gòu)模型參數(shù)試驗樣本點，使響應面擬合函數(shù)在整個參數(shù)空間內(nèi)能最有效地描述參數(shù)和響應之間的隱式關(guān)系。試驗設(shè)計方法包括析因?qū)嶒炘O(shè)計、中心復合設(shè)計 (CCD)以及Box-Behnken設(shè)計等。建立參數(shù)與響應間的顯式關(guān)系即對樣本進行擬合。常規(guī)的多項式擬合一般采用2階完全多項式，表達式如下:

(1)

其中α0、αi、αii、αij為多項式系數(shù)。為了更適應于樣本特性，提升擬合效果，可以通過判斷單個參數(shù)對響應的最優(yōu)擬合階次，構(gòu)造一種區(qū)分參數(shù)擬合階次的多項式響應面，具體理論如下。

1.1 最優(yōu)階次判定

最優(yōu)階次是指樣本數(shù)據(jù)中單參數(shù)對響應的最佳擬合階次。設(shè)有m個參數(shù)與n個響應，共判斷m×n個最優(yōu)階次。單參數(shù)x對單響應y的最優(yōu)階次，通過逐階增添回歸變量中冪次項的多元回歸對比分析。

設(shè)階次判定上限為k，從第2階起，冪次項每增加1階，建立共k-1個回歸方程。各階回歸方程見式(2)。

(2)

其中a0、a1、ai為系數(shù)。共得到k-1個回歸函數(shù)，即單變量對單響應的2-k階回歸函數(shù)組。

最佳擬合函數(shù)的選擇采用擬合優(yōu)度與顯著F檢驗結(jié)合的方法。擬合優(yōu)度指回歸函數(shù)對觀測值的擬合程度，由可決系數(shù)R2來度量。R2的值越接近1，回歸函數(shù)對樣本值的擬合程度越好。顯著性F檢驗則為判定回歸假設(shè)為0的概率，并用p值來量化。p值越小，則越拒絕，即回歸越可靠。為保證整體回歸的精度和回歸系數(shù)的顯著，將兩種指標結(jié)合，取判定指標δ：

δ=R2-p

(3)

其中：R2為可決系數(shù)；p為顯著性概率，可通過擬合工具軟件計算。對每階回歸函數(shù)進行相應δ值計算，可建立k-1維數(shù)組D：

D=[δ1;δ2;δ3;…;δk-1]

(4)

D最大值對應的擬合函數(shù)回歸精度最高，回歸系數(shù)最顯著，其冪次項最高階次即為單參數(shù)對單響應的最優(yōu)擬合階次τ：

δi=Dmax

(5)

τ=i+1 (τ≥2)

(6)

通過對每個參數(shù)和響應進行階次τ判斷，組成最優(yōu)階次矩陣Ord：

(7)

Ord針對不同參數(shù)對不同響應的階次做到了逐個區(qū)分，但由于試驗設(shè)計獲取的樣本數(shù)據(jù)通常沒有進行單變量控制，最優(yōu)階次的判定并不能反映參數(shù)與響應之間的物理特性關(guān)系。

1.2 多項式響應面構(gòu)造

由于不同參數(shù)之間階次存在差異，高階交叉項的組合也變得十分復雜，固定形式的多項式不再適用。為了保證多項式構(gòu)造方法的簡便性，在式(1)的基礎(chǔ)上，增添所有參數(shù)的2階及2階以上冪次項，形式如下：

(8)

其中：α0、αi、αij為多項式系數(shù)；βiε為參數(shù)i的第ε冪次項系數(shù)。若Ord各項均為2，則最優(yōu)階次多項式即為2階完全多項式。

根據(jù)式(7)-式(8)，共構(gòu)造n個響應的多項式函數(shù)，代入樣本點計算得相應系數(shù)值后，該組函數(shù)表達式即為最優(yōu)階次多項式響應面模型Y：

Y=[y1;y2;y3;…;yn]

(9)

最優(yōu)階次多項式實現(xiàn)了參數(shù)間階次的提升與區(qū)分，避免了高階交叉項的復雜組合，突破了固定的階次和變量形式，在提高精度的同時控制了計算量。初步建立響應面模型后需進行精度檢驗。

1.3 精度檢驗

響應面模型可以從兩方面檢驗精度，一是統(tǒng)計意義上，可將擬合優(yōu)度可決系數(shù)R2或均方誤差MSE作為檢驗標準。R2值越接近1，得到的響應面模型越接近有限元模型；MSE則正好相反：

(10)

(11)

二是從力學角度進一步驗證，這部分可以結(jié)合其在有限元模型修正中的應用來進行。

最優(yōu)階次多項式響應面在構(gòu)建過程中非常依賴樣本數(shù)據(jù)，不同數(shù)據(jù)可能得到不同階次判定結(jié)果。當精度檢驗結(jié)果不夠理想時，應返回重新進行試驗設(shè)計。

1.4 模型修正中的應用

應用響應面法進行模型修正時，模型修正可看作優(yōu)化問題：

(12)

其中：X為待修正參數(shù)；Ye為試驗值；Y(X)為響應面計算值；VLB和VUB為X的設(shè)計范圍。求解如下：

R=SΔX

(13)

其中：R為響應殘差；S為靈敏度；ΔX為參數(shù)殘差。有限元軟件中靈敏度采用Nelson法計算[6]，而響應面模型建立后可跳過有限元軟件的調(diào)用，利用擬合函數(shù)的偏導計算靈敏度，求解如下：

(14)

修正后若誤差收斂，則修正完成，否則需重新進行試驗設(shè)計構(gòu)造響應面。將修正后參數(shù)值代入初始有限元模型中，獲得修正后的有限元模型。最優(yōu)階次多項式響應面法應用于有限元模型修正的流程總結(jié)如圖1所示。

圖1 最優(yōu)階次多項式響應面法應用于有限元模型修正流程圖

2 算例

2.1 功能函數(shù)

功能函數(shù)[7]定義如下:

其中：X1∈(0,1)；X2∈(0,1)。

功能函數(shù)響應面如圖2所示。

圖2 功能函數(shù)響應面

通過拉丁超試驗設(shè)計，獲取20組樣本點，以此進行改進前后兩種響應面的擬合。經(jīng)過式(1)計算，建立2階完全多項式響應面，如圖3所示。

圖3 2階完全多項式回歸響應面

矩陣中顯示，參數(shù)1對響應的階次提高到4，參數(shù)2對響應的階次提高到5。接下來根據(jù)式(8)-式(9)建立最優(yōu)階次多項式響應面Y，如圖4所示。

圖4 最優(yōu)階次多項式回歸響應面

兩種響應面初步建立后，可以從圖像上直觀看出，最優(yōu)階次多項式更好地還原了功能函數(shù)的趨勢分布，但在邊界角點處仍存在可見差異。分析認為是由于拉丁超試驗設(shè)計方法中遵循分層概率區(qū)間取樣，樣本點落在角點邊界附近較少，影響了邊界處的擬合精度。

以R2和MSE為統(tǒng)計指標對兩種響應面的精度進行檢驗，相關(guān)計算見式(10)-式(11)。數(shù)據(jù)對比見表1。

表1 兩種多項式響應面模型精度對比

表1數(shù)據(jù)顯示，最優(yōu)階次多項式響應面的擬合優(yōu)度非常逼近于1，相較于2階完全多項式提升了29.1%，且MSE更加趨近于0。表明該改進的響應面具有更高的精度，更真實地還原了功能函數(shù)。由于樣本點處均達到了良好擬合，所以角點的誤差并未影響整體精度。

2.2 GARTEUR飛機

GARTEUR飛機模型為法國國家航天研究院(ONEDR)設(shè)計制造,被定為評估有限元模型修正技術(shù)的基礎(chǔ)模型。該飛機具有機身、機翼、平尾和垂尾4個主體部分，材料主體為鋁，翼展長2 m、機身長1.5 m，總質(zhì)量38.3 kg。GARTEUR 模型照片如圖5所示，建立GARTEUR初始有限元模型如圖6所示。

圖5 GARTEUR模型

圖6 GARTEUR有限元模型

初始有限元模型共46個二維梁單元，包含4種梁屬性，分別為機身的梁單元、機翼的梁單元、平尾的梁單元以及垂尾的梁單元。另有8個彈簧元，包含5種彈簧元屬性，如圖6所示，包括在機翼與機身連接處1，用于模擬左、右機翼對機身繞x軸方向抗彎強度的彈簧單元；在垂尾與機身的連接處2，分別用于模擬垂尾繞x軸和繞z軸方向抗彎強度的兩種彈簧單元；在平尾與垂尾的連接處3，分別用于模擬垂尾繞x軸和繞y軸方向抗彎強度的兩種彈簧單元。有限元模型共52個節(jié)點，312個自由度。對GARTEUR模型進行有限元分析，其前5階固有頻率見表2。

表2 初始有限元模型前5階固有頻率 單位：Hz

其中第1階主要為機翼繞x向的一彎，第2階主要為垂尾繞x向一彎，第3階主要為垂尾繞y向的一彎，第4階主要為機翼繞x向的反對稱二彎和平尾繞z向轉(zhuǎn)動的組合，第5階主要為平尾繞z向的轉(zhuǎn)動。

根據(jù)實測試驗結(jié)果，模型的實測前5階固有頻率與有限元計算值存在誤差，目標值與初始值的相關(guān)數(shù)據(jù)對比見表3。

表3 前5階固有頻率目標值與初始值對比

為修正頻率誤差，應選取合適參數(shù)對有限元模型進行修正。考慮到連接剛度不易準確建模，在此選取機翼與機身連接處繞x向轉(zhuǎn)動的連接剛度pelas11、平尾和垂尾連接處繞y向轉(zhuǎn)動的連接剛度pelas15、與垂尾與機身連接處繞x向轉(zhuǎn)動的連接剛度pelas16三個彈簧單元剛度作為三個待修正參數(shù)，位置如圖6所示，初始建模數(shù)據(jù)如表4所示。

表4 彈簧單元剛度初始建模值

依據(jù)圖1流程，對GARTEUR模型進行基于最優(yōu)階次多項式響應面法的模型修正。下文中將把3種彈簧單元簡稱為3個參數(shù)，5階固有頻率簡稱為5個響應。

設(shè)定3個參數(shù)合理的設(shè)計范圍，分別為：[1.8,4.2]、[1.56,3.64]、[1.02,2.38]。以5階固有頻率為響應，進行中心復合試驗設(shè)計，獲得15組樣本數(shù)據(jù)。結(jié)合工程實際情況，設(shè)定階次判定上限k為5。經(jīng)過式(2)-式(7)計算，獲得最優(yōu)階次矩陣Ord：

結(jié)果顯示，參數(shù)1對響應2、響應3的階次提高到4，參數(shù)2對響應4的階次提高到3，參數(shù)3對響應1、響應4、響應5的階次分別提高到4、4、3。接下來根據(jù)式(8)-式(9)建立最優(yōu)階次多項式響應面Y。

在此選取第4階固有頻率作為典型響應進行回歸響應面繪制，見圖7、圖8。響應面模型計算第4階固有頻率范圍為[25，29] Hz，真實試驗數(shù)據(jù)中第4階固有頻率為26.62 Hz，說明參數(shù)上下限取值范圍合理。

圖7 第4階固有頻率對參數(shù)1、參數(shù)3回歸響應面

圖8 第4階固有頻率對參數(shù)2、參數(shù)3回歸響應面

根據(jù)式(14)，對初步建立響應面模型進行參數(shù)靈敏度計算，結(jié)果如圖9所示。

圖9 響應面模型靈敏度計算結(jié)果

響應面初步建立后，以R2為指標對最優(yōu)階次多項式響應面的精度進行檢驗，相關(guān)計算見式(10)。同時以同組樣本數(shù)據(jù)建立2階完全多項式響應面，將其精度作為對照，數(shù)據(jù)對比如表5。

表5 兩種多項式響應面模型精度對比

表5數(shù)據(jù)顯示，最優(yōu)階次多項式響應面的5個響應擬合優(yōu)度都非常逼近于1，表明該響應面完成了良好的數(shù)據(jù)擬合，且在部分階次中最優(yōu)階次多項式模型較2階完全多項式擬合優(yōu)度略有提升，表明該響應面模型具有更高的精度。由于該算例模型簡單，參數(shù)較少，運用2階完全多項式擬合也可達到較高精度，且兩種響應面間參數(shù)階次差異較小，因此精度提升空間并不高。

接下來利用真實試驗數(shù)據(jù)對最優(yōu)階次多項式響應面模型進行修正。相關(guān)數(shù)據(jù)代入式(12)-式(14)，經(jīng)過4次迭代計算后誤差收斂，獲得修正后結(jié)果。修正前后參數(shù)值對比見表6，修正前后響應值對比見表7。

表6 修正前后參數(shù)值對比

表7 修正前后響應值對比

表6數(shù)據(jù)顯示，參數(shù)修正后變化幅度在20%以內(nèi)，沒有超出參數(shù)設(shè)計變化范圍。表7數(shù)據(jù)顯示，頻率的平均誤差從修正前的3.52%下降到修正后的1.79%，說明最優(yōu)多項式響應面模型可以在模型修正中進行有效的模型替代和良好的優(yōu)化計算。其中響應1和響應5的誤差略有放大，分析是由于部分參數(shù)靈敏度趨近0，可能引起矩陣奇異、誤差不收斂的情況。

為進一步檢驗最優(yōu)階次多項式響應面在模型修正中的結(jié)果精度，在此進行3種修正情況的對比：不使用響應面的經(jīng)典模型修正，使用2階完全多項式響應面的模型修正以及使用最優(yōu)階次多項式響應面的模型修正。達到誤差收斂時3種修正方法分別迭代了8次、2次和4次。相關(guān)數(shù)據(jù)對比見表8。

表8 3種修正情況的結(jié)果誤差對比 單位：%

表8數(shù)據(jù)顯示，初始有限元模型修正前的平均誤差為3.52%，不使用響應面修正的平均誤差為2.81%，使用2階完全多項式響應面模型修正的平均誤差為3.30%，而使用最優(yōu)階次多項式響應面模型修正的誤差達到1.79%，表明3種修正方法都可以有效修正初始模型頻率的誤差。3種修正方法中都存在少數(shù)響應誤差放大的情況，因為部分參數(shù)靈敏度接近于0，會產(chǎn)生矩陣奇異誤差不收斂的情況。對比可得：1)修正中使用響應面與不使用響應面相比，無需調(diào)用有限元軟件，且迭代次數(shù)下降。但2階多項式響應面修正結(jié)果誤差較不使用響應面的修正結(jié)果有所提高，分析認為是代理模型的精度不夠而引發(fā)了誤差。2)使用最優(yōu)階次多項式響應面與完全2階多項式響應面相比，迭代次數(shù)少量增加，但修正結(jié)果更接近真實數(shù)據(jù)，并且少數(shù)響應誤差放大情況得到了改善，計算效率和修正結(jié)果精度更高。

3 結(jié)語

本文提出了一種最優(yōu)階次多項式響應面法的構(gòu)建方法，解決了常規(guī)多項式響應面法階數(shù)限制的問題，區(qū)分了多項式中不同參數(shù)間的階次特性，提高了響應面模型的精度。將改進的響應面法應用于模型修正中，降低了經(jīng)典修正過程中反復調(diào)用有限元軟件靈敏度模塊消耗的計算量，減少了由于代理模型精度不夠引發(fā)的修正結(jié)果誤差。數(shù)值算例表明該改進響應面法可以提高傳統(tǒng)響應面的精度；GARTEUR算例表明該響應面法可以成功應用于模型修正領(lǐng)域，并通過提高代理模型的精度而有效減少修正誤差，提高修正效率，具有良好的應用前景。