【MATLAB第110期】#保姆级教学 | 基于MATLAB的PAWN全局敏感性分析方法（无目标函数）含特征变量置信区间分析

news2024/12/27 17:45:56

【MATLAB第110期】#保姆级教学 | 基于MATLAB的PAWN全局敏感性分析方法（无目标函数）含特征变量置信区间分析

一、介绍

PAWN（Probabilistic Analysis With Numerical Uncertainties）是一种基于密度的全局敏感性分析（Global Sensitivity Analysis, GSA）方法，它通过累积分布函数来评估模型参数对输出的影响，考虑输出的分布函数在参数取值不同时的差异，利用Kolmogorov-Smirnov（KS）统计量来表征参数的灵敏度，特别适合于处理输出具有复杂分布特性（非对称或多峰）的情况。

Sobol方法是基于方差的全局敏感性分析方法，它通过计算模型输出方差的贡献来评估参数的重要性。Sobol方法的核心在于方差分解，即认为输出的总方差可以分解为各个参数的方差贡献和参数间交互作用的方差贡献。

PAWN方法是一种基于矩独立的全局敏感性分析方法，它不依赖于输出分布的特定矩（如方差），而是使用整个输出分布来表征不确定性。PAWN通过比较无条件输出分布（所有参数自由变化时）和条件输出分布（固定一个参数时）的差异来度量敏感性。

与Sobol方法相比，优势体现在：

1、适用范围更广：PAWN方法适用于更广泛的模型，特别是那些输出分布非对称或具有多峰性的模型，而Sobol’方法适用于那些输出分布接近正态分布的模型，若输出的分布为非对称时，会产生相应偏差。

2、需要样本量更少：PAWN方法需要的样本数量相对较少，能够更有效地区分参数的相对重要性，而Sobol方法可能需要更多的样本来获得准确的参数排序和灵敏度指数。

3、参数灵敏度区别明显：PAWN方法能够更明显地区分参数的灵敏度，提供了更细致的参数重要性排序，而Sobol方法计算得到的参数灵敏度相对更小，且不同参数间的灵敏度差异也更小，可能不能很明显地对参数的灵敏度进行划分。

二、运行步骤

在这里插入图片描述
输出满足多峰分布。

1、建立代理模型
（1）数据设置：常用的案例数据，103*8 ，前7列代表输入变量M，最后1列代表因变量。
（2）数据变量上下限
X1-X7的最小值为：[137 0 0 160 4.4 708 650]
X1-X7的最大值为：[374 193 260 240 19 1049.90 902]
（3）代理模型选择：BP

2、建立模型

（1）有/无条件输出分布样本数设置：

无条件输出分布：所有参数自由变化时 （传统方法）
有条件输出分布：固定其中一个参数时，其余参数自由变化

NU=150;%无条件CDF的样本数，即当所有输入都变化时
NC=150;%有条件CDF的样本数， 当一个或多个输入固定时
n=10；%分布水平，有条件样本，固定一个输入变量时，选取的值为上下限空间采样中的的10个值（通过lhs采样）

（2）无条件输出分布数据采样
数据抽样采用拉丁超立方抽样lhs方法，但因为一般lhs抽样分布在【0，1】之间，为了解决此问题，采用连续均匀逆累积分布函数解决。

p = lhcube(Nu,M)% 拉丁超立方抽样
Xu=unifinv（p，a，b）
Yu = model_evaluation(bpfun,Xu) ;

p为lhs抽样后的矩阵，区间范围为【0,1】。
a和b分别为变量取值的最小值和最大值。
Xu为以p中相应的概率计算参数空间范围内的均匀cdf的逆。
Yu为以BP模型为代理模型情况下，输入为Xu计算的值

均匀cdf的逆为：
在这里插入图片描述
标准均匀分布时，a=0，b=1。

（3）有条件输出分布数据采样
每个变量在范围空间里抽取n=10个值，分别固定对应变量的值，其余值自由抽样。
举例：针对X1[137,374]进行抽样10个值，得到101507矩阵，其余6个变量同理。

240.769544971060
310.436967493479
145.692630610127
165.239449093029
283.868192881610
208.624632109356
345.374320503857
259.777657520337
362.342087092143
191.097454825086

11507矩阵
在这里插入图片描述
21507矩阵

101507矩阵

XX=unifinv（p，a，b） %同理，只不过矩阵结构是7*10*150*7
YY = model_evaluation(bpfun,XX) ;%同理，只不过矩阵结构是7*10*150*1

（4）估算无条件和有条件的CDF值

[ YF, Fu, Fc  ] = cdfs(Yu,YY) ;
%YF=通过分析Yu和YY中的输出最大最小值，重新对Y值进行等距分布。
%FU=根据Yu估算的经验无条件CDF值
%FC=根据YY{i，k}估算的经验有条件CDF值

取YY （矩阵结构7×10×150×1 ）值和Yu值（矩阵结构150×1）样本中的最大Ymax和最小Ymin值进行分析，考虑整体样本输出分布。

取间隔值为：

deltaY = (Ymax-Ymin)/4000 ;   %  间隔值

则YF公式为：

YF =  [ (Ymin-std(Yu)/10) : deltaY : (Ymax+std(Yu)/10) ]' ;

FU的计算与Yu和YF值相关，FU的矩阵结构与YF相同，

当Nc=150时，F = [1:Nc]/Nc ，即将1等比例分为150份，累计增加。
在N=1:Nc情况下，Yu中小于YF（N）的部分对应序列的值 ，赋对应的F值。

FC则是与YY和YF值相关，矩阵结构为710（YF*1）
与FU计算同理。

（5）CDF值绘图：
灰色的图：YF为输入，FC为输出
红色的图：YF为输入，FU为输出
在这里插入图片描述

（6）KS统计量计算

计算经验无条件CDF和不同条件值下的条件CDF之间的Kolmogorov-Smirnov（KS）统计量：

KS(k,i) = max | F(y) - F(y|x(i)=x_ik) |

每个变量的KS统计值：
在这里插入图片描述
（7）PAWN灵敏性值计算
根据KS的max值进行计算

（8）变量灵敏性值区间范围计算
通过bootstrapping自举法评估鲁棒性，使用越来越多的样本重复计算，以评估收敛性，重复收敛分析以得出置信区间
重采样次数选为100.
在这里插入图片描述
（9）收敛性分析
即对不同样本数量进行分析，选取的采样空间进行扩充，
如NCb = [ NC/10 NC/2 NC ] ;
NUb = [ NU/10 NU/2 NU ] ;NN = NUb+n*NCb ;
即有条件和无条件输出分布样本数量为【165 ，825 ，1650】
纵坐标为不同变量对应的灵敏性值
在这里插入图片描述
（10）考虑输出阈值情况的敏感性分析
可参考【MATLAB第109期】基于MATLAB的带置信区间的RSA区域敏感性分析方法，无目标函数
将PAWN应用于限定输出范围的区域，计算对应的PAWN灵敏性值，
即对应的输出要高于给定阈值，如阈值为10，则10以下的输出要筛除。