基于偏导数的全局灵敏度指标的高效求解方法

冯凯旋，吕震宙, ，蒋献

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。可以看出，与MC方法相比，本文所提方法的计算量显著降低，且随输入变量维数的增加呈线性增长。

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).