基于偏導(dǎo)數(shù)的Sobol’總測(cè)度指標(biāo)的上下限分析

宋述芳　周　桐　呂震宙

1.西北工業(yè)大學(xué),西安,710072　　2.香港理工大學(xué),香港,999077

?

基于偏導(dǎo)數(shù)的Sobol’總測(cè)度指標(biāo)的上下限分析

宋述芳1周桐2呂震宙1

1.西北工業(yè)大學(xué),西安,7100722.香港理工大學(xué),香港,999077

摘要：以Sobol’主測(cè)度指標(biāo)Si作為總測(cè)度指標(biāo)i的下限，建立并推導(dǎo)了基于偏導(dǎo)數(shù)的測(cè)度指標(biāo)作為i的新上限?；诜汉虴uler-Lagrange等式，進(jìn)行了不同變量分布形式下(均勻、正態(tài)、指數(shù)、Beta、三角分布等)，Sobol’總測(cè)度指標(biāo)i的基于偏導(dǎo)數(shù)的上限分析，并給出了新上限詳細(xì)的推導(dǎo)過(guò)程和具體的計(jì)算公式。通過(guò)簡(jiǎn)單數(shù)值和工程算例，驗(yàn)證了新上限的精度及效率，為更準(zhǔn)確地界定總測(cè)度指標(biāo)i的取值區(qū)間提供了參考。

關(guān)鍵詞：總測(cè)度指標(biāo)；主測(cè)度指標(biāo)；基于偏導(dǎo)數(shù)的測(cè)度指標(biāo)；Euler-Lagrange等式

0引言

全局靈敏度(又稱(chēng)重要測(cè)度)指標(biāo)可以全面反映模型輸入變量的不確定性對(duì)模型輸出響應(yīng)不確定性的貢獻(xiàn)程度，它在工程設(shè)計(jì)及概率安全評(píng)估中具有很重要的作用[1-6]。依據(jù)全局靈敏度指標(biāo)的大小，為基本變量進(jìn)行重要性排序，進(jìn)而在設(shè)計(jì)和優(yōu)化中優(yōu)先或者重點(diǎn)考慮重要性程度高的基本變量或者忽略重要性程度低的基本變量，對(duì)系統(tǒng)工程設(shè)計(jì)和優(yōu)化具有重要的指導(dǎo)作用。

1Sobol’的測(cè)度指標(biāo)

文獻(xiàn)[1-2,6]利用方差分析(analysis of variance, ANOVA)給出了基于方差的重要性測(cè)度指標(biāo)。因而，一般的ANOVA分解被認(rèn)為是基于方差的重要性分析的基礎(chǔ)。ANOVA分解的前提是假設(shè)各子項(xiàng)期望為零且各子項(xiàng)正交，在這種假設(shè)下ANOVA的分解唯一，即功能函數(shù)Y=f(x)的唯一分解式為

(1)

(2)

fij(xi,xj)表征兩變量(xi,xj)的相互作用，它可以通過(guò)下式求得:

fi(xi)-fj(xj)-f0

(3)

類(lèi)似地，可得到s(s=3,4,…,n)個(gè)變量的分解項(xiàng)fi1i2…is(xi1,xi2,…,xis)，也可稱(chēng)之為s階交叉項(xiàng)，可通過(guò)對(duì)除這s個(gè)變量外的所有變量求數(shù)學(xué)期望后減去這s階變量的任意子集影響及常數(shù)項(xiàng)f0得到。

考慮式(1)各分解項(xiàng)的正交特點(diǎn)，功能函數(shù)Y=f(x)的方差:

可以由各分解項(xiàng)方差之和表示，即

(4)

(5)

定義與變量xi有關(guān)的子項(xiàng)之和為ui(x)[14]，即

f1…n(x1,…,xn)=f(x)-∫Rf(x)ρXi(xi)dxi

(6)

(7)

(8)

Di=var[E(Y|xi)]

(9)

(10)

其中，E(·)為期望算子，var(·)為方差算子，下標(biāo)“-i”表示除xi以外的變量，即x-i=(x1,…,xi-1,xi+1,…,xn)，從而可得到輸入變量xi的重要性測(cè)度兩個(gè)指標(biāo)的計(jì)算表達(dá)式分別為

(11)

(12)

求解Sobol’測(cè)度指標(biāo)的方法有很多，如基于MC或QMC模擬的方法(如式(11)、式(12))，隨機(jī)抽樣的高維模型替代法(HDMR)等(如式(5))，這些方法多涉及雙重抽樣分析，計(jì)算量較大。我們更期望找到較簡(jiǎn)單的計(jì)算方式以獲得總測(cè)度指標(biāo)的區(qū)間，進(jìn)而確定變量的重要性排序。

