快速交叉驗證改進的運載火箭近似建模方法

文謙,楊家偉,武澤平,楊希祥,趙海龍,王志祥

國防科技大學 空天科學學院,長沙 410073

運載火箭總體設計是涵蓋多學科的復雜系統工程,在其研制過程中需綜合考慮氣動、載荷、結構、動力等多個因素。隨著計算機技術的飛速發展,智能化設計方法逐漸成為主流,也對設計效率及精細化程度提出更高要求。同時工程仿真中常用的計算流體力學(Computational Fluid Dynamics, CFD)和有限元分析在復雜構型下計算耗時長,資源消耗多。因此亟需發展一套高精度,高效率的近似建模方法來縮短整箭的研制周期。近似建模是根據離散的訓練樣本,采用插值或擬合的方法對模型在全變量域的輸出進行預測的技術。在復雜工程問題進行分析和優化的過程中,通常采用近似模型作為原復雜耗時模型的代替,達到減少仿真次數,提升設計效率的目標。

徑向基函數(Radial Basis Function,RBF)模型由Hardy在1971年首次提出,是以待測點與樣本點之間的歐氏距離為自變量的函數,通過簡單基函數的線性組合來模擬復雜非線性模型的響應特性,通過求解線性方程組得到基函數的系數,實現對設計空間內任意輸入參數的預測。RBF近似模型由于其原理簡單、操作方便,已成為眾多復雜耗時的工程設計問題中廣泛使用的近似模型之一。Jin等通過14個代表不同問題的算例對包括RBF在內的4種方法進行了對比研究,結果表明RBF的精度和魯棒性最為可靠。為快速而又準確地獲取可信的近似模型,通常需要對基函數的形狀參數進行合理確定,以提升RBF近似模型預測精度。常用的形狀參數確定方法有按經驗直接確定法和優化形狀參數法等。直接確定法受樣本點在空間內的分布影響,Kitayama和Yamazaki提出的非均勻形狀參數確定法是常用方法之一。優化形狀參數法通過定義一定的精度評估指標(通常為交叉驗證誤差),將各基函數的形狀參數作為設計變量,尋求預測精度最高的形狀參數組合,能夠得到最合理的形狀參數取值,但由于每個基函數都有一個形狀參數,隨著訓練樣本的增多,優化問題的規模也相應增加。

交叉驗證是一種統計學和機器學習領域估計模型泛化能力的方法,利用交叉驗證可以解決過擬合和欠擬合的問題,避免使用同一數據集既進行模型訓練,又進行泛化誤差估計,能夠有效合理的評估近似模型的預測精度。交叉驗證方法由Wold首次提出,其基本思想是將原始數據分組,一組作為訓練集訓練模型,另一組為測試集來對模型做出評價,以驗證模型精度和泛化能力。在模型存在未知參數時,通過優化未知模型參數,使得交叉驗證誤差最小化,可實現預測模型與訓練數據的匹配,從而提升預測模型性能。Van Gestel等證明了將交叉驗證法用于支持向量機(Support Vector Machine, SVM)和最小二乘支持向量機分類器(Least Squares SVM classifiers, LS-SVMs)的超參數選擇,兩者可以獲得相似的結果,為交叉驗證方法在近似建模領域的應用奠定了堅實的理論基礎。

