晏 良，段曉君，劉博文，徐 琎

(國防科技大學 文理學院， 湖南 長沙 410073)

在工程研究中，利用計算機通過仿真試驗可對真實復雜系統(tǒng)進行高精度的模擬，從而分析該系統(tǒng)的相關(guān)特性[1]。然而，單次高精度仿真試驗通常面臨計算效率低、時間成本大等問題，因此諸如可靠性評估、失效概率驗證等需要大規(guī)模仿真的任務在實際操作中很難完成。代理模型由于其計算復雜度低、數(shù)學表達明晰等特點，在計算機仿真試驗中得到了廣泛應用[2-3]。實際上，代理模型重點關(guān)注輸入和輸出之間的數(shù)學關(guān)系。給定試驗設計及觀測數(shù)據(jù)，代理模型通過回歸或插值等方法擬合觀測數(shù)據(jù)，可給出相應的預測均值、方差等統(tǒng)計量。

常用的代理模型構(gòu)造方法有很多，例如：多項式響應曲面(Polynomial Response Surface, PRS)[4]、多元自適應回歸樣條(Multivariate Adaptive Regression Splines, MARS)[5]、高斯過程(Gaussian Process, GP)[6]，等等。針對不同的模型特點(如光滑、線性、非線性)和數(shù)據(jù)特征(噪聲、異常值)，通常會選用不同的代理模型。

事實上，除了以上文獻所提供的經(jīng)驗，基于一些準則可以對代理模型進行篩選。經(jīng)典的模型選擇準則[7]包括：基于擬合優(yōu)度(Goodness of Fit, GoF)、基于預測誤差(Prediction Error, PE)或者模型誤差(Model Error, ME)、基于概率分布距離(Distributional Discrepancy, DD)、基于后驗概率(Posterior Probability, PP)等。然而，基于上述準則需要對每類候選模型進行計算，會增加較多的計算量；另外，從所有候選模型中僅挑選單個最優(yōu)模型，不考慮其他模型也在一定程度上造成了信息的浪費。出于上述考慮，一些研究[8-10]也側(cè)重于利用所有的候選模型，即構(gòu)造加權(quán)代理模型。

此外，構(gòu)造加權(quán)代理模型有利于減小不確定性。加權(quán)模型由于考慮了所有不同類型的代理模型，因此在一定程度上增加了該模型的穩(wěn)健性。同時，加權(quán)模型在構(gòu)造的同時充分考慮到各模型的優(yōu)勢，因此其精度也能得到保證。

經(jīng)典的加權(quán)模型方法往往從兩個指標考慮：一是考慮模型預測的方差大小，通常情況下，預測方差越小，則認為模型越準確；二是考慮模型的預測精度，即預測均值與真實值之間的偏差，加權(quán)的主要目的在于降低加權(quán)模型與真實模型之間的偏差。然而在實際應用中，預測精度高(即偏差小)的代理模型，其預測方差可能較大，或者預測方差小的代理模型，其與真實值之間的偏差比較大。本文基于以上問題，提出了基于KL(Kullback-Leibler)距離離散度的加權(quán)方法，直接對預測的分布入手，同時考慮預測的均值和方差，在兼顧預測精度的同時，降低預測模型與真實模型之間的分布差異。本文考慮兩類子代理模型，一類為線性回歸模型，另一類為高斯回歸模型。

1 代理模型簡介

1.1 線性回歸模型

線性回歸模型假設響應y是一類基函數(shù)(basis)的線性組合：

y=[f(x)]Tβ+

(1)

其中，x=(x1,…,xd)T為d維變量，f(x)為p維向量函數(shù)，其中每一項是關(guān)于x的基函數(shù)。通常在一類模型中基函數(shù)是給定的，因此構(gòu)造線性回歸模型重點在于求解未知參數(shù)β。一般假設為零均值的隨機誤差，因此[f(x)]Tβ實際上表示響應y的期望值。

給定N個觀測(Xi,Yi),i=1,…,N，即試驗設計(X,Y)，根據(jù)式(1)可以得到：

Y=Fβ+

(2)