2文獻(xiàn)中的Sobol’總測(cè)度指標(biāo)的上下限

2.1下限(lowerbound,LB)

2.2上限(upper bound, UB)

Sobol’等[12]提出的基于偏導(dǎo)數(shù)的測(cè)度指標(biāo)(DGSM)定義為對(duì)模型輸出偏導(dǎo)數(shù)平方的積分。假設(shè)模型函數(shù)Y=f(x)，其輸入變量x=(x1,x2,…,xn)是相互獨(dú)立的隨機(jī)向量，其聯(lián)合概率密度函數(shù)和累積分布函數(shù)分別為ρX(x)和FX(x)，如果?f/?xi存在并且平方可積，變量xi的DGSM指標(biāo)υi為

(13)

表1　滿(mǎn)足log-concave概率分布情況下的

定理1在以下假設(shè)條件下：①隨機(jī)輸入變量相互獨(dú)立；②函數(shù)f(x)為實(shí)函數(shù)；③f(x)的一階導(dǎo)數(shù)為實(shí)函數(shù)；④隨機(jī)變量xi的分布是Boltzman概率測(cè)度，有

(14)

其中，F(xiàn)Xi(xi)為隨機(jī)變量xi的累積分布函數(shù)。

定理2在以下假設(shè)條件下：①隨機(jī)輸入變量相互獨(dú)立；②函數(shù)f(x)為實(shí)函數(shù)；③f(x)的一階導(dǎo)數(shù)為實(shí)函數(shù)；④隨機(jī)變量xi的分布是log-concave概率測(cè)度, 有

(15)

3.1均勻分布

(16)

并且當(dāng)且僅當(dāng)u是常數(shù)的時(shí)候等號(hào)成立。

(17)

其中，Hn為n維的超立方體空間，ui為與變量xi有關(guān)的項(xiàng)。

根據(jù)變量代換可得出(a,b)區(qū)間上的均勻分布變量的總測(cè)度指標(biāo)的新上限為

(18)

3.2其他變量分布類(lèi)型

假設(shè)對(duì)于其他分布，構(gòu)造如下形式的上限：

(19)

假設(shè)Φ=Φ[u]是一個(gè)關(guān)于u(t)的泛函，即

(20)

將積分不等式轉(zhuǎn)化為變分法求泛函極值的問(wèn)題，當(dāng)u(t)滿(mǎn)足∫Ru(t)ρ(t)dt=0時(shí)使得泛函Φ[u]最小。假設(shè)泛函有極值函數(shù)u*=t-μ，其中μ是隨機(jī)變量t的平均值，此極值函數(shù)需滿(mǎn)足Euler-Lagrange方程：

(21)

表2　滿(mǎn)足log-concave概率分布情況下的新上限A(xi)函數(shù)

表2中的常數(shù)K1、K2、K3和K4分別為

K2≥0

4算例及分析

將文獻(xiàn)[12-13]提出的基于偏導(dǎo)數(shù)的測(cè)度指標(biāo)DGSM稱(chēng)為γUB1，本文所提出的新上限稱(chēng)為γUB2，通過(guò)對(duì)不同分布、不同類(lèi)型的功能函數(shù)進(jìn)行測(cè)試，驗(yàn)證其正確性與有效性。

4.1數(shù)值算例

4.1.1均勻分布算例

算例1～3的函數(shù)及測(cè)度指標(biāo)結(jié)果顯示列于表3，其中所有輸入變量均服從(0,1)區(qū)間上的均勻分布，變量之間相互獨(dú)立。

表3　算例1～3的函數(shù)及測(cè)度指標(biāo)結(jié)果對(duì)比

圖1　算例3的隨α的變化趨勢(shì)圖

4.1.2指數(shù)分布算例

算例4～6的函數(shù)及測(cè)度指標(biāo)結(jié)果參見(jiàn)表4，其中輸入變量x1,x2,x3,x4分別服從均值為1,2,3,4的指數(shù)分布，變量之間相互獨(dú)立。

表4　算例4～6的函數(shù)及測(cè)度指標(biāo)結(jié)果對(duì)比

圖2　算例5的隨α的變化趨勢(shì)圖

工程算例1懸臂梁結(jié)構(gòu)在自由端受到載荷P作用，以自由端位移不超過(guò)0.004 m建立極限狀態(tài)方程為