目前近似建模領域常用的交叉驗證方法主要有簡單交叉驗證(Hold-out Cross-Validation, Hold-out CV)法、K-折交叉驗證方法(K-fold Cross-Validation, K-fold CV)法和留一交叉驗證(Leave-One-Out Cross-Validation, LOOCV)法。Hold-out CV法是將數據集分為2個子集,一個子集為測試集,而另一個子集則為訓練集。該方法對數據的處理較為簡單,但對于原始數據進行隨機分組,得到的結果不精確,且不同分組情況所得到的結果波動很大,驗證結果不可復制。為了確保驗證結果可復制,降低驗證過程隨機性,提高樣本利用率,可選用LOOCV法。LOOCV法是將數據集中的個樣本輪流作為測試集,余下的-1個樣本作為訓練集,按照此方法將得到個模型,計算次誤差。LOOCV法的優點是幾乎將所有的樣本用于擬合模型,得到的結果近似無偏,不受數據集劃分方法的影響。該方法雖然是對樣本利用率最高的方法,但其計算成本過高,面對較大的數據集時,缺點明顯。通常在樣本容量較小時才會選用LOOCV法。近年來國內的張英堂、劉學藝、李佳等逐漸開展了針對留一交叉驗證法求解效率提升的研究。但LOOCV法僅是K-折交叉驗證法的特例,面對不同問題的響應特性和不同的樣本點散布,通用K-折交叉驗證方法是近似建模中不可回避的問題。目前針對任意的K-fold CV法求解效率問題鮮有研究。

本文利用樣本點局部密度實現將多個形狀參數確定轉換為縮放系數單一參數確定,建立通用快速交叉驗證方法實現縮放系數的高效求解,應用于運載火箭智能設計過程中,提高近似模型求解效率及預測精度。

1 基于樣本局部密度的徑向基函數形狀參數計算方法

1.1 徑向基函數近似模型

基本徑向基函數的數學模型為

(1)

式中：為待預測樣本；為采樣點個數；為每個基函數的權系數；()是以未知樣本到已知樣本的歐式距離為自變量的基函數；表示樣本點距樣本點的歐氏距離,即

(2)

如式(3)所示,基函數()選用高斯函數

(3)

式中：為形狀參數。的值確定后,通過訓練集合及其輸出代入式(1)中,求解線性方程組即可計算出當前訓練樣本對應的基函數權系數矩陣

=

(4)

(5)

基函數系數決定后,針對任意樣本,利用式(1)可以得到相應的預測輸出。任意不為0的代入上述過程,均可得到經過所有訓練樣本的近似模型,但不合理的會導致模型泛化能力差,所以需要合理確定。本文在分析形狀參數直觀意義的基礎上,提出了基于局部密度的形狀參數歸一化表征方法。

1.2 高斯徑向基函數形狀參數歸一化表征方法

RBF的核心過程在于確定形狀參數的值,直接影響數值試驗的誤差結果。而Rippa指出交叉驗證是估算形狀參數的一種有效方法,形狀參數的合理選擇能夠有效提高近似模型的精度,減小訓練樣本個數,指導下一步的采樣優化。李響和王曉鵬指出形狀參數越大,基函數越光滑,近似模型也越光滑,有利于提高近似模型精度。Wu等的研究表明,隨著值的不斷增大,模型的精度呈現先上升后下降的趨勢,且由于不同真實模型性質不同,因此也難以通過選擇統一的形狀參數保證不同近似模型的精度。

對于RBF而言,未知樣本點的預測結構由所有已知樣本點的輸出共同影響決定,而形狀參數則控制著當前樣本點在設計空間內向周圍影響衰減的快慢。如圖1所示,形狀參數越小,樣本點對周圍空間的影響衰減越快,反之,則對周圍的影響能夠傳播更遠。

圖1 不同形狀參數下Gauss基函數曲線Fig.1 Gaussian radial basis function curves with different shape parameters

樣本點在空間內并不是均勻分布的,所以不能采用單一的形狀參數。應根據樣本點分布情況,選擇合適的形狀參數。在樣本密集區域,應采用較小的形狀參數,使每個樣本的影響相對較小；在樣本稀疏區域,應采用較大的形狀參數,適當增大每個樣本的影響范圍。采用樣本局部密度表征樣本疏密程度,形狀參數的選擇就應與樣本局部密度相關,樣本局部密度的表達式如式(6)和式(7) 所示

(6)

(7)

式中：與形狀參數類似,表示樣本點對樣本局部密度響應的影響程度。

在將樣本歸一化到單位超立方體后,通常可取

(8)

式中：為樣本數量；為設計維度。

