崔慶安

(鄭州大學(xué)管理工程學(xué)院，鄭州 450001)

0 引言

對(duì)于產(chǎn)品實(shí)現(xiàn)復(fù)雜且成本較高的自動(dòng)化加工制造過程，計(jì)算機(jī)仿真是實(shí)現(xiàn)產(chǎn)品優(yōu)化設(shè)計(jì)的主要手段之一。其基本原理是根據(jù)物理或化學(xué)的先驗(yàn)知識(shí)，建立實(shí)際制造過程的數(shù)學(xué)模型，而后利用有限元分析等數(shù)值計(jì)算方法，通過運(yùn)行計(jì)算機(jī)程序代碼對(duì)制造過程進(jìn)行仿真模擬［1］。由于計(jì)算機(jī)仿真模型的計(jì)算復(fù)雜度高，計(jì)算時(shí)間長，為提高優(yōu)化設(shè)計(jì)的效率，研究者引入了計(jì)算機(jī)實(shí)驗(yàn)(Computer experiments，CE)的方法［2］。即采用實(shí)驗(yàn)設(shè)計(jì)(Design of Experiments，DOE)的方法［3］，對(duì)計(jì)算機(jī)仿真過程的輸入輸出進(jìn)行采樣，擬合出該過程的統(tǒng)計(jì)學(xué)代理模型，并通過對(duì)代理模型的尋優(yōu)來實(shí)現(xiàn)優(yōu)化設(shè)計(jì)。計(jì)算機(jī)實(shí)驗(yàn)已被廣泛地應(yīng)用于汽車制造、機(jī)械電子加工、生物醫(yī)藥合成、航空航天設(shè)計(jì)等領(lǐng)域［4-5］，取得了顯著的成效。

如何得到預(yù)測(cè)性能較好的統(tǒng)計(jì)代理模型，是計(jì)算機(jī)實(shí)驗(yàn)的關(guān)鍵問題。為提高模型對(duì)于復(fù)雜過程的全局性描述能力，多采用 Kriging［6］、ANN［7］、SVR［8］、B-樣條函數(shù)等結(jié)構(gòu)較為靈活的非參數(shù)模型作為基本代理模型，其中應(yīng)用最廣泛的是Kriging模型。此方面的研究較成熟，眾多研究者對(duì)于代理模型的性能對(duì)比、Kriging模型的優(yōu)勢(shì)等等進(jìn)行了深入研究［9-10］。模型形式確定后，實(shí)驗(yàn)設(shè)計(jì)方式的選擇就成為關(guān)鍵。目前常用的是以超拉丁方抽樣［11］(Latin hypercube sampling，LHS)、均勻設(shè)計(jì)［12］(Uniform design，UD)為代表的空間填充設(shè)計(jì)。LHS與UD的實(shí)驗(yàn)點(diǎn)個(gè)數(shù)與因子個(gè)數(shù)無關(guān)，并且實(shí)驗(yàn)點(diǎn)在各因子維度上的投影是均勻和不重疊的，因此可以將少量的實(shí)驗(yàn)點(diǎn)均勻地散布在可行域內(nèi)。但是當(dāng)因子個(gè)數(shù)較多時(shí)，少量的實(shí)驗(yàn)點(diǎn)對(duì)于可行域空間的填充程度不夠，而且極有可能出現(xiàn)實(shí)驗(yàn)點(diǎn)較為稀疏的子區(qū)域，極大影響了代理模型的性能。為提高代理模型的全局性預(yù)測(cè)性能，不得不采用較大的樣本量。對(duì)于空間填充設(shè)計(jì)來說，如何實(shí)現(xiàn)樣本量與模型的全局性預(yù)測(cè)性能的之間平衡，有價(jià)值的研究尚不多見。

事實(shí)上，對(duì)于計(jì)算機(jī)仿真過程來說，由于輸入因子與輸出特性之間存在著復(fù)雜的非線性關(guān)系，因而輸出特性在可行域內(nèi)的變化是不均勻的，一味地追求實(shí)驗(yàn)點(diǎn)的均勻分布，并不一定能夠反映這種變化特征。如果能夠根據(jù)輸出特性的變化情況，有針對(duì)性地調(diào)整空間填充設(shè)計(jì)實(shí)驗(yàn)點(diǎn)的分布，在輸出特性變化較大的空間安排較多的實(shí)驗(yàn)點(diǎn)，而在輸出特性變化較小的空間安排較少的實(shí)驗(yàn)點(diǎn)，將會(huì)在降低實(shí)驗(yàn)設(shè)計(jì)樣本量的同時(shí)提高代理模型的預(yù)測(cè)性能。

