基于偏導數的全局靈敏度指標的高效求解方法

馮凱旋，呂震宙, ，蔣獻

1. 西北工業大學 航空學院，西安 710072 2. 中國飛行試驗研究院 飛機所，西安 710089

隨著計算機科學技術和數值計算方法的快速發展，在生物技術、環境科學、飛機設計、車輛工程等科學或工程領域，研究者們提出了許多復雜的計算模型[1-2]。這些計算模型通常包括成千上萬個輸入變量，而這些變量的不確定性會對系統的輸出性能產生不同程度的影響[3]。在大多數情況下，只有很少一部分輸入變量對輸出產生重要的影響。因此，如何確定輸入變量對輸出性能影響程度的大小，篩選出重要變量，從而簡化或優化模型是十分必要的。

局部靈敏度定義為模型響應函數在名義值處對輸入變量的偏導數，其值僅反映了響應函數在名義點處的靈敏度信息。全局靈敏度分析方法旨在研究輸入變量在整個取值域內的不確定性對輸出性能的綜合影響，其分析結果可廣泛應用于模型簡化、優化設計和模型確認等方面，更適用于解決上述篩選重要變量的問題。目前，全局靈敏度分析方法可大致分為以下幾類：微分法[2,4]、掃描法[5-6]、方差靈敏度分析[7-9]、矩獨立靈敏度分析[10-13]以及隨機森林[14]等。在上述分析方法中，基于方差的全局靈敏度分析理論應用最為廣泛。在該理論中，輸入變量對模型響應方差的影響被分解為一階影響和高階影響，并在此基礎上提出了兩類靈敏度指標(Sobol 指標)：方差主指標和方差總指標。其中，方差主指標反映了單個輸入變量獨自作用對輸出方差的貢獻，而方差總指標則反映了輸入變量對輸出方差的總貢獻，其中既包含其獨自貢獻，還包括由于與其他變量的交互作用對輸出方差的貢獻?；谄珜档娜朱`敏度分析方法隸屬于微分法，該方法建立在求解響應函數對輸入變量在多個點處的偏導數的基礎上，可以看作是局部靈敏度在全局范圍內的擴展，相比于局部靈敏度分析理論，具有更為豐富的內涵。與此同時，Lamboni等[2]提出的基于偏導數的全局靈敏度指標與方差總指標之間存在直接的定量不等式關系，因此該指標既可以看作是局部靈敏度指標的一種擴展，又可以視為方差總指標一種很好的近似。

目前，基于偏導數的全局靈敏度指標的計算方法討論較少，應用最多的仍是經典的數字模擬方法-Monte Carlo (MC)方法。該方法雖然具有較高的計算精度和廣泛的模型適用性，但其需要大量的樣本才能得到穩定的結果，其計算量在實際工程問題中往往是無法接受的。因此，本文針對基于偏導數的全局靈敏度指標，提出了一種高效的求解方法。該方法首先利用乘法降維公式將模型響應函數展開為連乘積的形式，從而將靈敏度指標中的高維積分轉化為多個一維積分的連乘積，然后利用高斯積分公式對一維積分進行近似求解。與此同時，在求解基于偏導數的靈敏度指標時，需要求解特定點的偏導數，以往大多采用向前差分法或中心差分法，此類方法的求解精度在很大程度上依賴于差分步長的選取，當步長選取不合適時，有可能得出錯誤的結果。因此，本文采用復數步長方法進行偏導數的求解，提高了靈敏度指標的計算精度。數值算例和工程算例驗證了本文所提方法的準確性和高效性。

1 兩類全局靈敏度分析理論

1.1 基于方差的全局靈敏度分析理論

設模型響應函數為Y=g(X),Y為一維輸出，X=[X1X2…Xn]為n維輸入隨機變量。根據高維模型展開(High Dimensional Model Representation，HDMR)理論，當各輸入變量相互獨立時，g(X)可唯一展開為

g1,2,…,n(X1,X2,…,Xn)

(1)

式中:常數g0為函數g(X)的均值;i(i=1,2,…,n)為輸入變量的次序;gi(Xi)為依賴于Xi的單變量函數;i1、i2(i1=1,2,…,n；i2=1,2,…,n)為輸入變量的次序;gi1,i2(Xi1,Xi2)為依賴于Xi1、Xi2的雙變量函數，以此類推。

對式(1)兩邊同時求方差，可得

(2)

式中:V為g(X)的方差；Vi為gi(Xi)的方差；Vi1,i2為gi1,i2(Xi1,Xi2)的方差，以此類推。

據此，輸入變量Xi的方差主指標Si可定義為

(3)

輸入變量Xi的方差總指標STi定義為

(4)

式中:S～i為除與Xi相關項外所有方差分量之和與總方差之比。式(3)中的方差主指標Si反映了輸入變量Xi獨自作用對輸出方差的貢獻。式(4)中的方差總指標STi度量了Xi對輸出方差總的貢獻，其中除了包含Xi獨自的貢獻外，還包括Xi與其他所有變量的交互作用對輸出方差的貢獻。

1.2 基于偏導數的全局靈敏度分析理論

局部靈敏度Ei(x*)定義為響應函數在名義值處對輸入變量的偏導數，即

(5)

(6)

式中:vi為局部靈敏度的平方在全域內的積分，其表達式為

(7)

其中：j(j=1,2,…,n)表示輸入變量的次序，fXj(xj)為輸入隨機變量Xj的概率密度函數；Ci為Cheeger常數[2,15]，其值與Xi的分布類型和分布參數有關，具體表達式為

(8)

式中：υi(xi)為與fXi(xi)相關的測度函數;mi為υi(xi)的中位數。表1給出了滿足log-concave概率分布情況下不同分布類型下的Cheeger常數。

由式(6)～式(8)可以看出，靈敏度指標γi可以看做是局部靈敏度在全局范圍內的平均，并綜合考慮了變量的分布形式及分布參數。除此之外，文獻[2,4]還證明了靈敏度指標γi與方差總指標STi之間存在的定量關系：

(9)

由式 (9)可以看出，γi/V是方差總指標STi的一個上限。此外，實際算例表明了γi和STi確定的輸入變量重要性排序基本一致。

表1 υi(xi)函數、中位數mi及Cheeger常數Ci

2 幾種求解導數的數值方法

對于靈敏度指標γi的求解，無論是用文獻[2]中所提出的MC方法還是其他高效的求解算法都不可避免地需要求解類似式(5)中特定點的偏導數。在實際工程問題中，響應函數Y=g(X)往往為隱函數，無法直接解析求得其偏導數函數，此時就需要采用數值方法來近似求出特定點的偏導數值。由于在求解功能函數在特定點對某一維輸入變量的偏導數時，其他維度的輸入變量固定為常數，因此該點偏導數的求解可轉化為：將其他變量在該點的值代入響應函數中，此時響應函數轉化為只含所求變量的單變量函數，然后利用常用的單變量數值求導方法進行求解。因此，為簡化表達，本節以單變量函數Y=g(x)為例來說明幾種數值求導方法的原理及過程。

根據Taylor展開式得到的有限差分法是工程中常用的數值求導方法。該方法常用的兩種形式為向前差分法和中心差分法。其表達式為

向前差分法：

(10)

中心差分法：

(11)

式(10)和式(11)中：h均為差分步長。在實際應用該方法時，經常會面臨“步長困境[16]”的問題，即試圖通過選擇盡可能小的步長來提高估算精度，但當步長過小時，又會因計算機的減消誤差而得到錯誤的結果。實際算例表明，有限差分法由于減消誤差的存在，其計算精度對差分步長比較敏感，合適的差分步長需要經過多次嘗試才能得到，這無疑增加了計算過程的復雜性和計算結果的不穩定性。因此，本節將采用一種新的數值求導的方法-復數步長方法以克服該問題。

將函數y=g(x+ih)在x處進行Taylor展開(i為虛數單位)，即有

(12)

對式(12)左右兩端同時取虛部后整理可得

(13)

式中：Im[·]表示復數的虛部。略去高階項，得到采用復數步長方法的一階導數計算公式為

(14)

與有限差分法相比，復數步長方法在求解函數Y=g(x)的一階導數時只需計算函數在點x+ih處函數值的虛部即可，相較于有限差分法減少了一次函數調用次數。與此同時，該方法在計算過程中并沒有引入由于兩個數值相近的函數值相減所帶來的減消誤差，因此在理論上，采用復數步長方法計算一階導數的精度將隨著步長h的減小而提高。同時，若函數Y=g(x)存在間斷點，采用有限差分法計算其導數時，若差分步長h范圍內包含間斷點，則需要對有限差分法的導數計算公式進行修正。而采用復數步長方法后，即使步長h中包含間斷點，采用該方法只需知道目標函數在間斷點的函數值即可，所求得的導數也無需進行人工修正。因此，在計算基于偏導數的靈敏度指標中涉及偏導數運算時，本文均采用復數步長方法。

3 基于偏導數的全局靈敏度指標的高效近似算法

3.1 基于偏導數的全局靈敏度指標的計算

求解基于偏導數的全局靈敏度指標γi的關鍵在于準確地估算靈敏度指標vi。由式(7)可知，vi為n維積分，直接對其求解較為困難。因此，首先根據乘法降維方法將響應函數Y=g(X)近似表示為[17]

g(X)≈[g(c)]1-n·

(15)

式中:c=[c1c2…cn]為輸入變量的均值向量。式(15)兩邊同時對Xi求偏導，可得

(16)

式(16)左右兩邊同時平方，可得

(17)

對式(17)兩邊同時積分，可得

fXj(xj)dxj·

(18)

步驟2記a為均值向量，即

a=[c1c2…ci…cn]

(19)

令j=1。

(20)

(21)

(22)

(23)

(24)

j=j+1。

步驟4重復步驟3，直至j>5。根據五點高斯積分公式估算以下兩個一維積分：

(25)

(26)

式中:ai為積分系數，當輸入變量Xi服從均勻分布時，ai=0.5；其他情況下，ai=1。

步驟5重復步驟1～4，估算出所有2n個一維積分，并根據式(18)估算出n個靈敏度指標vi(i=1,2,…,n)。

步驟6根據輸入變量的分布形式及分布參數，計算出對應的Cheeger常數Ci，再根據式(6)計算出n個基于偏導數的全局靈敏度指標γi(i=1,2,…,n)。

3.2 所提方法的計算量

對于n個基于偏導數的全局靈敏度指標γi(i=1,2,…,n)的計算，若采用MC方法，其所需計算量為N0=nN，其中N為所需樣本點數目。由于MC方法的根本依據是概率論中的大數定律，需要大量的樣本模擬才能得到穩定收斂的解，因此N值常取104～106，計算量較大。而本文所提方法利用乘法降維公式，將n維積分問題近似轉化為n個一維積分乘積的形式，而對于一維積分問題只需要很少的計算量就能得到較高精度的解。采用本文所提方法進行計算時，所需計算量為N1=nN′+nN′+1=2nN′+1，其中N′為估算每個一維積分所需要的樣本點數目，若采用五點高斯積分，則N′=5，總的計算量為N1=10n+1?？梢钥闯觯cMC方法相比，本文所提方法的計算量顯著降低，且隨輸入變量維數的增加呈線性增長。

4 算例驗證

下面以一個數值算例和兩個工程算例來驗證文中所提方法的準確性和高效性, 并對基于方差的全局靈敏度指標STi和基于偏導數的全局靈敏度指標γi的計算結果加以分析。

4.1 數值算例

B 函數[18]由于其強的非線性和非單調性，在全局靈敏度分析中被廣泛的作為驗證算例。其數學表達式為

(27)

式中:Xi(i=1,2,…,n)為n個相互獨立且服從正態分布的隨機變量；m=n/2。在本例中，選擇n=10,m=5，各輸入變量分布參數如表2所示。分別使用MC方法和本文所提方法計算靈敏度指標vi，結果如圖1所示，并在此基礎上計算出基于偏導數的全局靈敏度指標γi，結果見圖2。其中，使用MC方法所需的計算量為N0=10×104=105，使用本文所提方法所需計算量為N1=10×10+1=101。此外，使用準MC方法計算出各輸入變量的方差總指標STi，結果如圖3所示，所需計算量為N2=106。

表2 B函數輸入隨機變量的分布參數

從圖1、圖2可以看出，本文所提方法的計算結果與MC方法計算的結果非常接近，從而可以得出本文所提方法的準確性。對比兩者的計算量，可見本文所提方法的高效性。對比圖2、圖3可以看出，靈敏度指標γi和STi的排序是一致的，驗證了γi對STi優良的近似特性。

4.2 工程算例

工程算例1懸臂管結構(圖4)[19]具有6個輸入隨機變量，分別為外加載荷F1、F2、P和T，管壁厚度t和外徑d，其分布類型及分布參數見表3。在實際應用中，工程人員常關注的一個輸出量是頂面中心處的等效應力σmax，其表達式為

(28)

式中:

(29)

其中：θ1和θ2分別為F1和F2與豎直方向夾角。

表3 懸臂管結構輸入變量的分布類型及分布參數

分別使用MC方法和本文所提方法計算靈敏度指標vi和γi，結果見圖5和圖6。兩種方法所需的計算量分別為N0=6×104、N1=61。使用準MC方法計算各輸入變量的方差總指標STi，結果見圖7，所需計算量為N2=106。

通過與大量樣本下的MC方法計算結果對比可以看出，本文所提方法在計算基于偏導數的全局靈敏度指標的準確性和高效性。對比圖6和圖7可以看出，兩種靈敏度指標γi和STi計算結果的排序是一致的，且輸入變量d的靈敏度指標遠大于其他輸入變量，即d的不確定性對頂面中心處等效應力σmax的不確定性貢獻最大。因此，設計者若想有效的減小σmax的不確定性，最經濟有效的方法是降低外徑d的不確定性。

工程算例2考慮文獻[20]中某型民機單側襟翼不對稱運動故障樹模型。襟翼傳動機構及其連接關系如下：內襟翼由1、2號作動器驅動，且它們均沒有設置監控內襟翼傾斜角的傳感器，故不能單獨監控內襟翼的傾斜角度；外襟翼由3、4號作動器驅動，且它們均設置了角度傳感器，因而可以單獨監控外襟翼的傾斜角度。襟翼傳動機構最外側的扭力管處裝有位置傳感器，用于監控襟翼所處位置。襟翼位置控制系統是冗余的，由1、2號襟翼控制單元組成，控制系統可以自動隔離故障控制單元的信號，采用正常的控制裝置進行控制。襟翼控制裝置根據各傳感器的監控信號采取相應的控制行為，若監控到系統傾斜或非對稱，則通過動力驅動裝置使襟翼停止運動，從而將傾斜或非對稱控制在安全范圍內。

該系統故障樹模型如圖8所示，頂事件A為該民機“單側襟翼不對稱運動”，包括8個中間事件M1～M8和12個底事件X1～X12，各事件的意義見文獻[20]。模型輸出為頂事件A的發生概率，輸入為12個底事件的發生概率，假設12個輸入變量均服從正態分布，其分布參數如表4所示。在實際工程中，底事件發生概率的不確定性來源于認知不足，因而可通過收集數據而減小。本算例的任務是討論12個輸入變量對模型輸出的全局貢獻以高效減小頂事件發生概率的不確定性。

分別使用MC方法和本文所提方法計算靈敏度指標vi和γi，結果見圖9和圖10。兩種方法所需計算量分別為N0=1.2×105、N1=121。使用準MC方法計算12個輸入變量的方差總指標STi，結果見圖11，所需計算量為N2=106。圖9和圖10的計算結果驗證了本文所提方法在計算基于偏導數的全局靈敏度指標的準確性和高效性。由圖10和圖11可以看出，僅有X1和X2計算所得靈敏度指標值較大，由此可知X1和X2為重要變量，而其余變量均為不重要變量。因此，設計者若想高效的減小頂事件發生概率的不確定性，應該多搜集X1和X2事件的數據以減小其發生概率的不確定性。

表4 故障樹底事件分布參數Table 4 Distribution parameters of basic events

5 結 論

1) 基于偏導數的全局靈敏度指標γi既可以看作是局部靈敏度指標的一種擴展，又可以視為方差總指標STi的一種近似。鑒于全局靈敏度指標γi這種優良的特性，本文提出了一種高效估算γi指標的方法。該方法利用乘法降維公式將靈敏度指標γi中的高維積分問題轉化為多個一維積分連乘積的形式，然后利用高斯積分公式對一維積分進行近似求解。

2) 在求解基于偏導數的靈敏度指標時，需要求解特定點的偏導數，以往方法大多采用有限差分法。此類方法的求解精度對積分步長的大小十分敏感，當積分步長選取不當時，有可能得出錯誤的計算結果。因此，本文采用復數步長方法進行偏導數的求解，其求解精度理論上隨步長的減小而提高。

3) 本文所提方法可適用于復雜的非線性響應函數，通過3個算例進行了驗證。

4) 由于本文所提方法的計算量隨輸入變量維數的增加呈線性增長，因此，可適用于高維問題。

參 考 文 獻

[1] 阮文斌, 劉洋, 熊磊. 基于全局靈敏度分析的側向氣動導數不確定性對側向飛行載荷的影響[J]. 航空學報, 2016, 37(6): 1827-1832.

RUAN W B, LIU Y, XIONG L. Influence of side aerodynamic derivate uncertainty on side flight load based on global sensitivity analysis[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(6): 1827-1832 (in Chinese).

[2] LAMBONI M, IOOSS B, POPELIN A L, et al. Derivative-based global sensitivity measures: General links with Sobol’s indices and numerical tests[J]. Mathematics and Computers in Simulation, 2013, 87: 45-54.

[3] PATELLI E, PRADLWARTER H. Monte Carlo gradient estimation in high dimensions[J]. International Journal for Numerical Methods in Engineering, 2010, 81: 172-188.

[4] SOBOL I M, KUCHERENKO S. Derivative global sensitivity measures and their link with global sensitivity indices[J]. Mathematics and Computers in Simulation, 2009, 79(10): 3009-3017.

[5] MORRIS M D. Factorial sampling plans for preliminary computational experiments[J]. Technimetrics, 1991, 33(2): 161-174.

[6] CAMPOLONGO F, CARIBONI J, SALTELLI A. An effective screening design for sensitivity analysis of large models[J]. Environment Model Software, 2007, 22(10): 1509-1518.

[7] SOBOL I M. Sensitivity analysis for non-linear mathematical, methods[J]. Mathematical Modelling and Computational Experiment, 1993, 1: 407-414.

[8] HOMMA T, SALTELLI A. Importance measures in global sensitivity analysis of nonlinear models[J]. Relia-bility Engineering and System Safety, 1996, 52(1): 1-17.

[9] 鞏祥瑞, 呂震宙, 左鍵巍. 兩種基于方差的全局靈敏度分析W指標改進算法[J]. 航空學報, 2016, 37(6): 1888-1898.

GONG X R, LU Z Z, ZUO J W. Two importance methods for variance-based global sensitivity analysis’ W-indices[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(6): 1888-1898 (in Chinese).

[10] BORGONOVO E. A new uncertainty importance measure[J]. Reliability Engineering and System Safety, 2007, 92(6): 771-784.

[11] LI L Y, LU Z Z, FENG J, et al. Moment-independent importance measure of basic variable and its state dependent parameter solution[J]. Structural Safety, 2012, 38: 40-47.

[12] WEI X B, WEI Y J, HUANG Q X. Experimental investigation of critical ventilating coefficient of ventilated supercavity in water tunnel[J]. Journal of Harbin Institute of Technology, 2007, 39(5): 797-799.

[13] WEI P F, LU Z Z, HAO W R, et al. Efficient sampling methods for global reliability sensitivity analysis[J]. Computer Physics Communications, 2012,183(8): 1728-1743.

[14] BREIMAN L. Random forest[J]. Machine Learning, 2001, 45(1): 5-32.

[15] BOBKOV S G. Isoperimetric and analytic inequalities for Log-concave probability measures[J]. The Annals of Probability, 1999, 27(4): 1903-1921.

[16] MARTINS J, KROO I, ALONSO J. An automated method for sensitivity analysis using complex variables: AIAA-2000-0689[R]. Reston, VA: AIAA, 2000.

[17] YUN W Y, LU Z Z, ZHANG K C, et al. An efficient sampling method for variance-based sensitivity analysis[J]. Structural Safety, 2017, 65: 74-83.

[18] ROCQUIGNY E, DEVICTOR N, TARANTOLA S. Uncertainty in industrial practice[M]. Hoboken: Wiley, 2008.

[19] JIANG C, LI W X, HAN X, et al. Structural reliability analysis based on random distributions with interval parameters[J]. Computers and Structures, 2011, 89(23-24): 2292-2302.

[20] 魏鵬飛. 結構系統可靠性及靈敏度分析研究[D]. 西安: 西北工業大學, 2014: 79-81.

WEI P F. Research on the reliability and sensitivity analysis of structure systems[D]. Xi’an: Northwestern Polytechnical University, 2014: 79-81 (in Chinese).