采用式(6)～式(8)計算出樣本對應的局部密度,得到如圖2所示結果。圖中局部密度曲線表明,樣本越密集,局部密度值越大,樣本越稀疏則值越小。表明該方法可以合理地對離散樣本點的局部密度進行量化。

理想的形狀參數應在樣本密集處取值相對小、樣本稀疏處較大,因此得到樣本點局部密度之后,將徑向基函數形狀參數視為該樣本點影響半徑,令每個樣本點影響范圍(影響體積)與該點局部密度成反比,即

(9)

(10)

(11)

圖2 一維樣本局部密度示意圖Fig.2 Diagram of one-dimensional sampling points local density

1.3 縮放系數選擇方法

(12)

式中：為縮放系數。

如圖3所示,當取較大值時,每個基函數的影響范圍被放大,適用于目標函數變化較平滑的情況；取較小值時,基函數的影響范圍縮小,適用于目標函數局部變化劇烈的情況。圖中實心點為采樣點,細實線代表基函數,虛線區域代表基函數影響范圍,粗實線是疊加后的輸出預測。

圖3 不同縮放系數影響范圍示意圖Fig.3 Diagram of influence range of different scale factors

為合理確定,建立如式(13)所示的優化問題,并采用黃金分割法等一維優化算法對其進行求解以獲得最優縮放系數。

min()

(13)

式中：(·)為縮放系數對應的交叉驗證誤差,用來評估所構建的近似模型精度。

交叉驗證基本思想是將樣本集中的樣本輪流作為測試樣本,以評估近似模型的預示精度。相對于Hold-out CV法和LOOCV法,K-折交叉驗證是最通用的方法,其基本思想是將數據集分為個子集,每個子集輪流作為測試集,其余-1個子集則作為訓練集。LOOCV法是K-fold CV法中=時的特殊情況。所有的樣本都作為訓練集和測試集,且每個樣本都被驗證了一次。同時,K-fold CV法中選取不同的值可能會得到具有差異的結果。其迭代過程如圖4所示。

圖4 K-fold CV迭代圖Fig.4 Iteration diagram of K-fold CV

將原始系數矩陣的行數據劃分為組,每組內有行數據。K-fold CV法誤差Δ為近似模型預測值與真實值之差,即

Δ=pre,-

(14)

式中：pre,為針對-個訓練樣本構建的近似模型在個測試樣本處的預測值；為測試樣本真實輸出。

(15)

由式(15)可見,在K-fold CV中每折都需要訓練一個新近似模型,即對一個-階矩陣進行求逆,整個過程需求解個-階矩陣的逆。和越大、越小,計算速度越慢。線性方程組式(15) 的求解占據了交叉驗證過程大部分的計算時間,消耗了大量的計算資源。

為提升K-fold CV運算效率,節約計算成本,本文提出了逆矩陣快速求解算法和交叉驗證誤差快速求解算法,有效解決了反復求解高階矩陣的逆的問題。采用K-fold CV進行縮放系數的確定,由于一次交叉驗證誤差評估需要對訓練集進行次近似建模,增加了RBF近似建模的計算耗時。針對該問題,提出快速K-折交叉驗證(Fast K-fold Cross-Validation, FK-fold CV)方法以提升近似建模效率與精度。

2 快速交叉驗證方法

2.1 逆矩陣快速求解算法

針對K-折交叉驗證計算效率低下的問題,本文提出的FK-fold CV法通過矩陣初等變換與分塊矩陣求逆來改進常規K-fold CV法,降低計算復雜度,提升直接求解交叉驗證誤差的效率。

如圖5所示,K-fold CV法通過將樣本劃分為組,每次取一組樣本進行測試,剩下的所有樣本作為訓練樣本,循環次之后即得到近似模型在每一組樣本點上的誤差。虛線框中將第組與第組交換的步驟是FK-fold CV法中獨有的。

圖5 快速K-折交叉驗證誤差求解示意圖Fig.5 Flow chart of FK-fold CV

FK-fold CV法以K-fold CV法為基礎,將數據集分為個子集,求解每一子集的交叉驗證誤差時,通過初等變換將測試集置于訓練數據集的末端,變換過程如圖6所示。