其中，彈性模量E=200 GPa。將梁的長(zhǎng)L、寬b和高h(yuǎn)看作隨機(jī)變量，其均值分別為0.5 m，0.02 m，0.05 m，標(biāo)準(zhǔn)差為0.05 m，0.002 m，0.005 m，其中，L和h為對(duì)數(shù)正態(tài)分布變量，b為正態(tài)分布變量，載荷P為(200，400)N的均勻分布變量。該算例的測(cè)度指標(biāo)對(duì)比表參見(jiàn)表5。

表5　工程算例1的測(cè)度指標(biāo)結(jié)果對(duì)比

從表5可以看出，兩種上限給出的變量重要性排序是一致的，且都能夠很好地找出對(duì)輸出響應(yīng)影響小的兩個(gè)參數(shù)b和h，均勻分布時(shí)指標(biāo)γUB2比γUB1好，正態(tài)分布時(shí)指標(biāo)γUB1≤γUB2，當(dāng)且僅當(dāng)為正態(tài)線性函數(shù)時(shí)，等號(hào)成立。

工程算例2基于材料的蠕變和疲勞試驗(yàn)數(shù)據(jù)，并考慮一級(jí)載荷水平，文獻(xiàn)[17]采用如下非線性極限狀態(tài)方程來(lái)定義失效與安全的邊界線：

g(Nc,Nf,nc,nf,θ1,θ2)=

Dc=nc/NcDf=nf/Nf

其中，θ1和θ2為從試驗(yàn)數(shù)據(jù)中得到的兩個(gè)參數(shù)；Nc與Nf分別為蠕變壽命和疲勞壽命；nc和nf分別為蠕變載荷和疲勞載荷作用的實(shí)際周次。

文獻(xiàn)[17]依據(jù)試驗(yàn)分析，假定上述極限狀態(tài)方程中的基本隨機(jī)變量Nc、Nf、nc、nf均服從對(duì)數(shù)正態(tài)分布，θ1和θ2服從正態(tài)分布(表6)，則該算例的測(cè)度指標(biāo)對(duì)比參見(jiàn)表7。

表6　基本變量分布形式

表7　工程算例2的測(cè)度指標(biāo)結(jié)果對(duì)比

工程算例3單層單跨結(jié)構(gòu)(圖3)的極限狀態(tài)方程為

f(x)=0.01 m-u3(A1,A2,P)=

圖3　單層單跨結(jié)構(gòu)

SiStotiγUB1γUB2A10.996370.997621.023521.19883A20.002290.002330.002340.00239P0.000360.000380.002410.00038

5結(jié)語(yǔ)

本文圍繞Sobol’總測(cè)度指標(biāo)和基于偏導(dǎo)數(shù)的測(cè)度指標(biāo)進(jìn)行了如下工作：建立了(0,1)區(qū)間上均勻分布變量的偏導(dǎo)數(shù)測(cè)度指標(biāo)，并將其作為Sobol’總測(cè)度指標(biāo)的新上限，在此基礎(chǔ)上基于泛函和Euler-Lagrange等式，推導(dǎo)出了正態(tài)分布、指數(shù)分布、Beta分布、三角分布等分布形式的偏導(dǎo)數(shù)測(cè)度指標(biāo)公式。采用QMC方法，在輸入變量服從不同分布，不同的功能函數(shù)情況下，計(jì)算出不同測(cè)度指標(biāo)，并比較它們之間的關(guān)系，從理論和數(shù)值上證實(shí)了所提指標(biāo)的可靠性。通過(guò)算例比較分析得出，對(duì)于均勻分布來(lái)說(shuō)，大多數(shù)情況下，γUB2優(yōu)于γUB1，對(duì)于指數(shù)分布來(lái)說(shuō)，當(dāng)函數(shù)非線性程度不高時(shí)，γUB2的優(yōu)越性還是很明顯的，但隨著函數(shù)非線性程度的提高，γUB1的優(yōu)勢(shì)逐漸顯現(xiàn)出來(lái)。此外，γUB2和γUB1確定的變量重要性排序和總測(cè)度指標(biāo)基本一致。

參考文獻(xiàn)：

[1]呂震宙,李璐祎,宋述芳,等. 不確定性結(jié)構(gòu)系統(tǒng)的重要性分析理論與求解方法[M].北京:科學(xué)出版社，2015.

[2]Sobol’IM.GlobalSensitivityIndicesforNonlinearMathematicalModelsandTheirMonteCarloEstimates[J].MathematicsandComputerinSimulation, 2001, 55(1):221-280.

[3]BorgonovoE.MeasuringUncertaintyImportance:InvestigationandComparisonofAlternativeApproaches[J].RiskAnalysis, 2006, 26(5): 1349-1361.

[4]BorgonovoE.ANewUncertaintyImportanceMeasure[J].ReliabilityEngineering&SystemSafety, 2007, 92(6):771-784.

