基于Stochastic Kriging模型的不確定性序貫試驗設計方法

王 波， GEA Haechang， 白俊強， 張玉東， 宮 建， 張衛民

(1. 中國航天空氣動力技術研究院 研發中心， 北京 100074； 2. 新澤西州立大學 機械宇航學院，新澤西 Piscataway， 08854； 3. 西北工業大學 航空學院， 陜西 西安 710072)

基于Stochastic Kriging模型的不確定性序貫試驗設計方法

王 波1， GEA Haechang2， 白俊強3， 張玉東1， 宮 建1， 張衛民1

(1. 中國航天空氣動力技術研究院 研發中心， 北京 100074； 2. 新澤西州立大學 機械宇航學院，新澤西 Piscataway， 08854； 3. 西北工業大學 航空學院， 陜西 西安 710072)

不確定性研究中需要計算大量重復樣本，這無疑對計算量較大的數值模擬提出了巨大的挑戰.通過試驗設計方法可以有效地減少不確定性研究中的計算量，然而，目前考慮不確定性的試驗設計方法研究大多仍專注于傳統試驗設計方法.針對這一問題，為了通過更為合理的計算資源分配得到更精準的不確定性評估，基于有限樣本的Stochastic Kriging模型提出了針對不確定性問題的三階段序貫試驗設計方法.首先，通過特定位置的采樣對IMSE進行簡化，構建了預選步進信息選取策略，通過預選增量樣本總個數以及各取樣位置處的分布信息，達到隨機代理模型目標精度要求；同時，基于IMSE構建了基于步進信息的單輪選點試驗設計準則，以同時考慮設計變量的取樣位置及其分布信息.由算例與傳統方法的對比分析可知，所建立方法通過等量的采樣得到了精度更高的隨機代理模型，驗證了其在不確定性問題中的可行性和優勢.

試驗設計方法； 不確定性； 代理模型； 均方差積分法； 序貫設計

在科學研究和工程技術中，幾乎所有的變量和參數都是非精確、非完全的，幾乎所有問題中都包含著不確定性，例如大氣狀態的擾動、機械的制造誤差、主觀判斷的經驗偏差等都是不確定性的表現.近年來，隨著計算機科學的飛速發展，為了使設計結果更具可靠性和穩健性，基于數值模擬的不確定性研究越來越多地引起各界關注[1-5].然而，在很多數值模擬的研究中，即使對控制方程單狀態的解算都需要消耗較多的計算資源，毫無疑問，對于需要大量重復采樣的不確定性研究而言，計算量更是難以承受.

由于代理模型可以通過有限次采樣得到整個空間可靠的響應信息，因此其在計算量較大的數值模擬研究中已被廣泛應用.目前應用較多的是確定性(deterministic)代理模型[6]，其可以通過已知的確定性樣本為輸入輸出建立映射關系，進而取代傳統的數值模擬，這類模型包括多項式響應面、支持向量機、Kriging和RBF等.而隨機(stochastic)代理模型是專門針對不確定性問題進行建模的.對于不確定性問題而言，由于隨機代理模型具有對問題針對性強、響應可靠等優勢[7]，在不確定性研究中受到了廣泛關注，其通過不確定性統計矩信息作為輸入建立不確定性輸入輸出間的映射關系，可直接對不確定性問題進行分析處理.由于不確定性研究仍屬于新興領域，因此現有隨機代理模型理論相對于確定性代理模型而言還不多，其中包括多項式混沌方法[8]、Stochastic Kriging[9]和Stochastic RBF[10].Kriging理論具有數學背景強、插值非線性程度高的特點，其作為確定性代理模型理論是非常具有吸引力的，而Stochastic Kriging是Kriging理論在隨機空間的拓展.因此，本文將基于Stochastic Kriging對不確定性的試驗設計方法展開研究.