圖6 FK-fold CV迭代圖Fig.6 Iteration diagram of FK-fold CV

通過左乘變行,右乘變列的原則對原系數矩陣進行行列交換,將應作為測試集的第組與第組換位后,得到新的系數矩陣。(,)定義為階單位陣第組與第組交換后得到的矩陣,與關系如式(16)所示。

(,)(,)=

(16)

(17)

系數矩陣被劃分為組后,設當前測試組內有個樣本,將第組與第組數據交換位置,式(17)矩陣中的每個元素代表一個×的小矩陣,例,表示數據集中第組數據。對式(16) 求逆得

(18)

由(,)(,)=,可以得到式(19)。

(19)

()=

(20)

(21)

、、、滿足如下數學關系時,可以計算得出

(22)

其中:是(-)階方陣;是(-)×階矩陣;是×(-)階矩陣;是階方陣。采用式(23)即可直接計算出(-),運用階矩陣求逆替代了(-)階矩陣直接求逆的過程,由于?,因此轉化后的矩陣求逆計算量會大幅降低。

(23)

式中：-表示矩陣去掉最后行、列所得矩陣的逆矩陣。

上述推導表明,當LOOCV法是K-fold CV法的特殊情況,即等于樣本點的數量,每次取1個樣本進行測試,其余作為訓練樣本集。當采用留一法進行推導時,式(21)中為一常數,此時近似模型的系數矩陣的逆矩陣可通過式(24) 計算

(24)

2.2 交叉驗證誤差快速求解算法

根據式(15)可得,FK-fold CV法誤差Δ為

Δ=pre,-=,-(-)--

(25)

式中：,-為系數矩陣中測試樣本對應的行和前-列元素；-為訓練樣本預測輸出。由式(23)和式(25),可得式(26)

(26)

(27)

式中：為所有樣本的真實輸出,且測試樣本的個輸出對該項沒有任何貢獻；,為系數矩陣中測試樣本對應的行元素。式(27)展開得

(28)

根據逆矩陣定義有

(29)

(30)

因此,式(28)轉化為

(31)

根據RBF的定義可知

(32)

將式(31)、式(32)代入式(25),可得

(33)

式中：為測試樣本點對應的個權系數。

顯然,采用FK-fold CV法計算預測誤差,僅需要對小矩陣進行求逆,極大地節省了計算量。當每折樣本個數為1時,交叉驗證誤差按照以上過程推導得

(34)

本文針對LOOCV法計算誤差的推導結果,Rippa所推導的快速誤差計算公式是本文取時的特殊情況,表明本方法具有更強的通用性。在計算時間上,對采樣點個數為的近似建模問題,普通K-折交叉驗證的時間復雜度為,快速K-折交叉驗證的時間復雜度為()。采用500個樣本點的數值算例分別對K-fold CV法和FK-fold CV法求解效率進行測試,在不同折數下時間節省幅度如圖7所示。折數越多時間節省幅度越大,尤其在與常規留一交叉驗證對比時,節省了超過90%的時間,算例表明快速K-折交叉驗證具有更高的求解效率。

圖7 時間節省幅度Fig.7 Graph of time saving percentage

3 數值算例校驗

選取4個多峰非線性函數驗證本文近似建模方法的效率、精度和穩定性,算例形式如表1所示,對應的函數圖像如圖8～圖11所示。

圖8 測試函數1圖像Fig.8 Graph of test function 1

圖9 測試函數2圖像Fig.9 Graph of test function 2

圖10 測試函數3圖像Fig.10 Graph of test function 3

圖11 測試函數4圖像Fig.11 Graph of test function 4

為驗證本方法的有效性,將本文方法與Kitayama、Rippa等提出的近似建模方法以及Kriging方法在相同條件下進行性能對比。采用4種近似模型訓練方法分別針對4個測試函數進行訓練,選擇樣本預測輸出與樣本真實輸出的均方根誤差(Root Mean Square Error, RMSE)來表征近似模型的精度。RMSE指標定義如下：

