其他的不用多说，能找到这篇博客说明读者已经基本了解了一、二阶灵敏度的定义是啥，但苦于其他博客的博大精深，很难用python代码写出自己模型的灵敏度测试。
这里直接给出代码，读者只需要修改模型的函数表达式即可。
from SALib.sample import saltellifrom SALib.analyze import sobolfrom pylab import *import matplotlib.pyplot as pltimport pandas as pdmpl.rcParams['font.sans-serif'] = ['SimHei']# 修复中文乱码plt.rcParams['axes.unicode_minus'] = False# 用来正常显示负号np.seterr(divide='ignore', invalid='ignore')# 消除被除数为0的警告# 上面不用管problem = {'num_vars': 4,# 模型变量的数量'names': ['T', 'PT', 'Pr', 'i'],# 模型的变量'bounds': [[2000, 10000], [0, 0.2], [0.1, 0.5], [1, 30], ]# 指定变量的范围，一一对应}def evaluate(xx):# 进行灵敏度分析的模型return np.array([x[0] * np.power(1 - (1 - np.power(1 - x[1], 1.4263 * np.power(x[1], -0.426))) * x[2], x[3] - 1) * np.power(x[1], 0.426) / 1.4263 * 5 + x[0] * np.power(1 - (1 - np.power(1 - x[1], 1.4263 * np.power(x[1], -0.426))) * x[2], x[3] - 1) * ( 1 - np.power(1 - x[1], 1.4263 * np.power(x[1], -0.426))) * 16 for x in xx])'''注意返回的是np.array([function for x in X]) function是函数表达式比如function(T,PT,Pr,i)=T+PT+Pr+i 那么function就写成x[0]+x[1]+x[2]+x[3]很显然，一一对应定义模型中的变量，注意列表下标从0开始'''# 下面不用管samples = 128param_values = saltelli.sample(problem, samples)print('模型运行次数', len(param_values))Y = evaluate(param_values)Si = sobol.analyze(problem, Y)print()print('ST:', Si['ST'])# 总灵敏度print('S1:', Si['S1'])# 一阶灵敏度print("S2 Parameter:", Si['S2'][0, 1])# 二阶灵敏度# 一阶灵敏度与总灵敏度图片df_sensitivity = pd.DataFrame({"Parameter": problem["names"],"一阶灵敏度": Si["S1"],"总灵敏度": Si["ST"]}).set_index("Parameter")fig, axes = plt.subplots(figsize=(10, 6))df_sensitivity.plot(kind="bar", ax=axes, rot=45, fontsize=16)# 二阶灵敏度图片second_order = np.array(Si['S2'])pd.DataFrame(second_order, index=problem["names"], columns=problem["names"])figs, axes = plt.subplots(figsize=(8, 10))ax_image = axes.matshow(second_order, vmin=-1.0, vmax=1.0, cmap="RdYlBu")cbar = figs.colorbar(ax_image)ax_image.axes.set_xticks(range(len(problem["names"])))ax_image.axes.set_xticklabels(problem["names"], rotation=45, fontsize=24)r = ax_image.axes.set_yticklabels([""] + problem["names"], fontsize=24)plt.show()读者只需要修改
problem = {'num_vars': 4,# 模型变量的数量'names': ['T', 'PT', 'Pr', 'i'],# 模型的变量'bounds': [[2000, 10000], [0, 0.2], [0.1, 0.5], [1, 30], ]# 指定变量的范围，一一对应}def evaluate(xx):# 进行灵敏度分析的模型return np.array([x[0] * np.power(1 - (1 - np.power(1 - x[1], 1.4263 * np.power(x[1], -0.426))) * x[2], x[3] - 1) * np.power(x[1], 0.426) / 1.4263 * 5 + x[0] * np.power(1 - (1 - np.power(1 - x[1], 1.4263 * np.power(x[1], -0.426))) * x[2], x[3] - 1) * ( 1 - np.power(1 - x[1], 1.4263 * np.power(x[1], -0.426))) * 16 for x in xx])具体的含义注释已经给出了，需要注意的是，这个代码仅适用于有明确函数表达式的（单个）、多参数的模型，类似于规划问题之类的约束模型不太适用（我太菜了不会写）
不要忘记修改模型变量的数量