試驗設計可以有效地減少不確定性研究中的計算量，其目的是通過合理采樣，充分地運用有限的計算資源.傳統的試驗設計方法[11]，比如拉丁超立方法[12]、均勻分布法等，由于在資源配置過程中具有不依賴現有的樣本信息、充滿設計空間等優勢，已被較多地用于不確定性初始采樣過程中.目前考慮不確定性的試驗設計方法研究大多仍專注于這類傳統的試驗設計方法，但對不確定性問題而言，理論上通過傳統試驗設計準則基本不可能得到最優的試驗分布[13-14].序貫試驗設計已被證明是比傳統試驗設計更行之有效的方法[15-16]，相比傳統試驗設計而言，序貫試驗設計在采樣過程中很大程度上利用了已有信息的價值，其不僅需要較少的計算資源，而且與真實試驗分階段、階次產生樣本過程比較接近.因此，對不確定性問題開展序貫設計的研究具有重要意義.同時值得注意的是，要進行序貫試驗設計，就需要基于已知信息額外建立序貫采樣的判斷標準[17-18].

在基于代理模型進行試驗設計的研究中，即使是很多確定性問題，其初始采樣都較為耗時，毫無疑問對于需要在各設計位置處進行重復取樣的不確定性問題更是如此.因此，基于針對不確定性問題的隨機代理模型開發高效的試驗設計方法，以對有限的計算資源進行合理分配，對不確定性問題的研究具有重要意義.在基于確定性Kriging模型的序貫試驗設計研究中，就判斷標準而言，大部分工作都是根據Kriging估計的均方差來開展的.其中具有代表性的包括最大均方差法(MMSE)和均方差積分法(IMSE)[17].IMSE是理論上直接表示整個預測模型精度的概念，因而在試驗設計中得到廣泛應用.

目前在基于隨機代理模型開展的試驗設計研究方面，大部分研究仍采用傳統的試驗設計方法，對于序貫試驗設計方法的研究還不多，其中有代表性的工作包括：Van等通過Bootstrapping對估計方差進行處理[14]，Ankenman等通過對IMSE優化問題的求解得到了采樣處重復樣本個數的估計[19].然而，這些方法還未能同時對不確定性問題的取樣位置和重復樣本個數進行處理.因此，本文在借鑒前人工作的基礎上，將基于作者前期建立的針對不確定性問題的有限樣本Stochastic Kriging模型[9]來開發序貫試驗設計方法，構建序貫試驗設計準則，通過模型初始化、預選樣本選取和單輪最終選點準則的構建，對取樣位置和各位置上的重復樣本點的個數進行選取.

本文將分為如下4個部分進行闡述：第1部分首先對Stochastic Kriging模型進行簡要敘述；第2部分將基于Stochastic Kriging模型提出分為3個階段的序貫試驗設計標準與方法；第3部分將通過基礎算例驗證所提出試驗設計方法的有效性；最后，對全文進行總結和展望.

1 Stochastic Kriging簡述

本小節將首先對Stochastic Kriging理論進行簡要闡述，詳細的推導可參見參考文獻[18].空間中的隨機樣本將通過如下形式表示：

(1)

各個位置處樣本的均值表示為

(2)

Stochastic Kriging方法采用最優線性估計的形式進行預測，以得到未知位置的均值：

(3)

(4)

而后通過一系列統計學推導，均方差MSE變為如下的形式：

(5)

隨后，Stochastic Kriging模型優化問題的形式可表示為

(6)

通過對這一優化問題的求解可得：

w=(CM+Cε)-1CM(x,·)，

(7)

w0=β-CM(x,·)T(CM+Cε)-1lnβ.

(8)

因此，均值和均方差也即：

(9)

(10)

2 序貫試驗設計方法

選取試驗設計標準是為了提高預測模型在整個設計域H內的估計精度，所以本節將基于Stochastic Kriging模型構建IMSE.該方法大致流程如下：

階段1:

1)通過拉丁超立方采樣(LHS)對設計變量的初始樣本{x1,x2,…,xn}進行選擇，獲得充滿設計域的樣本位置；