本文提出一種基于正交設(shè)計(jì)和LHS的兩階段實(shí)驗(yàn)設(shè)計(jì)策略，首先以正交設(shè)計(jì)初步探知輸出特性方差在可行域內(nèi)的變化規(guī)律，根據(jù)其變化程度，識(shí)別出輸出特性變化較大的子區(qū)域，而后在此子區(qū)域內(nèi)采用LHS來安排實(shí)驗(yàn)點(diǎn)，再利用Kriging模型建立起過程的全局性模型。

1 理論簡(jiǎn)介

1.1 計(jì)算機(jī)實(shí)驗(yàn)的基本原理

設(shè)計(jì)算機(jī)仿真過程的輸入因子為x=［x1，…xm］'，其可行域?yàn)閄∈Rm，輸出特性為 y，且 y與 x之間存在某種確定型關(guān)系，即：

式(1)通常高度復(fù)雜，不存在明確的解析表達(dá)式，一般采用有限元分析等數(shù)值計(jì)算的方法來表達(dá)。這樣一來，當(dāng)仿真模型的精度較高時(shí)，程序運(yùn)行需要耗費(fèi)大量時(shí)間。例如對(duì)于汽車車身結(jié)構(gòu)的優(yōu)化設(shè)計(jì)，用于模擬車身正面碰撞過程的有限元模型的單元可達(dá)數(shù)萬，完成一次車身結(jié)構(gòu)的優(yōu)化需要數(shù)百小時(shí)。此時(shí)，建立計(jì)算機(jī)仿真過程的統(tǒng)計(jì)學(xué)代理模型以簡(jiǎn)化計(jì)算就顯得十分有必要。即首先選擇某種形式的代理模型(例如Kriging模型)，而后根據(jù)某種準(zhǔn)則生成實(shí)驗(yàn)設(shè)計(jì)集

分別以x1，…，xn為輸入，運(yùn)行計(jì)算機(jī)仿真過程，獲取對(duì)應(yīng)的輸出特性值y1，…，yn，形成樣本集：

再根據(jù)該樣本集擬合出式(1)的具體統(tǒng)計(jì)學(xué)代理模型：

由于g(x)一般較為簡(jiǎn)單，且具有明確的解析表達(dá)式，因此用其來近似的代替計(jì)算機(jī)仿真過程將會(huì)顯著地縮短計(jì)算時(shí)間，而且可以通過對(duì)g(x)的尋優(yōu)近似地實(shí)現(xiàn)產(chǎn)品的優(yōu)化設(shè)計(jì)。為保證優(yōu)化結(jié)果的精度，希望在可行域X內(nèi)，g(x)對(duì)于f(x)的預(yù)測(cè)誤差越小越好，即：

1.2 超拉丁方抽樣

與傳統(tǒng)的因子設(shè)計(jì)、中心復(fù)合設(shè)計(jì)相比，空間填充設(shè)計(jì)的實(shí)驗(yàn)點(diǎn)的安排沒有特定的位置要求，只需滿足某些設(shè)計(jì)準(zhǔn)則即可。其實(shí)驗(yàn)點(diǎn)安排較為靈活，而且實(shí)驗(yàn)次數(shù)n與參數(shù)個(gè)數(shù)m無關(guān)，不會(huì)出現(xiàn)所導(dǎo)致的“維數(shù)災(zāi)難”，較適用于復(fù)雜作用關(guān)系過程的實(shí)驗(yàn)設(shè)計(jì)與建模，代表性的方法包括LHS、均勻設(shè)計(jì)等。

假設(shè)輸入因子的個(gè)數(shù)為m，實(shí)驗(yàn)次數(shù)為n(n＞m)，LHS首先將可行域均勻分成nm個(gè)超立方體的子區(qū)域，選取其中的n個(gè)子區(qū)域，使所有子區(qū)域都能夠不重疊地投影在每一個(gè)維度上;而后在各個(gè)子區(qū)域內(nèi)隨機(jī)選取一個(gè)點(diǎn)組成n個(gè)實(shí)驗(yàn)點(diǎn)。LHS具有操作簡(jiǎn)單、樣本量小、不易導(dǎo)致維數(shù)災(zāi)難等優(yōu)點(diǎn)。但LHS是一種分層隨機(jī)抽樣方法，在某些情況下，生成實(shí)驗(yàn)點(diǎn)的隨機(jī)性與均勻性不高，實(shí)驗(yàn)設(shè)計(jì)的效果有待提高。

