基于因子函數(shù)的α－IHS抽樣方法

方 磊，許銀龍，汪 進(jìn)，王 芳，胡麗琴，1，＊

（1.中國(guó)科學(xué)技術(shù)大學(xué)核科學(xué)技術(shù)學(xué)院，安徽合肥 230027；2.中國(guó)科學(xué)院核能安全技術(shù)研究所，中國(guó)科學(xué)院中子輸運(yùn)理論與輻射安全重點(diǎn)實(shí)驗(yàn)室，安徽合肥 230031）

基于因子函數(shù)的α－IHS抽樣方法

方 磊1，2，許銀龍1，2，汪 進(jìn)2，王 芳2，胡麗琴2，1，＊

（1.中國(guó)科學(xué)技術(shù)大學(xué)核科學(xué)技術(shù)學(xué)院，安徽合肥 230027；2.中國(guó)科學(xué)院核能安全技術(shù)研究所，中國(guó)科學(xué)院中子輸運(yùn)理論與輻射安全重點(diǎn)實(shí)驗(yàn)室，安徽合肥 230031）

核電廠概率安全分析中，由于數(shù)據(jù)源存在不確定性，導(dǎo)致無(wú)法進(jìn)行準(zhǔn)確評(píng)估，因此需開(kāi)展不確定性分析。樣本在空間分布的均勻特性是不確定性分析的關(guān)鍵因素，不同的樣本分布導(dǎo)致不確定分析結(jié)果差異較大。傳統(tǒng)的拉丁超立方抽樣方法在樣本空間分布均勻性方面未進(jìn)行優(yōu)化，改進(jìn)分布式超立方抽樣（IHS）方法通過(guò)保持樣本點(diǎn)之間的最優(yōu)距離來(lái)實(shí)現(xiàn)空間均勻分布，但其最優(yōu)距離只能在理想分布下達(dá)到最優(yōu)。為改進(jìn)IHS設(shè)計(jì)上的缺陷，提出了基于因子函數(shù)的α－IHS方法，利用修正因子α來(lái)優(yōu)化IHS方法中的最優(yōu)距離。結(jié)果表明，該方法較IHS方法具有更均勻的空間分布，提高了抽樣效率。

IHS方法；最優(yōu)距離；因子函數(shù)；α－IHS方法

Key words：IHS method；optimal distance；factor function；α－IHS method

概率安全評(píng)價(jià)（PSA）方法能對(duì)系統(tǒng)進(jìn)行定量安全評(píng)估，并為系統(tǒng)優(yōu)化提供指導(dǎo)。目前PSA方法已在核能、航天、電子等領(lǐng)域得到了廣泛應(yīng)用［1－5］。但在PSA分析中，由于數(shù)據(jù)源存在不確定性，因而需進(jìn)行不確定分析。不確定分析主要是分析各種系統(tǒng)單元參數(shù)的不確定性對(duì)系統(tǒng)的影響，而在整個(gè)樣本空間上對(duì)參數(shù)進(jìn)行均勻的抽樣是不確定分析的關(guān)鍵。

目前研究者已提出了多種抽樣方法，如常見(jiàn)的蒙卡抽樣、LHS（Latin hypercube sampling）以及優(yōu)化的LHS方法等。蒙卡抽樣是一種普遍應(yīng)用的抽樣方法，其實(shí)現(xiàn)較方便，對(duì)各種復(fù)雜系統(tǒng)均有很強(qiáng)的適用性。但該方法需大量抽樣，時(shí)間消耗大，抽樣效率低，不具有無(wú)坍塌性（non－collapsing）［6］。在蒙卡抽樣的基礎(chǔ)上，Mc Kay等［7］提出了LHS抽樣方法。采用分層抽樣的策略，有效減少蒙卡抽樣方法估計(jì)的被積函數(shù)的誤差，并提高了抽樣效率［8］。但LHS抽樣并未解決樣本的空間填充性（spacefilling）［6］。Beachkofski等［9］提出了改進(jìn)分布式超立方抽樣（IHS）方法，利用最優(yōu)距離的思想，改進(jìn)LHS在空間分布不均勻的缺陷。但I(xiàn)HS方法在選取最優(yōu)距離時(shí)，從理想的超立方體出發(fā)，得到的抽樣結(jié)果并不是最優(yōu)的。