2)考慮隨機因素的影響，在初始樣本處進行n0次數值模擬，得到各初始位置n0個重復樣本；

3)基于所得采樣得到Stochastic Kriging模型不確定性相關性矩陣中的初始參數.

階段2:

1)對每次序貫試驗設計確定模型精度的目標水平r；

2)基于LHS獲取額外采樣位置；

3)通過額外樣本構建基于IMSE的試驗設計準則，以獲得增量樣本個數ΔN和采樣分布信息.

階段3:

1)建立基于IMSE的試驗設計標準，以同時對取樣位置和分布信息進行考慮；

2)運用粒子群優化算法進行尋優.

根據所得結果對Stochastic Kriging模型參數進行更新，并循環返回到第2階段的試驗設計中，直至達到所需精度.

2.1 樣本增量的選取

本部分將針對試驗設計階段2中獲取樣本增量個數的問題進行研究.如階段1所述，通過初始樣本選取，已經得到Stochastic Kriging模型中的參數.假設要在額外k個位置進行取樣，通過LHS選取X0={xn+1,xn+2,…,xn+k}.這時需要得到的是樣本增量數ΔN和樣本分布n={Nn+1,Nn+2,…,Nn+k}.為了通過對樣本分布n的選取獲得最高的模型精度，對優化問題的構建可以表示為

(11)

式中：Ni≥0(n+1≤i≤n+k)，

并且，

(12)

此后

(13)

為了對樣本增量的分配進行估計，同時因為相對于非固有不確定性而言，固有不確定性是較小量，而隨著采樣數量的增多及重復樣本個數的增加，固有不確定將趨于零，所以這里近似取C(n)≈CM.隨后，運用拉格朗日法對優化問題進行求解，ΔN的點數分布為

(14)

至此已得到相對于樣本增量總數在固定位置最佳重復樣本分配近似值，每個樣本增量都會有相應的最佳分布與其對應.所以為了得到目標的預測模型精度r，樣本增量的總個數可以通過公式IMSE≤r得到.

2.2 階段3的試驗設計

本部分將基于階段2的預選信息提出同時考慮采樣位置和重復采樣分布的試驗設計標準，將此標準轉化為優化問題可表示為

(15)

式中：Ni≥0( n+1≤i≤n+k)，

MSE(x0;X0;n)=CM(x0,x0;X0)-

隨后，通過推導IMSE的計算公式可以轉化為下面的形式：

(16)

從中可見，目標函數變量為待取樣位置X0和各位置上的取樣個數向量n.為獲得優化問題的解，此處將基于權重改進的混合PSO算法進行優化.該算法更新粒子位置和速度可表示為

(17)

式中：d=1，2，…，k；m為粒子個數，1≤i≤m；慣性因子是從0.8到0.2的遞減函數；加速因子c1=c2=2；pbestid為d維i粒子的最優值；gbestd為當前全局最優值.

3 算例分析

本部分將分別通過方差各處非均等和各處均等的排隊論問題來研究所建立的序貫試驗設計方法.所建立方法在預選步進信息時，需要對具體問題計算資源、代理模型的目標精度與步進信息權衡處理.例如，為了更精確地描述問題，對某實際問題設置了代理模型的目標精度，該問題隨即轉變為選定序貫設計多少次進行采樣較為合適.對于這一問題，可根據模型目標精度和現有精度之間的差值得到，同時依據步進信息由公式(14)得到所需步進的樣本數量.因此，本部分將重點對預選準則的有效性進行驗證，并對比分析試驗設計方法的最終結果.此外，為便于分析不確定性問題，所得結果精度將被近似到同一個數量級.

1)不確定性各處非均等的排隊論.

