analysis - 如何使用 SALib 对 ODE 执行灵敏度分析？

Question

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，初始条件)的影响。另外，我可以获得渐近解值吗？

Answer 1

以下代码是 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)