(35)

表1 數值測試函數Table 1 Numerical test function

為避免樣本選取帶來的影響,針對隨機取樣的情況獨立運行30次對方法性能進行測試。基于快速交叉驗證法選取最優縮放系數,對于一維函數在設計域內隨機選取20個樣本點,對于二維函數在設計域內隨機選取40個樣本點,通過樣本點局部密度計算不同樣本點相對形狀參數大小,然后采用一維搜索法針對縮放系數進行尋優,得到最優縮放系數及其對應的預測模型。之后采用相同的樣本點,用其他3種方法構造近似模型并計算對應的RMSE,不同算法得到的均方根誤差箱型圖分別如圖12～圖15所示,其中測試函數3、4由于其誤差較大,將其結果取對數處理后,再用于表征其趨勢。結果表明,基于FK-fold CV法構建的近似模型在4個測試函數中,都展示出比其他3種方法更優越的性能。

本文近似模型構建算法每一次測試都能達到與已知方法相近甚至更優的性能。針對測試函數1、2,如圖12、圖13所示,本文算法展現出針對復雜多峰問題近似建模的潛力,可以看出Rippa提出的FLOOCV法在處理多峰值問題時,沒有考慮樣本在設計空間內的分布,僅使用單一的形狀參數對所有基函數進行統一處理,其精度不夠理想。對于測試函數3,在處理角落處有尖峰出現的二維多極值問題時,Kitayama方法和Rippa方法均體現出一定的局限性。當測試函數存在大量極值時,例如測試函數4,Rippa方法的性能較差。

圖12 測試函數1均方根誤差箱型圖Fig.12 Boxplot of RMSE of test function 1

圖13 測試函數2均方根誤差箱型圖Fig.13 Boxplot of RMSE of test function 2

圖14 測試函數3均方根誤差箱型圖Fig.14 Boxplot of RMSE of test function 3

圖15 測試函數4均方根誤差箱型圖Fig.15 Boxplot of RMSE of test function 4

4 快速交叉驗證方法工程應用

針對工程仿真中常用的計算流體力學(Computational Fluid Dynamics, CFD)和有限元分析模型,以運載火箭減阻特性和和結構承載能力近似建模為例,驗證本文方法預測精度。

4.1 超聲速減阻桿阻力特性近似建模

針對火箭飛行過程中存在的激波阻力問題,常在火箭頭部安裝減阻桿以進行流場重構,減小氣動阻力。減阻桿作為被動減阻技術的一種,具有結構簡單、減阻效果好等優點,在實際工程系統中得到廣泛應用,如圖16、圖17所示的三叉戟潛射洲際彈道導彈及Igla單兵肩扛式防空導彈均采用減阻桿降低飛行器激波阻力。

圖16 三叉戟潛射洲際彈道導彈Fig.16 UGM-133A ballistic missile

圖17 Igla單兵肩扛式防空導彈Fig.17 Igla air defense missile

減阻桿的減阻效果受到包括減阻桿頂端半徑及減阻桿桿長等外形尺寸參數的影響,本節將采用FLUENT對圓盤形減阻桿與球形頭部結合減阻桿的減阻效果受到包括減阻桿頂端半徑及減阻桿桿長等外形尺寸參數的影響,本節將采用FLUENT對圓盤形減阻桿與球形頭部結合后外形的阻力系數進行數值計算,幾何外形如圖18 所示。設計參數為減阻桿頂端半徑和減阻桿桿長,見圖19。飛行器頭部半徑=100 mm,減阻盤厚度=10 mm。數值模擬中邊界條件設置為壓力遠場條件及壁面無滑移絕熱壁條件,來流馬赫數=2,氣體模型設置為完全氣體模型,氣流參數設置為標準海平面大氣參數,湍流模型選用SST-模型。由于該幾何外形為軸對稱外形,采用軸對稱條件進行計算可提高計算效率,因此選用軸對稱條件進行數值模擬。