在穩定狀態下排隊論的等待時間Y的均值具有解析解，對于負荷x(0

(18)

此處，各位置取樣點是有限的，采樣位置為[0.1 0.3 0.5 0.7 0.9].在各位置處，響應值為一系列重復樣本，表1是基于排隊論的樣本信息，其中Y是真實響應值，Ymean是各位置根據20個重復樣本獲得的均值，由公式(2)計算得到，Variance是采樣位置的方差.

表1 排隊論算例

基于表1中的數據，可得到Stochastic Kriging模型的初始參數，同時可得Stochastic Kriging的初始IMSE=0.468 78.為了得到樣本增量的總個數ΔN，并對方法的有效性進行驗證，這里假定步進的目標精度r=0.33，即需要IMSE≤0.33.此時結合公式(13)和(14)，可得樣本增量數為20.其間，首先通過拉丁超立方選取額外的采樣位置X0={0.2， 0.4， 0.6， 0.8}，進而由公式(14)可知，最佳分布比例是353∶471∶565∶611.當增量數是20時，根據公式(13)可知，IMSE=0.322 76.這是約束條件下的最佳分布結果，即20個樣本將被用來作為增量數進行階段3的分析.

在階段3中，根據對公式(16)的優化，所得結果如表2中“最優”部分所示.一般來說，為了使模型更精確，對不確定性問題常采用均勻分布，等分樣本空間，或者根據直覺在不確定性大的位置多取樣本.所以表2同時給出了均勻分布和直覺分布的IMSE結果.通過分析可知，所建立方法得到的樣本分布為非均勻分布，并且在不確定性相對較小的0.19和不確定性相對較大的0.81位置上重復取樣個數較多.也即：要使模型更精確，不僅需要在不確定性較大處采樣，同時也要在已較精確處添加樣本.

對基于Stochastic Kriging進行的這一不確定性問題研究而言，IMSE數值結果主要來源于固有不確定性η(x0)，因為η(x0)本身是個較大的量.然而，隨著樣本個數的增加，在IMSE數值中這一部分將漸趨恒定，其不會對優化結果產生過多影響.

表2 算例1試驗設計結果對比

此外，在階段2預估樣本增量時，所得IMSE的結果為0.322 76，稍低于所建立方法優化所得結果.這是由于在預估步中，為得到更合理的樣本增量值，代入公式(13)的是樣本分配的小數值，分別為3.53，4.71，5.65，6.11，但階段3的優化中所采用的試驗設計準則與階段2不同，所得到的是整數化的結果.

2)不確定性各處均等的排隊論.

為了便于與不確定性各處非均等的問題對比研究，本算例采用各處不確定性均等的設置中，對表1中的數據只改變方差列，將其都設為0.1.

與算例1相似，首先可通過已知信息建立Stochastic Kriging模型，得到初始參數，這時初始IMSE=0.262 57.為對比算例1，同時保證階段2預估點的有效性，此處將增量樣本的個數設置為20.通過公式(13)和(14)可知，點數在X0={0.2， 0.4， 0.6， 0.8}處的分布n={4.90， 5.10， 5.10， 4.90}.將n代入(13)，可得IMSE=0.115 79.隨后，根據20個樣本增量個數開展階段3的分析，對公式(16)進行優化，所得結果如表3所示.

表3 算例2試驗設計結果

通過分析可知，與算例1類似，采樣都分布在偏離均勻分布的位置，且重復樣本分布也趨于在0.19和0.81位置上多取.與算例1不同的是，此問題的結果在0.19和0.81處擁有相同的分布.相對于原始模型的精度而言，兩者之間具有極大的可比性，顯示了預估步的實效性.

3)飛行器機翼氣動擾動的不確定性問題.

通常情況下對飛行器氣動問題的考量都是在給定設計狀態下進行的，屬于確定性問題.然而實際上，飛行器氣動問題常常由于各種不確定因素的擾動，確定性的考量結果偏離真實響應，這極可能使飛行器性能急劇變差.因此，目前在飛行器氣動問題的研究中，越來越多的研究工作開始考慮氣動擾動下的不確定性問題.