本文提出一種基于修正因子的IHS方法，對(duì)IHS中設(shè)定的最優(yōu)距離增加一修正因子，使IHS具有更好的空間均勻填充性能。不同系統(tǒng)的需求不同，因而參數(shù)個(gè)數(shù)及樣本大小都會(huì)存在一定的差異，1個(gè)修正因子并不能解決所有的情況。針對(duì)大量的抽樣數(shù)據(jù)，分析各種情況下的修正因子，擬合出修正因子函數(shù)，設(shè)計(jì)α－IHS方案，并驗(yàn)證方案的可行性。

1 傳統(tǒng)IHS方法

1.1 方法簡(jiǎn)介

研究者提出了一系列評(píng)價(jià)標(biāo)準(zhǔn)來(lái)反映抽樣方法的空間填充性，本文選擇廣為適用的G（L）作為評(píng)價(jià)標(biāo)準(zhǔn)［10－12］。G（L）的定義為：

其中：N為樣本數(shù)；Xi為樣本；L為抽樣方法，如果G（L1）＜G（L2），則表明L1較L2有更均勻的分布。

圖1為簡(jiǎn)單隨機(jī)抽樣（SRS）、LHS、IHS的比較，其中橫軸和縱軸分別表示樣本參數(shù)的空間，可看到SRS出現(xiàn)坍塌，IHS、LHS沒(méi)有；IHS的G（L）最小，其抽樣效果最好。圖1中，M為空間維數(shù)。

在IHS方法中，把樣本空間看成1個(gè)超立方體，由有限個(gè)相同的小超立方體組合而成，這樣采樣點(diǎn)應(yīng)均勻分布在這些小超立方體中。設(shè)超立方體的體積為NM，要抽取N個(gè)樣本，就要把整個(gè)超立方體分成N個(gè)小超立方體，而每個(gè)小超立方體也是M維的，這樣設(shè)定的最優(yōu)距離就是每個(gè)小超立方體的邊長(zhǎng)。在IHS方法中，對(duì)樣本點(diǎn)間的最小距離作出限制：最小距離須接近于理論上的每個(gè)小超立方體之間的最優(yōu)距離。最優(yōu)距離計(jì)算公式為：

圖1 3種抽樣方法在N＝50、M＝2下的抽樣分布Fig.1 Sampling distribution of three sampling methods at N＝50 and M＝2

具體實(shí)現(xiàn)是，每當(dāng)抽取1個(gè)樣本時(shí)，首先會(huì)抽取多個(gè)備選點(diǎn)，之后從這些備選點(diǎn)中選取最優(yōu)的點(diǎn)作為本次抽樣的結(jié)果，其中選點(diǎn)標(biāo)準(zhǔn)是基于樣本點(diǎn)之間的最優(yōu)距離opt。

常見(jiàn)的輸液報(bào)警器缺乏對(duì)輸液滴速的自動(dòng)控制,比如輸液快輸完了,需要停止輸液,系統(tǒng)自動(dòng)讓滴速控制機(jī)械傳動(dòng)裝置來(lái)讓輸液器停止輸液。滴速調(diào)節(jié)裝置設(shè)計(jì)圖如圖4所示,相應(yīng)的滴速控制機(jī)械傳動(dòng)裝置實(shí)物如圖5所示。輸液滴速控制的核心部分是通過(guò)步進(jìn)電機(jī)的工作來(lái)控制的,圖6是步進(jìn)電機(jī)經(jīng)典的控制電路。

1.2 方法不足

在IHS方法中，選定opt即小超立方體的邊長(zhǎng)作為最優(yōu)距離。但由于實(shí)際采樣中無(wú)坍塌性要求，樣本間的最小距離并不都是opt。