2 方法研究

2.1 基本思路

由上述分析可知，LHS的實(shí)驗(yàn)點(diǎn)的個(gè)數(shù)與在各因子維度上的投影較為均勻，但這并不意味著該實(shí)驗(yàn)設(shè)計(jì)方式對(duì)于復(fù)雜的計(jì)算機(jī)仿真過程具有良好代表性。首先，對(duì)于因子個(gè)數(shù)為m的LHS設(shè)計(jì)，若在nm個(gè)子區(qū)域中隨機(jī)選擇n個(gè)子區(qū)域安排實(shí)驗(yàn)點(diǎn)，則共有nm－n個(gè)子區(qū)域未被填充。當(dāng)m較大時(shí)，有nm－n?n，這樣一來，雖然樣本點(diǎn)較為均勻分散，但是相對(duì)稀疏，仍有大量的未安排實(shí)驗(yàn)點(diǎn)的子區(qū)域。其次，由于計(jì)算機(jī)仿真過程的輸出特性與輸入因子之間存在復(fù)雜非線性作用關(guān)系，因而輸出特性在整個(gè)可行域內(nèi)的分布是不均勻的，既有輸出特性的變化較為平坦的子區(qū)域，也有變化較大的子區(qū)域，還有存在局部極值的子區(qū)域。顯然，實(shí)驗(yàn)點(diǎn)的均勻分散并不能夠刻畫輸出特性的不均勻變化，樣本點(diǎn)的利用效率不高。

事實(shí)上，為提高實(shí)驗(yàn)設(shè)計(jì)的效率，應(yīng)該根據(jù)子區(qū)域內(nèi)輸出特性的變化情況有針對(duì)性地安排實(shí)驗(yàn)點(diǎn)的數(shù)目，對(duì)于輸出特性變化較大的子區(qū)域多安排實(shí)驗(yàn)點(diǎn)，而對(duì)于輸出特性變化較小的子區(qū)域少安排實(shí)驗(yàn)點(diǎn)。進(jìn)一步地，在沒有過程先驗(yàn)知識(shí)的情況下，可以采用對(duì)稱(即各因子水平數(shù)均相等)的正交設(shè)計(jì)來初步探知輸出特性在可行域內(nèi)的變化情況。根據(jù)正交設(shè)計(jì)的理論，對(duì)于一個(gè)因子數(shù)為m，各因子水平數(shù)均為q，實(shí)驗(yàn)次數(shù)為n的對(duì)稱正交設(shè)計(jì)，有n=qm－p，1≤p≤m－2，因此實(shí)驗(yàn)點(diǎn)在每一個(gè)因子的每一個(gè)水平上的投影均有重復(fù)，且重復(fù)次數(shù)為n/q。那么根據(jù)這n/q次重復(fù)即可估計(jì)出在該水平處過程輸出特性的方差，而方差變化較大的水平處則提示該因子水平附近過程輸出特性的變化較大。據(jù)此識(shí)別出過程輸出特性變化較大的子區(qū)域，在該子區(qū)域內(nèi)采用樣本量較小的LHS實(shí)驗(yàn)設(shè)計(jì)，即可在降低樣本量同時(shí)，提高實(shí)驗(yàn)設(shè)計(jì)效率及模型性能。

2.2 實(shí)現(xiàn)步驟

根據(jù)上述分析，這里提出基于正交設(shè)計(jì)與LHS的計(jì)算機(jī)實(shí)驗(yàn)設(shè)計(jì)及建模方法如下：

Step1：根據(jù)輸入因子的個(gè)數(shù)m及問題規(guī)模，確定正交設(shè)計(jì)的各因子的水平數(shù)q，利用正交表，選擇合適的正交設(shè)計(jì)集POD：

Step2：根據(jù)正交設(shè)計(jì)POD，運(yùn)行計(jì)算機(jī)仿真代碼，獲得初始樣本集：

