数学建模灵敏性分析（一阶、二阶灵敏度）python代码+懒人专用版-Toy模板网

这篇具有很好参考价值的文章主要介绍了数学建模灵敏性分析（一阶、二阶灵敏度）python代码+懒人专用版。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

其他的不用多说，能找到这篇博客说明读者已经基本了解了一、二阶灵敏度的定义是啥，但苦于其他博客的博大精深，很难用python代码写出自己模型的灵敏度测试。
这里直接给出代码，读者只需要修改模型的函数表达式即可。

from SALib.sample import saltelli
from SALib.analyze import sobol
from pylab import *
import matplotlib.pyplot as plt
import pandas as pd

mpl.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 = 128
param_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])