圖2說(shuō)明了具體的情況，假設(shè)需要抽取4個(gè)樣本，每個(gè)樣本是2維的，每一維都服從均勻分布。按照IHS方法：樣本空間首先被分為4等分，樣本點(diǎn)即小超立方體中心，這些樣本點(diǎn)之間的最小距離即為opt＝2。這樣圖2a是一種滿足要求的最佳分布，但樣本有坍塌性。圖2b則是實(shí)際中能實(shí)現(xiàn)的一種較為理想的分布，這時(shí)樣本點(diǎn)間最小距離為2.236，并不是IHS方案中設(shè)定的理論值2。因而可猜測(cè)，有較理論值2更合適的opt使樣本空間分布更加均勻。

圖2 IHS抽樣理想（a）和實(shí)際（b）條件下的分布Fig.2 Ideal（a）and practical（b）distributions of IHS

2 α－IHS方法

2.1 α因子

假定對(duì)于opt，存在一修正因子α，使得：

將opt′作為新的最佳距離，來(lái)優(yōu)化IHS。為驗(yàn)證該優(yōu)化思想，選擇N＝100、M＝10為例，結(jié)果列于表1。從表1可看出，當(dāng)α為1.1～1.5時(shí)，G（L）均有所改進(jìn)，尤其是當(dāng)α＝1.3時(shí)，G（L）出現(xiàn)了最小值，且方差和均值的收斂效果很好，與0～1均勻分布上的均值和方差理論值0.5、1／12逼近。圖3示出了改進(jìn)方法α－IHS的實(shí)現(xiàn)流程。

表1 N＝100、M＝10時(shí)的因子結(jié)果Table 1 Result ofαat N＝100 and M＝10

圖3 α－IHS實(shí)現(xiàn)流程Fig.3 Implementation process ofα－IHS

修正因子的引入改善了IHS抽樣結(jié)果。但表1中的因子α＝1.3僅是一個(gè)例子，并不適用所有情況下的抽樣。進(jìn)而分析N＝100、M＝20的情況，發(fā)現(xiàn)最優(yōu)α＝1.7。為使修正因子的IHS抽樣方法得到普遍的適用，需找出滿足所有抽樣情況的因子函數(shù)指導(dǎo)實(shí)際應(yīng)用。

2.2 α因子函數(shù)

通過(guò)大量的實(shí)驗(yàn)數(shù)據(jù)分析來(lái)擬合因子函數(shù)。IHS原方法中參數(shù)dup值默認(rèn)為5。在實(shí)驗(yàn)中，N從20～200，依次增加20，M從10～100，依次增加10。為避免偶然性，對(duì)于每種抽樣情況的G（L），采取30次實(shí)驗(yàn)求平均值的策略。表2列出了不同抽樣情況下的最優(yōu)因子α，表3、4分別列出了α＝1和α在最優(yōu)情況下的部分G（L）值。通過(guò)表3、4比較表明，在同樣情況下，α－IHS方法的G（L）值得到了很大的改進(jìn)，尤其是在抽樣維數(shù)低的情況下，效果特別明顯。

表2 不同抽樣情況的最優(yōu)因子αTable 2 Optimal factorαfor different sampling situations

表3 α＝1時(shí)IHS方法的G（L）值Table 3 G（L）value of IHS method atα＝1

表4 α最優(yōu)情況下α－IHS的G（L）值Table 4 G（L）value ofα－IHS method at optimal factorα

表2列出了實(shí)驗(yàn)情況下的所有最優(yōu)因子α，但在實(shí)際應(yīng)用中抽取的樣本可能更多。因此，基于表2的結(jié)果，需擬合出因子函數(shù)。使用Matlab擬合出如下的α因子函數(shù)表達(dá)式：

圖4a示出了M＝10、50、100下最優(yōu)因子α隨N的變化，可看到N對(duì)α影響不大；圖4b示出了N＝60、120、200下α隨M的變化，相對(duì)于N，M對(duì)α影響更為顯著。為進(jìn)一步確認(rèn)因子函數(shù)的可用性，將用實(shí)驗(yàn)結(jié)果進(jìn)行驗(yàn)證。