Step4：將各因子維度上過程輸出特性變化較大的區(qū)間進(jìn)行組合，確定需要進(jìn)行LHS設(shè)計(jì)的子區(qū)域X*;

Step5：在子區(qū)域內(nèi)X*安排適當(dāng)數(shù)目的LHS實(shí)驗(yàn)點(diǎn)，運(yùn)行計(jì)算機(jī)仿真代碼，獲取補(bǔ)充樣本集：

Step6：根據(jù)樣本集SOD∪S，擬合計(jì)算機(jī)仿真過程的Kriging模型。

3 算例研究

3.1 仿真算例

為考察方法的性能，這里給出了一個(gè)較為復(fù)雜的算例。設(shè)計(jì)算機(jī)仿真過程的輸入輸出之間存在如下關(guān)系：

圖1給出了在x3= －3，0，3時(shí)，y關(guān)于x1，x2的等高線圖，可以看出，這是一個(gè)典型的非線性復(fù)雜作用關(guān)系過程。

圖1 過程的等高線圖

為說明方法的性能，選擇與實(shí)驗(yàn)設(shè)計(jì)樣本集不同的測(cè)試樣本集：

令yi代表測(cè)試樣本集的實(shí)際輸出，^yi代表Kriging模型對(duì)于測(cè)試樣本集的預(yù)測(cè)輸出，則預(yù)測(cè)誤差可表示為：

采用如下指標(biāo)評(píng)價(jià)模型的性能：

(1)預(yù)測(cè)的均方誤差

(2)預(yù)測(cè)誤差的標(biāo)準(zhǔn)差

(3)預(yù)測(cè)的最大絕對(duì)偏差：

3.2 計(jì)算機(jī)實(shí)驗(yàn)的結(jié)果

按照3.2所示步驟，首先選擇L50(51121)正交表的前3列作為正交設(shè)計(jì)的實(shí)驗(yàn)點(diǎn)集，如表1所示。表1還給出了正交實(shí)驗(yàn)點(diǎn)及對(duì)應(yīng)的輸出特性值。圖2給出了各因子維度上5個(gè)不同水平處的過程輸出特性的樣本方差，在x1、x3的維度上各有一個(gè)局部極值，而在x2的維度上有兩個(gè)局部極值。進(jìn)一步，可以根據(jù)式(8)確定了各因子維度的包含過程輸出特性方差極值的區(qū)間，分別為：

表1 第一階段的正交設(shè)計(jì)實(shí)驗(yàn)點(diǎn)及輸出特性值

圖2 各因子水平處的輸出特性的樣本方差

圖3 兩階段設(shè)計(jì)與傳統(tǒng)LHS的實(shí)驗(yàn)點(diǎn)對(duì)比(樣本量 =80)

而后在上述兩子區(qū)域內(nèi)添加樣本量大小相同的LHS實(shí)驗(yàn)點(diǎn)。為說明方法的性能，選取了不同樣本量大小(n2=10+10，n2=15+15，n2=20+20)的LHS設(shè)計(jì)，與原有的正交設(shè)計(jì)組成樣本集，而后采用Kriging方法進(jìn)行建模。作為對(duì)比，單純地采用LHS設(shè)計(jì)生成樣本量相同的全部實(shí)驗(yàn)點(diǎn)集進(jìn)行Kriging建模。在可行域內(nèi)隨機(jī)生成5000個(gè)測(cè)試樣本點(diǎn)，利用式(14)至式(16)考察各類模型的預(yù)測(cè)性能，其對(duì)比結(jié)果如表2所示。圖3給出了在每個(gè)子區(qū)域LHS樣本量大小為15的實(shí)驗(yàn)設(shè)計(jì)點(diǎn)集與正交設(shè)計(jì)點(diǎn)集的分布，圖3還給出了作為對(duì)比的樣本量大小為80的LHS設(shè)計(jì)的實(shí)驗(yàn)點(diǎn)集的分布。

表2 不同樣本量條件下方法的建模性能對(duì)比

(續(xù)表)

3.3 結(jié)果討論

(1)關(guān)于模型預(yù)測(cè)性能。從表2可以看出，采用本文所提出的兩階段的實(shí)驗(yàn)設(shè)計(jì)，在樣本量相同的情況下，所擬合的Kriging模型的預(yù)測(cè)性能均優(yōu)于傳統(tǒng)的LHS。具體來說，在三次對(duì)比中，所提方法的Kriging模型的MSE的均值比傳統(tǒng)LHS的MSE的均值降低了52.4%;MaxE的均值降低了19.2%;StdE的均值則降低了28.8%。