其中，F(xiàn)∈RN×p為設計矩陣，=(1,…,N)T為觀測誤差。由于～N(0,σ2I)，可得參數(shù)β的最優(yōu)線性無偏估計(Best Linear Unbiased Estimate, BLUE)為：

(3)

(4)

(5)

1.2 高斯回歸模型

高斯回歸模型是一類典型的非線性模型。與線性回歸模型指定基函數(shù)不同，高斯回歸模型建立在響應y服從某一個確定的高斯過程基礎之上，即：y～GP(0,k(x,x)+σ2)，其中k(x,x)為核函數(shù)。具體來說，給定試驗設計(X,Y)，則有如下似然函數(shù)：

Y|X～N(0,k(X,X)+σ2I)

(6)

記K=k(X,X)∈RN×N，由高斯過程假設以及式(6)，可得預測均值和方差分別為：

(7)

(8)

其中，K*=k(X,x)∈RN×1,K**=k(x,x)∈R。實際上，高斯過程的預測能力由核函數(shù)所決定。一般來說，針對似然函數(shù)式(6)進行優(yōu)化可以得到關(guān)于核函數(shù)參數(shù)的相關(guān)估計。

2 加權(quán)代理模型

不同的代理模型其性能也有所不同，因此在加權(quán)代理模型中的權(quán)重也會隨之改變。加權(quán)代理模型的一般表達式為：

(9)

2.1 基于預測誤差的權(quán)函數(shù)構(gòu)造

(10)

由此可求解權(quán)函數(shù)如下：

(11)

2.2 基于交叉驗證均方誤差和均方根誤差的權(quán)函數(shù)構(gòu)造

Goel[9]等選擇泛化交叉驗證均方誤差(Generalized Mean Square cross validation Error, GMSE)作為衡量各子代理模型精度的重要參數(shù)，構(gòu)造權(quán)函數(shù)如下：

(12)

其中，α<1、β<0為引入的尺度變換參數(shù)，Gi為第i個子代理模型的GMSE，即：

(13)

(14)

GMSE反映了代理模型在所對應試驗設計點上的平均誤差，GMSE越小，說明模型在該試驗設計上越穩(wěn)定，即模型與試驗設計相適應；而RMSE反映了代理模型在測試集上的預測與真實值的誤差，與GMSE相比，其在一定程度上更關(guān)注模型在全局上的預測性能，RMSE越小，其預測精度越高。同時注意到，與式(11)不同，基于GMSE(RMSE)的權(quán)函數(shù)在模型定義域上是常值，因此權(quán)函數(shù)代表各子模型整體預測性能的比重。

此外，由于GMSE對每個子模型都需要計算N次，因此會增加較大的計算量。一般情況下，可以進行并行處理或者進行k(≤10)折交叉驗證。然而，針對線性回歸模型，GMSE的計算可以直接進行簡化，參見文獻[7]第12章。對于RMSE來說，則需要增加Nv個觀測，對于復雜系統(tǒng)來說，這會成為主要制約因素。

2.3 基于優(yōu)化加權(quán)模型誤差的權(quán)函數(shù)構(gòu)造

Acar[10]等在前述誤差的基礎上，提出了基于優(yōu)化誤差的權(quán)函數(shù)構(gòu)造方法，即求解如下最優(yōu)化問題：

2012年黨的十八大以來，以習近平同志為核心的黨中央堅定不移地堅持高舉改革開放旗幟、推進全面深化改革、擴大對外開放的戰(zhàn)略決策部署，展現(xiàn)了黨中央將改革開放進行到底的政治魄力和堅定決心。

(15)

其中，Err[·]表示誤差函數(shù)，如RMSE或GMSE。對于式(15)的求解分為兩步，首先針對試驗設計構(gòu)造各子代理模型，其次將各子代理模型代入式(15)求解該約束優(yōu)化問題。與基于PE(RMSE, GMSE) 的方法相比，其計算量主要增加在求解最優(yōu)化問題式(15)，當子模型數(shù)量M較大時，很難求得全局最優(yōu)解。

3 Kullback-Leibler距離離散度加權(quán)方法

以上三類經(jīng)典加權(quán)模型方法基于誤差函數(shù)，往往只考慮代理模型預測的某一方面，如預測方差或偏差。本文從代理模型的預測分布入手，同時考慮這兩類因素，利用Kullback-Leibler距離[11]構(gòu)造權(quán)函數(shù)，從而降低預測總體分布的不確定性。