為了對所構建方法開展應用研究，選用DLR-F4的翼身組合體，采用結構網格，劃分的網格量為1.2×106個，DLR-F4翼身組合體表面網格如圖1所示.采用RANS方程、S-A湍流模型進行流場計算，單狀態計算時間為15 min.通過對Ma=0.75，Cl=0.50，Re=3×106的實驗結果進行對比檢驗，結果表明所采用流場計算方法結果較好地反映了實驗結果的變化趨勢，各典型剖面處的壓力分布對比如圖2所示，其中X/C為采用歸一化處理的相對弦長，Cp為壓力系數.

圖1 DLR-F4翼身組合體表面網格Fig.1 The surface grid of DLR-F4 wing-body configuration

圖2 DLR-F4翼身組合體計算和實驗結果對比Fig.2 The comparison between computation and experiment of DLR-F4 wing-body configuration

此處對機翼在馬赫數擾動情況下的不確定性問題進行研究，馬赫數的擾動服從正態分布N(0，0.1).與驗證算例類似，馬赫數各位置處采樣點個數是有限的，不確定性初始采樣個數為20個，采樣位置為{0.69 0.71 0.73 0.75 0.77}.據此可得模型的初始參數，同時可得模型初始IMSE=0.128 3×10-5.而后，通過設定步進的目標精度0.10×10-5，并結合公式(13)和(14)，得到樣本增量的總個數ΔN=16.在階段3中，將ΔN作為增量數進行分析，根據對公式(16)的優化，所得結果如表4所示.

表4 算例3試驗設計結果

4 結 語

本文基于Stochastic Kriging模型構建了適用于不確定性研究的序貫試驗設計方法，建立了預選點步進信息的選取策略和基于步進信息的單輪選點準則.通過算例分析驗證了所建立試驗設計方法的有效性和可行性.未來可對高維不確定性問題開展進一步的工作.

[1] ZANG T A， HEMSCH M J， HILBURGER M W. et al. Needs and opportunities for uncertainty-based multidisciplinary design methods for aerospace vehicles[R/OL]. [2016-11-02]. http:// www.cs.odu.edu/～mln/ltrs-pdfs/NASA-2002-tm211462.pdf.

[2] SLOTNICK J， KHODADOUST A, ALONSO J, et al. CFD vision 2030 study: a path to revolutionary computational aerosciences[R/OL]. [2016-11-02]. https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20140003093.pdf.

[3] BLATTNIG S R， LUCKRING J M， JOSEPH H M， et al. NASA standard for models and simulations: philosophy and requirements overview[J]. Journal of Aircraft, 2012， 50(1): 20-28．

[4] Editorial policy statement on numerical and experimental accuracy[J/OL].[2015-07-06]. http://servidor.demec.ufpr.br/CFD/bibliografia/erros_numericos/AIAA_Journals_NumericalAccuracy.pdf.

[5] MURTHY J Y， MATHUR S R. Computational heat transfer in complex systems: a review of needs and opportunities[J]. Journal of Heat Transfer, 2012， 134(3): 031016．

[6] RAZAVI S， TOLSON B A， BURN D H. Review of surrogate modeling in water resources[J]. Water Resources Research， 2012， 48(7): 107-116．

[7] XUE Z， MARCHI M， PARASHAR S， et al. Comparing uncertainty quantification with polynomial chaos and metamodels-based strategies for computationally expensive CAE simulations and optimization applications[R/OL]. [2016-11-02]. http://papers.sae.org/2015-01-0437/.

[8] GHANEM R， SPANOS P. Stochastic finite elements: a spectral approach[M]. New York: Courier Dover Publications, 2003．

[9] WANG Bo， BAI Jun-qiang. GEA Haechang. Stochastic kriging for random simulation metamodeling with finite sampling[C]. 39th ASME Design Automation Conference, Portland， Oregon, Aug. 5-8， 2013．

