analysis - 如何使用 SALib 对 ODE 进行敏感性分析?
问题描述
def iron(XYZ,t,a12,a21,a23,a32,b13,b31,I):
X1,X2,X3=XYZ
dX1=-a12*(X1)+a21*(X2)-b13*(X1)+b31*(X3)
dX2=-a23*(X2)-a21*(X2)+a12*(X1)+a32*(X3)
dX3=-a32*(X3)-b31*(X3)+a23*(X2)+b13*(X1)-I
return dX1,dX2,dX3;
a12=0.0005
a21=0.00001
a23=0.0003
a32=0.0002
b13=0.0001
b31=0.000001
I=0.001
XYZ0=[1000.,30.,10.]
X10=1000.
X20=50.
X30=30.
t=linspace(0,100,1000) #(start,stop,num samples to generate)
XYZ=odeint(iron,XYZ0,t,args=(a12,a21,a23,a32,b13,b31,I))
是否可以使用 SALib 对这个 ODE 系统进行敏感性分析?我想研究模型输入的影响(参数 a 和 b,初始条件)。另外,我可以获得渐近解值吗?
解决方案
以下代码是 Sobol 分析问题的示例实现。每个参数的界限必须根据问题进行调整;我为这个例子假设了一个范围。如果存在通常可以作为参数 (a12, a21,...,b13,b31,I) 的函数找到的渐近解,您可以遵循类似的过程。确定渐近解决方案可能最好作为一个单独的问题发布。
时间结束时的 X1、X2 和 X3 值必须通过 Sobol 方法单独分析。可以为所有时间步的每个 X1、X2 和 X3 计算灵敏度指数,但它需要保存循环每次迭代的所有输出。它还需要多次运行 Sobol 分析。
以下代码示例的一些示例输出是:
====X2 Sobol output====
Parameter S1 S1_conf ST ST_conf
a12 0.409635 0.083264 0.411180 0.049683
a21 0.000002 0.000095 0.000001 0.000000
a23 -0.000626 0.002955 0.000471 0.000057
a32 0.000068 0.000504 0.000017 0.000002
b13 0.000045 0.000232 0.000004 0.000001
b31 0.000000 0.000000 0.000000 0.000000
x1_0 0.430008 0.078269 0.434074 0.053487
x2_0 0.169098 0.051591 0.162944 0.018678
x3_0 -0.000038 0.000335 0.000007 0.000001
示例代码实现
# importing packages
from scipy import integrate as sp
import numpy as np
import SALib
from SALib.sample import saltelli
from SALib.analyze import sobol
# definition of the system of ODEs
def iron(XYZ,t,a12,a21,a23,a32,b13,b31,I):
X1,X2,X3=XYZ
dX1=-a12*(X1)+a21*(X2)-b13*(X1)+b31*(X3)
dX2=-a23*(X2)-a21*(X2)+a12*(X1)+a32*(X3)
dX3=-a32*(X3)-b31*(X3)+a23*(X2)+b13*(X1)-I
return dX1,dX2,dX3;
# default parameter values
a12=0.0005
a21=0.00001
a23=0.0003
a32=0.0002
b13=0.0001
b31=0.000001
I=0.001
# initial condition
XYZ0=[1000.,30.,10.]
X10=1000.
X20=50.
X30=30.
# tmie steps
t=np.linspace(0,100,1000) #(start,stop,num samples to generate)
# example single calculation
XYZ=sp.odeint(iron,XYZ0,t,args=(a12,a21,a23,a32,b13,b31,I))
### Sobol analysis ###
# defining problem
# can add the 'I' parameter
# assumed that range for each parameter is 80-120% of value assumed above
# can be changed
problem = {
'num_vars': 9, #a's, b's and initial condition
'names': ['a12', 'a21','a23','a32','b13','b31','x1_0','x2_0','x3_0'],
'bounds': np.column_stack((np.array([a12, a21,a23,a32,b13,b31,XYZ0[0],XYZ0[1],XYZ0[2]])*0.8,np.array([a12, a21,a23,a32,b13,b31,XYZ0[0],XYZ0[1],XYZ0[2]])*1.2))
}
# Generate samples
vals = saltelli.sample(problem, 500)
# initializing matrix to store output
Y = np.zeros([len(vals),1])
# Run model (example)
# numerically soves the ODE
# output is X1, X2, and X3 at the end time step
# could save output for all time steps if desired, but requires more memory
Y = np.zeros([len(vals),3])
for i in range(len(vals)):
Y[i][:] = sp.odeint(iron,[vals[i][6],vals[i][7],vals[i][8]],t,\
args=(vals[i][0],vals[i][1],vals[i][2],vals[i][3],vals[i][4],vals[i][5],I))[len(XYZ)-1]
# completing soboal analysis for each X1, X2, and X3
print('\n\n====X1 Sobol output====\n\n')
Si_X1 = sobol.analyze(problem, Y[:,0], print_to_console=True)
print('\n\n====X2 Sobol output====\n\n')
Si_X2 = sobol.analyze(problem, Y[:,1], print_to_console=True)
print('\n\n====X3 Sobol output====\n\n')
Si_X3 = sobol.analyze(problem, Y[:,2], print_to_console=True)
推荐阅读
- vba - 计算相同值之间有多少个单元格
- visual-studio-code - 保存格式文件上的格式更漂亮,但未保存更改
- nginx - TypeError: 'module' 对象在使用 Flask-SocketIO、uWSGI+gevent 和 Nginx 时不可调用。谁能帮我摆脱困境?
- c++ - 在getline之后在C ++中浏览一行内的单词
- php - return $pdf->stream 不加载徽标 - DomPDF Laravel
- sql - 选择和内连接返回重复项
- javascript - 如何返回对函数的引用而不是返回值
- plot - 从 InfuxDB 在 Grafana 中绘制单个传感器值
- java - 创建一个组合并行数组中的项目的流
- python - 如何创建具有多个标签和掩码的自定义 Pytorch 数据集?