圖4 最優(yōu)因子α隨N（a）和M（b）的變化Fig.4 Optimal factorαvs.N（a）and M（b）

3 實(shí)驗(yàn)結(jié)果

在PSA軟件進(jìn)行故障樹(shù)分析［13－15］時(shí)，通過(guò)基本事件的失效率來(lái)分析頂事件的失效率，而基本事件的失效率則是由諸多參數(shù)來(lái)決定。在對(duì)頂事件失效率進(jìn)行不確定分析時(shí)，需對(duì)所有參數(shù)進(jìn)行抽樣分析。如在中國(guó)鉛基研究堆CLEAR－Ⅰ的1個(gè)子系統(tǒng)中就用到了152個(gè)基本事件，39個(gè)參數(shù)［16－17］。表5列出了實(shí)際應(yīng)用中可能使用到的抽樣情況，根據(jù)因子函數(shù)計(jì)算得到的α。表6列出了具體的G（L），可看出，提出的α－IHS方法在實(shí)際應(yīng)用中是可行的，與IHS方法相比，G（L）均有改善，在低維時(shí)，改進(jìn)效果接近10%。其中，修正比

表5 實(shí)際應(yīng)用中根據(jù)因子函數(shù)計(jì)算得到的最優(yōu)因子αTable 5 Optimal factorαcalculated by factor function in practical application

表6 實(shí)際應(yīng)用中IHS（α＝1）和α－IHS（α最優(yōu)）在G（L）上的比較Table 6 Comparison of G（L）for IHS andα－IHS in practical application

4 結(jié)論

本文通過(guò)對(duì)IHS方法的研究和分析，針對(duì)IHS設(shè)計(jì)中最優(yōu)距離的問(wèn)題進(jìn)行改進(jìn)，提出了基于修正因子的α－IHS方法。本文通過(guò)提出修正因子來(lái)改進(jìn)IHS方法設(shè)計(jì)中的缺陷，優(yōu)化IHS方法中的最優(yōu)距離，并結(jié)合大量實(shí)驗(yàn)數(shù)據(jù)，擬合出了因子函數(shù)，最終通過(guò)實(shí)際系統(tǒng)驗(yàn)證了α－IHS方法相對(duì)于IHS方法的優(yōu)勢(shì)：在確保樣本點(diǎn)均值、方差均具有良好收斂特性的前提下，α－IHS方法使樣本點(diǎn)在整個(gè)樣本空間具有更均勻的分布。

衷心感謝中國(guó)科學(xué)院核能安全技術(shù)研究所FDS團(tuán)隊(duì)其他成員的大力支持和幫助。

［1］QIU L，WU Y，XIAO B，et al.A low aspect ratio tokamak transmutation system［J］.Nuclear Fusion，2000，40：629－633.

［2］WU Y，CHENG M，F(xiàn)DS Team.Development and application of the Chinese female computational phantom rad－HUMAN［C］∥SNA＋MC 2013－Joint International Conference on Supercomputing in Nuclear Applications＋Monte Carlo.Paris：EDP Sciences，2014.

［3］WU Y，JIANG J，WANG M，et al.The fusiondriven hybrid system and its material selection［J］.Journal of Nuclear Materials，2002，307－311：1 629－1 636.

［4］WU Y，F(xiàn)DS Team.CAD－based interface programs for fusion neutron transport simulation［J］.Fusion Engineering and Design，2009，84：7－11.

［5］HUANG Q，LI J，CHEN Y.Study of irradiation effects in China low activation martensitic steel CLAM［J］.Journal of Nuclear Materials，2004，329：268－272.

［6］JANSSEN H.Monte－Carlo based uncertainty analysis sampling efficiency and sampling convergence［J］.Reliability Engineering and System Safety，2013，109：123－132.

［7］Mc KAY M D，BECKMAN R J，CONOVER W J.A comparison of three method for selecting values of input variables in the analysis of output from a computer code［J］.Technometrics，1979，42（1）：55－61.

［8］KUCHERENKO S，ALBRECHT D，SALTELLI A.Comparison of latin hypercube and quasi Monte Carlo sampling techniques［C］∥IBPSA conference.Sydney：［s.n.］，2011.