[10] VOLPI S， DIEZ M， GAUL N J， et al. Development and validation of a dynamic metamodel based on stochastic radial basis functions and uncertainty quantification[J]. Structural Multidisciplinary Optimization. 2015,51(2): 347-368.

[11] SANCHEZ S M. Work smarter， not harder: guidelines for designing simulation experiments[C]// Proceedings of the 2005 Winter Simulation Conference. Orlando, FLorida, Dec. 4-7, 2005: 69-82．

[12] KOEHLER J R， OWEN A B. Computer experiments[M]. Pennsylvania: Handbook of Statistics， 1996: 261-308．

[13] RIDGE E. KUDENKO D. Sequential experiment designs for screening and tuning parameters of stochastic heuristics[R/OL]. [2016-11-02]. http://www.imada.sdu.dk/～marco/EMAA/Papers/EMAA06-ridge.pdf.

[14] VAN Beers， KLEIJNEN JACK PC. Customized sequential designs for random simulation experiments: Kriging metamodeling and bootstrapping[J]. European Journal of Operation Research, 2008, 186(3): 1099-1113．

[15] PARK S， FOWLER J W， MACKULAK G T， et al. D-optimal sequential experiments for generating a simulation-based cycle time-throughput curve[J]. Operations Research， 2002， 50(6): 981-990．

[16] GHOSH B K， SEN P K. Handbook of sequential analysis[M]. New York: Marcel Dekker Inc.， 1991．

[17] SACKS J， WELCH W J， MITCHELL T J， et al. Design and analysis of computer experiments[J]. Statistical Science， 1989， 4(4): 409-423．

[18] WELCH W J， BUCK ITJ， Sacks J. Predicting and computer experiments[J]. Technometrics， 1992， 34(1): 15-25．

[19] ANKENMAN B E， NELSON B L， STAUM J. Stochastic kriging for simulation metamodeling[J]. Operations Research, 2010,58(2): 371-382．

The uncertainty-based sequential design of experiment method based on Stochastic Kriging metamodel

WANG Bo1， GEA Haechang2， BAI Jun-qiang3， ZHANG Yu-dong1，
GONG Jian1， ZHANG Wei-min1

(1. Research and Development Center， China Academy of Aerospace Aerodynamics， Beijing 100074, China; 2. Department of Mechanical and Aerospace Engineering， Rutgers， The State University of New Jersey， Piscataway NJ 08854; 3. School of Aeronautics， Northwestern Polytechnical University of China， Xi’an 710072， China)

The research on uncertainty requires many duplications and undoubtedly it puts forward a giant challenge to numerical simulations which is time-consuming. The amount of computation in the study of uncertainty can be effectively reduced through design of experiment method， but the current researches on design of experiment method about uncertainty mainly concentrate on traditional methods. Aiming at the problem, in order to address the problem and attain an accurate uncertainty assessment through reasonably allocating computational resources， the sequential design of experiment method with three stages was constructed based on the Stochastic Kriging metamodel with finite sampling. At the beginning， the criterion to choose the predetermined number and distribution of samples to attain certain accuracy of stochastic metamodel was proposed through the simplification of IMSE at specific sampling states. In addition， the criterion to obtain the optimum based on the predetermined information was also derived to simultaneously take the state and distribution of samples into account. Moreover， traditional methods were used to do the comparison with the proposed method， and the feasibility and advantages of proposed method were verified by examples with uncertainty， in which stochastic metamodel with more accuracy was achieved by using the same amount of sampling as traditional methods．

design of experiment method; uncertainty; metamodel; integration of mean square error; sequential design

2015-07-09.

本刊網址·在線期刊：http://www.zjujournals.com/gcsjxb

國家自然科學青年基金資助項目(11302213).

王波(1984—),男,河南內黃人,博士,從事神經網絡、隨機建模、飛行器設計等研究，E-mail:alexanbo@163.com. http://orcid.org//0000-0002-9913-5348

10.3785/j.issn. 1006-754X.2016.06.002

TP 391.9

A

1006-754X(2016)06-0530-07