[5]SaltelliA.SensitivityAnalysisforImportanceAssessment[J].RiskAnalysis, 2002, 22(3): 579-590.

[6]SaltelliA,AndresM,CampolongoF,etal.GlobalSensitivityAnalysis:thePrimer[M].NewYork:JohnWiley&SonsLtd., 2008.

[7]FeilB,KucherenkoS,ShahN.ComparisonofMonteCarloandQuasiMonteCarloSamplingMethodsinHighDimensionalModelRepresentation[C]//TheFirstInternationalConferenceonAdvancesinSystemSimulation.Porto,SIMUL, 2009:12-17.

[8]Sobol’IM,KucherenkoS.OnGlobalSensitivityAnalysisofQuasi-MonteCarloAlgorithms[J].MonteCarloMethodsandSimulation, 2005, 11(1): 1-9.

[9]LiGY,HuJS,WangSW,etal.RandomSampling-highDimensionalModelRepresentation(RS-HDMR)andOrthogonallyofItsDifferentOrderComponentFunctions[J].JournalofPhysicalChemistryA, 2006, 110(7): 2474-2485.

[10]FeilB,KucherenkoS,ShahN.VolatilityCalibrationUsingSplineandHighDimensionalModelRepresentationModels[J].Wilmott, 2009, 1(2): 179-195.

[11]KucherenkoS,SongS,ComparisonofDifferentNumericalEstimatorsforMainEffectGlobalSensitivityIndices[C]//1stECCOMASThematicConferenceonUncertaintyQuantificationinComputationalSciencesandEngineering.CreteIsland, 2015:22-53.

[12]Sobol’IM,KucherenkoS.DerivativebasedGlobalSensitivityMeasuresandTheirLinkwithGlobalSensitivityIndices[J].MathematicsandComputersinSimulation, 2009, 79(10): 3009-3017.[13]IoossB,PopelinAL,BlatmanG,etal.SomeNewInsightsinDerivative-basedGlobalSensitivityMeasures[C]//ProceedingsofPSAM11 &ESREL2012Conference.Helsinki, 2012: 1094-1104.

[14]KucherenkoS,SongS.Derivative-basedGlobalSensitivityMeasuresandTheirLinkwithSobol’SensitivityIndices[C]//MonteCarloandQuasi-MonteCarloMethods.KULeuven:SpringerInternationalPublishingSwitzerland, 2014: 23-38.

[15]BobkovSG.IsoperimetricandAnalyticInequalitiesforLog-concaveProbabilityMeasures[J].TheAnnalsofProbability, 1999, 27(4): 1903-1921.

[16]HardyGH,LittlewoodJE,PolyaG.Inequalities[M]. 2nded.Cambridge:CambridgeUniversityPress, 1973.

[17]MaoHY,MahadevanS.ReliabilityAnalysisofCreep-fatigueFailure[J].InternationalJournalofFatigue, 2000, 22(9): 789-797.

(編輯王艷麗)

收稿日期：2015-08-03

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(NSFC51308459)；中央高?；究蒲袠I(yè)務(wù)費(fèi)專(zhuān)項(xiàng)資金資助項(xiàng)目(310201401JCQ01014，3102015BJ(Ⅱ)CG009)

中圖分類(lèi)號(hào)：TP302.7

DOI：10.3969/j.issn.1004-132X.2016.13.014

作者簡(jiǎn)介：宋述芳,女,1982年生。西北工業(yè)大學(xué)航空學(xué)院副教授。主要研究方向?yàn)轱w行器設(shè)計(jì)、飛行器可靠性工程。發(fā)表論文30余篇。周桐，男,1993年生。香港理工大學(xué)機(jī)械工程學(xué)院博士研究生。呂震宙,女,1966年生。西北工業(yè)大學(xué)航空學(xué)院教授、博士研究生導(dǎo)師。

Analyses for Lower and Upper Bounds of Sobol’ Total Sensitivity Index Based on Derivative

Song Shufang1Zhou Tong2Lü Zhenzhou1

1.Northwestern Polytechnical University，Xi’an,710072 2.Hongkong Polytechnical University,Hongkong,999077

Abstract:A main sensitivity index Si was set as the lower bound of i and the new upper bound of i was built based on the derivative. On the basis of functional analysis and Euler-Lagrange equation, the new upper bound of i based derivative was analyzed for different variable distribution types, such as uniform, normal, exponential, triangular, Beta distribution etc. The derived process and formulas were presented in detail. Several numerical and engineering examples were used to verify the precision and efficiency of the presented bounds, which may provide the accurate bounds of i.

Key words：total global sensitivity index;main global sensitivity index;derivative based global sensitivity index;Euler-Lagrange equation