與式(9)直接假設加權(quán)模型為各模型預測值的加權(quán)不同，本文假設加權(quán)模型在各點的預測概率密度函數(shù)為各子代理模型預測概率密度函數(shù)的加權(quán)，即：

(16)

(17)

(18)

(19)

其中，DKL為通常意義下的Kullback-Leibler距離。對于本文中的兩類模型，其在各點的分布服從高斯分布，同時各模型之間相互獨立，因此概率密度函數(shù)p=N(μp,σp)和q=N(μq,σq)的Kullback-Leibler距離可以計算如下：

(20)

從式(20)可以看出，該距離包含預測均值和方差兩種參數(shù)。基于距離的權(quán)重函數(shù)構(gòu)造方法有多種，如基于反函數(shù)和減函數(shù)等特殊函數(shù)。本文擬采用高斯核函數(shù)，基于DSKL可定義權(quán)函數(shù)如下：

(21)

(22)

此外，尺度參數(shù)ρ可以選取為全體分布的平均離散度，即各分布之間的平均距離為：

(23)

將式(22)和式(23)代入式(21)中，并對其進行歸一化可以得到最終的權(quán)函數(shù)如下：

(24)

4 算例分析

本節(jié)首先針對兩個測試函數(shù)對基于Kullback-Leibler距離離散度的加權(quán)代理模型進行分析，其函數(shù)表達式如下：

f1(x)=(10cos2x+15-5x+x2)/50

(25)

(26)

其中，f1是Viana函數(shù)，f2為Rastrigin函數(shù)，且觀測誤差均為～N(0,0.1)。

其次，本文考慮5類子代理模型及其加權(quán)模型，其中2類屬于線性回歸模型，3類屬于高斯回歸模型。子代理模型的具體形式如表1所示。

因此，給定試驗設計(X,Y)，本文比較了9類代理模型，即5類子代理模型M1～M5(如表1所示)；基于PE的加權(quán)模型M6；基于GMSE的加權(quán)模型M7；基于優(yōu)化GMSE的加權(quán)模型M8；本文方法M9。同時，對于每個算例，均隨機生成包含10 000個樣本點的測試集(Xv,Yv)，并基于測試集比較兩種指標：RMSE和經(jīng)驗累積分布函數(shù)(Empirical Cumulative Distribution Function,ECDF)。具體計算結(jié)果如下。

表1 5類子代理模型

4.1 Viana函數(shù)

假設x～N(0,1)，則利用拉丁超立方體抽樣方法(Latin Hypercube Sampling, LHS)生成包含7個樣本的試驗設計(X,Y)。對每類代理模型總共進行50次獨立試驗，則其RMSE對比如圖1所示。

圖1 9類代理模型RMSE對比Fig.1 Comparison of RMSE of 9 surrogate models

由圖1可以看出，對于5類子代理模型來說，M4的逼近精度最高；對于4類加權(quán)模型來說，本文方法，即M9的逼近精度也是最高的；總體上，本文方法的精度與M4精度相近。圖2利用1000個樣本的測試集生成經(jīng)驗累積分布函數(shù)。

圖2中可以看到，對比于各子模型，加權(quán)模型的ECDF都能夠較好地反映真實響應的分布情形。同樣地，可以看出M9與真實響應之間的差距更小，在y>5的區(qū)域上基本與真實ECDF重疊。

圖2 9類代理模型ECDF對比Fig.2 Comparison of ECDF of 9 surrogate models

4.2 Rastrigin函數(shù)

假設xi～N(0,1),i=1,2，同樣利用拉丁超立方體抽樣方法(Latin Hypercube Sampling, LHS)生成包含20個樣本的試驗設計(X,Y)，并對每類代理模型總共進行50次獨立試驗，則其RMSE對比如圖3所示。

由圖3可以看出，對于5類子代理模型來說，除了M1以外，M2～M5逼近精度類似，其中M2,M3的整體預測性能更加穩(wěn)定；對于4類加權(quán)模型來說，其預測性能與M2～M5相當，M9的逼近精度略高，在某些試驗設計條件下具有最高的精度。圖4與圖2類似，同樣利用1000個樣本的測試集生成經(jīng)驗累積分布函數(shù)，其結(jié)果如圖4所示。