(2)從圖3可以看出，本文所提方法實(shí)驗(yàn)點(diǎn)的安排與傳統(tǒng)LHS有顯著區(qū)別。在輸出特性變化較大的區(qū)域，實(shí)驗(yàn)點(diǎn)密度較高，且分布較為均勻分散，而在其他區(qū)域，實(shí)驗(yàn)點(diǎn)的密度較低，但分布較為均勻整齊。這種實(shí)驗(yàn)點(diǎn)的分布結(jié)合了正交設(shè)計(jì)與LHS的優(yōu)點(diǎn)，使得實(shí)驗(yàn)點(diǎn)在整個(gè)可行域內(nèi)的分布具有較強(qiáng)的均勻性和針對(duì)性。與之對(duì)比，傳統(tǒng)LHS設(shè)計(jì)的實(shí)驗(yàn)點(diǎn)雖然在可行域全局范圍內(nèi)的較為均勻分散，但是并沒有針對(duì)性，導(dǎo)致有相當(dāng)數(shù)量的實(shí)驗(yàn)點(diǎn)安排在輸出特性變化不顯著的區(qū)域，這在很大程度上影響了模型性能。

4 結(jié)論

本文研究了計(jì)算機(jī)實(shí)驗(yàn)過程的小樣本實(shí)驗(yàn)設(shè)計(jì)及過程建模方法。理論研究與仿真分析表明，對(duì)于輸出特性與輸入因子之間作用關(guān)系復(fù)雜的計(jì)算機(jī)仿真過程，通過正交設(shè)計(jì)可以初步探知輸出特性波動(dòng)較大的子區(qū)域，而后再利用LHS可以有針對(duì)性地安排實(shí)驗(yàn)點(diǎn)。在樣本量相同的條件下，本文所提方法的實(shí)驗(yàn)設(shè)計(jì)效率與模型的預(yù)測(cè)性能均優(yōu)于傳統(tǒng)的LHS設(shè)計(jì)。

［1］馮吉路，姜增輝.基于abaqus高速切削ti-6al-4v切削狀態(tài)的有限元仿真［J］.組合機(jī)床與自動(dòng)化加工技術(shù)，2013(2)：47－53.

［2］Pronzato，L.，M?ller，W.，Design of Computer Experiments：Space Filling and Beyond［J］.Statistics and Computing，2012，22(3)：681 －701.

［3］Zhou，X.，Ma，Y.，Tu，Y.，et al.，Ensemble of Surrogates for Dual Response Surface Modeling in Robust Parameter Design［J］.Quality and Reliability Engineering International，2013，29(2)：173 －197.

［4］陳巍，周雄輝，張汝珍，等.基于kriging代理模型的注塑產(chǎn)品翹曲優(yōu)化［J］.上海交通大學(xué)學(xué)報(bào)，2010，44(4)：588－592.

［5］Deng，H.，Shao，W.，Ma，Y.，et al.，Bayesian Metamodeling for Computer Experiments Using the Gaussian Kriging Models［J］.Quality and Reliability Engineering International，2012，28(4)：455 －466.

［6］Can，B.，Heavey，C.，A Comparison of Genetic Programming and Artificial Neural Networks in Metamodeling of Discrete-Event Simulation Models［J］.Computers ＆amp;Operations Research，2012，39(2)：424 －436.

［7］崔慶安.面向多極值質(zhì)量特性的全局式序貫性實(shí)驗(yàn)設(shè)計(jì)方法［J］.系統(tǒng)工程理論與實(shí)踐，2012，32(10)：2143－2153.

［8］Gramacy，R.B.，Polson，N.G.，Particle Learning of Gaussian Process Models for Sequential Design and Optimization［J］.Journal of Computational and Graphical Statistics，2011，20(1)：102 －118.

［9］Dam，E.R.v.，Hertog，D.d.，Husslage，B.G.M.，et al.，Space-Filling Latin Hypercube Designs for Computer Experiments［J］.Optimization and Engineering，2011，12(4)：611－630.

［10］Fang，K.-T.，Lin，D.K.J.，Uniform Design in Computer and Physical Experiments［M］.in The Grammar of Technology Development，2008：105－125.