［9］BEACHKOFSKI B K，GRANDHI R.Improved distributed hypercube sampling［R］.USA：AIAA，2002.

［10］HUSSLAGE B G M，RENNEN G，van DAM E R，et al.Space－filling Latin hypercube designs for computer experiments［J］.Optimization and Engineering，2011，12（4）：611－630.

［11］STOCKI R.A method to improve design reliability using optimal Latin hypercube sampling［J］.Computer Assisted Mechanics and Engineering Sciences，2005，12（4）：393－414.

［12］LIEFVENDAHLA M，STOCKI R.A study on algorithms for optimization of Latin hypercubes［J］.Journal of Statistical Planning and Inference，2006，136：3 231－3 247.

［13］吳宜燦，劉萍，胡麗琴，等.大型集成概率安全分析軟件系統(tǒng)的研究與發(fā)展［J］.核科學(xué)與工程，2007，27（3）：69－75.WU Yican，LIU Ping，HU Liqin，et al.Development of an integrated probabilistic safety assessment program［J］.Chinese Journal of Nuclear Science and Engineering，2007，27（3）：69－75（in Chinese）.

［14］HUANG Q，LI C，WU Q，et al.Progress in development of CLAM steel and fabrication of small TBM in China［J］.Journal of Nuclear Materials，2011，417：85－88.

［15］CAO R，WU Y，PEI X，et al.Multi－objective optimization of inverse planning for accurate radiotherapy［J］.Chinese Physics C，2011，35（3）：313－317.

［16］吳宜燦，胡麗琴，龍鵬程，等.中國(guó)科研信息化藍(lán)皮書(shū)［M］.北京：科學(xué)出版社，2013：232－244.

［17］HU H，WU Y，CHEN M，et al.Benchmarking of SNAM with ITER 3D model［J］.Fusion Engineering and Design，2007，82：2 867－2 871.

α－IHS Sampling Method Based on Factor Function

FANG Lei1，2，XU Yin－long1，2，WANG Jin2，WANG Fang2，HU Li－qin2，1，＊
（1.School of Nuclear Science and Technology，University of Science and Technology of China，Hefei 230027，China；2.Key Laboratory of Neutronics and Radiation Safety，Institute of Nuclear Energy Safety Technology，Chinese Academy of Sciences，Hefei 230031，China）

In the probabilistic safety assessment of nuclear power plant，as for the uncertainty of data source，it’s essential to launch uncertainty analysis.As a key factor of uncertainty analysis，the distribution of sample space would have a great effect on results of uncertainty analysis.The Latin hypercube sampling method doesn’t make any improvement in space－filling of sample space，the improved distributed hypercube sampling method has improved it by keeping the optimal distance of sample points，but the optimal distance could only be achieved in ideal situation.In order to solve such a problem，α－IHS method based on factor function was put forward to optimize the optimal distance.The results prove that theα－IHS method has advantages over IHS method in a better distribution，and improves sampling efficiency.

TL364

：A

：1000－6931（2015）12－2299－06

10.7538／yzk.2015.49.12.2299

2014－09－21；

：2014－11－01

國(guó)際熱核聚變實(shí)驗(yàn)堆（ITER）計(jì)劃專項(xiàng)資助（2104GB112000）；中國(guó)科學(xué)院戰(zhàn)略性先導(dǎo)科技專項(xiàng)資助（XDA03040000）；國(guó)家自然科學(xué)基金資助項(xiàng)目（91026004）；中國(guó)科學(xué)院科技數(shù)據(jù)資源整合與共享工程“重點(diǎn)數(shù)據(jù)庫(kù)”項(xiàng)目資助（XXH12504－1－09）；中國(guó)科學(xué)院合肥物質(zhì)科學(xué)研究院青年火花基金資助項(xiàng)目（YZJJ201327）

方 磊（1990—），女，安徽樅陽(yáng)人，碩士研究生，核能科學(xué)與工程專業(yè)

＊通信作者：胡麗琴，E－mail：liqin.hu＠fds.org.cn