圖3 9類代理模型RMSE對比Fig.3 Comparison of RMSE of 9 surrogate models

圖4 9類代理模型ECDF對比Fig.4 Comparison of ECDF of 9 surrogate models

由圖4可以看出，除了M1,M9以外，其他代理模型對真實響應分布逼近并不理想。同時注意到，對這些模型，其數(shù)據(jù)分布集中在y=8.4附近；作為對比，M9仍然能夠較好地還原真實響應的分布。

最后對一個具體的案例進行測試，即懸臂梁(cantilever beam)問題[12]，如圖5所示。

圖5 懸梁臂問題Fig.5 Cantilever beam problem

其中，w和t分別代表橫截面的寬度和厚度，Y表示垂直載荷，X表示水平載荷，L代表長度。對于懸臂梁問題來說，總共有兩個響應，即位移D(x)和應力S(x)。具體數(shù)學表達式如下：

(27)

式中總共有三個自變量，即楊氏模量(Young′s modulus)E～N(2.9×107, 1.45×106)，垂直載荷Y～N(1000,100)，水平載荷X～N(500,100)。令w=4,t=2，同樣利用LHS方法生成包含36個樣本的試驗設計，并對9類代理模型總共進行50次獨立試驗，則對于應力響應來說，其RMSE對比如圖6所示。

由圖6可以看出，對于5類子代理模型來說，M4、M5的預測精度較低，其RMSE遠遠大于200，因此未在圖中展示。對于加權(quán)模型來說，M6、M7的效果同樣較差，M8、M9的逼近性能相當，與最佳子模型M2的精度也在同一量級。對于ECDF來說，同樣利用1000個樣本的測試集生成經(jīng)驗累積分布函數(shù)，其結(jié)果如圖7所示。

圖6 應力RMSE對比Fig.6 Comparison of RMSE of stress

事實上，由圖7可以看出，子模型M4、M5對本問題及試驗設計的應用效果較差，而M1～M3均能夠較好地估計應力的分布。對于加權(quán)模型來說，其對各子模型的“平均”作用都能使得模型更加穩(wěn)健，如M6、M7，一定程度上還原了真實分布，但在性能上有著較大差異。M8、M9能夠極大地利用“好”模型(M1～M3)的信息，從而其對ECDF的逼近效果也是最好的。對于位移D(x)來說，由于M3不能應用在該模型之上，因此本文僅比較其余8類代理模型。注意到此時加權(quán)模型為M1、M2、M4、M5四類模型的加權(quán)。使用與應力相同的試驗設計，得到關(guān)于位移的RMSE如圖8所示。

圖7 應力ECDF對比Fig.7 Comparison of ECDF of stress

由圖8可以看出，對于4類子代理模型來說，M4～M5其預測精度同樣較低，因此也未在圖中展示。對于加權(quán)模型來說，M6的效果最差，M7、M8、M9的逼近性能依次遞增，一般情況下，M9的效果要比單獨的代理模型都要好。利用1000個測試樣本生成相關(guān)的ECDF，其對比如圖9所示。

圖9與圖7類似，對比于各子模型，加權(quán)模型至少避免了類似M4、M5的模型選擇錯誤，其在一定程度上能夠較為穩(wěn)健地逼近原模型。在所有加權(quán)類型中，本文方法與M8對ECDF的逼近是最準確的。

圖8 位移RMSE對比圖Fig.8 Comparison of RMSE of Displacement

圖9 位移ECDF對比Fig.9 Comparison of ECDF of displacement

5 結(jié)論

與通常加權(quán)模型選擇降低加權(quán)模型精度不同，本文針對預測分布構(gòu)造權(quán)函數(shù)，引入Kullback-Leibler距離刻畫模型與模型之間預測分布的差異性。為了克服真實分布無法獲得的問題，本文利用各模型之間的離散度來表征各模型的預測能力。算例分析表明，本文方法在兼顧模型預測精度的同時，也能夠較好地還原真實響應的分布。