圖18 超聲速飛行器外觀及減阻結構示意圖Fig.18 Diagram of supersonic vehicle appearance and drag reduction structure

圖19 減阻結構設計參數示意圖Fig.19 Diagram of drag reduction structure design parameters

和的取值范圍如表2所示,采用10×10個均勻分布的樣本點作為測試集,進行數值模擬得到100 組不同頂端半徑及桿長的減阻桿與球形頭部組合外形的阻力系數。運用本文所提出的方法進行近似建模,實驗設計方法生成20個點作為快速交叉驗證的數據集。選用快速5-折交叉驗證,每折樣本容量為4,在單次驗證過程中,16個為訓練數據,4個為測試數據。

表2 減阻結構設計參數Table 2 Drag reduction structure design parameters

4種方法分別計算得出的均方根誤差值如圖20 所示。在樣本點不均勻的2維氣動工程問題中,由于Kitayama法中形狀參數為一個固定值,不會根據實際情況進行調參,其所構建近似模型精度缺點明顯,誤差達到了FK-fold CV法的7倍左右。

圖20 工程算例1均方根誤差值Fig.20 RMSE of engineering example 1

4.2 重型運載火箭級間段加筋殼結構承載能力近似建模

采用復雜加筋殼結構設計問題對本文方法的有效性進行進一步驗證。運載火箭及導彈艙段等航天結構中常用的薄壁加筋柱殼結構,具有較高的軸壓、彎曲和扭轉承載效率的特點。作為運載火箭的主要承力部段,薄壁加筋柱殼的輕量化設計可以在提高其運載能力的同時降低成本。蒙皮桁條結構是典型框桁加強薄壁加筋柱殼結構之一,軸壓承載性能優越,廣泛應用于運載火箭結構設計。

圖21 加筋柱殼結構及中間框布局Fig.21 Diagram of cylindrical stiffened shell structure and intermediate frame

上述復雜結構設計工程問題具有19維設計變量,各參數變化范圍如表3所示,輸出響應為結構軸向承載極限能力。采用優化拉丁超立方得到設計空間內均勻分布的400個樣本點,作為訓練樣本來構建近似模型。在選擇縮放系數時,采用快速10-折交叉驗證,每折中應有40個數據。在單次交叉驗證過程中,360個樣本點為訓練數據,40個樣本點為測試數據。同樣采用優化拉丁超立方得到設計空間內均勻分布的5 000個樣本點,仿真得到的極限載荷作為真實輸出,計算出真實輸出與近似模型預測輸出的RMSE值,再對比本文方法與常用近似建模技術所得結果的誤差大小。

表3 蒙皮桁條結構設計參數Table 3 Skinned purlin structure design parameters

4種方法分別計算得出的均方根誤差值如圖22 所示。在高維結構問題中,Kitayama法與Kriging法都有一定的局限性,其誤差分別為FK-fold CV法的194%和158%。

圖22 工程算例2均方根誤差值Fig.22 RMSE of engineering example 2

5 結 論

1) 提出基于樣本局部密度的形狀參數表征方法,將多個形狀參數確定轉化為縮放系數確定的單變量優化問題,降低徑向基函數形狀參數確定的計算復雜度,解決現有RBF形狀參數確定中待求變量多和計算復雜等問題。基函數中形狀參數的合理確定大幅提升了RBF近似模型預測精度。采用4種方法分別對運載火箭進行氣動力和結構承載能力近似建模,結果表明采用本文方法在復雜工程問題上預測精度遠高于現有方法且相對穩定,通用性強且能滿足工程實際需求。

3) 由于Rippa提出的FLOOCV法是本文所提出的FK-fold CV法的特例,2種方法思路大致相同但本方法普適性更強。在針對上述算例進行對比時發現,這2種方法應用于低維測試函數時差異較為明顯,所以面向不同問題時,應對的取值進行討論。對于FK-fold CV算法本身而言,綜合時間與近似模型精度考慮,選取的最優值,是下一步值得探索的問題。