Блог

Я провожу анализ чувствительности с использованием инструментов Sobol и SALib на Python. В настоящее время мои параметры рандомизированы между нулем и единицей (я не уверен, какие значения параметров у меня еще будут, но они, вероятно, изменятся), и я использую solve_ivp для решения моей системы из 6 ОДУ. Затем я хочу использовать анализ Sobol в своем решении, чтобы показать, какие параметры оказывают наибольшее влияние.

Изначально я использовал ODEINT, и он вернул nan для всех моих параметров. Было рекомендовано использовать solve_ivp, чтобы избавиться от этой жесткости, но я не уверен, как правильно его инициализировать с помощью анализа Sobol.

Вы увидите мой код ниже.

 from scipy import integrate as sp
import numpy as np
import matplotlib.pyplot as plt
import SALib
from SALib.sample import saltelli
from SALib.analyze import sobol
from scipy.integrate import solve_ivp
%matplotlib inline

def M(C,phi,chi):
    return (phi*C)/(chi C)

def H(A,epsilon,eta):
    return epsilon ((eta*(A**n)/((A**n) 1)))

#Z[0]=X, Z[1]=Y, Z[2]=C1, Z[3]=C2, Z[4]=A1, Z[5]=A2

def dualchamber(time,Z,phi1,phi2,chi1,chi2,nu1,nu2,delta,kappa, 
                epsilon1,epsilon2,eta1,eta2,xi1,xi2,omega,rho,zeta,psi):
    dZdt = [(M(Z[2],phi1,chi1)*Z[0])-(nu1*Z[0]),
            (M(Z[3],phi2,chi2)*Z[1])-(nu2*Z[1]),
            (-M(Z[2],phi1,chi1)*Z[0]) delta-Z[2] (Z[3]/kappa),
            (-M(Z[3],phi2,chi2)*Z[1]) delta-Z[3] (kappa*Z[2]),
            (H(Z[4],epsilon1,eta1)*xi1*Z[0])-(omega*Z[4])-(rho*Z[4]) ((zeta/psi)*Z[5]),
            (H(Z[5],epsilon2,eta2)*xi2*Z[1])-(omega*Z[5])-(rho*Z[5]) (zeta*psi*Z[4])]
    return dZdt

time=np.linspace(0,100,100)

Zinitial = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1]

n=2.2
phi1 = np.random.uniform(0,1)
phi2 = np.random.uniform(0,1)
chi1 = np.random.uniform(0,1)
chi2 = np.random.uniform(0,1)
nu1 = np.random.uniform(0,1)
nu2 = np.random.uniform(0,1)
delta = np.random.uniform(0,1)
kappa = np.random.uniform(0,1)
epsilon1 = np.random.uniform(0,1)
epsilon2 = np.random.uniform(0,1)
eta1 = np.random.uniform(0,1)
eta2 = np.random.uniform(0,1)
xi1 = np.random.uniform(0,1)
xi2 = np.random.uniform(0,1)
omega = np.random.uniform(0,1)
rho = np.random.uniform(0,1)
zeta = np.random.uniform(0,1)
psi = np.random.uniform(0,1)
#print(phi1,phi2,chi1,chi2,nu1,nu2,delta,kappa,epsilon1,
      #epsilon2,eta1,eta2,xi1,xi2,omega,rho,zeta,psi)

sol = solve_ivp(lambda time, Z: dualchamber(time,Z,phi1,phi2,chi1,chi2,nu1,nu2,delta,kappa, 
                epsilon1,epsilon2,eta1,eta2,xi1,xi2,omega,rho,zeta,psi),
               [time[0],time[-1]], Zinitial, t_eval=time)

problem = {
  'num_vars': 24,
  'names': ['phi1', 'phi2','chi1','chi2','nu1','nu2','delta','kappa',
            'eps1', 'eps2', 'eta1', 'eta2', 'xi1', 'xi2','o', 
            'r', 'zeta', 'psi', 'X_0', 'Y_0', 'C1_0', 'C2_0', 'A1_0', 'A2_0'],
  'bounds': [[0.0001,1],[0.0001,1],[0.0001,1],[0.0001,1],
            [0.0001,1],[0.0001,1],[0.0001,1],[0.0001,1],
            [0.0001,1],[0.0001,1],[0.0001,1],[0.0001,1],
            [0.0001,1],[0.0001,1],[0.0001,1],[0.0001,1],
            [0.0001,1],[0.0001,1],[0.0001,1],[0.0001,1],
            [0.0001,1],[0.0001,1],[0.0001,1],[0.0001,1]]
}

vals = saltelli.sample(problem, 1000, calc_second_order=True)

Y = np.zeros([vals.shape[0]])

for i in range(len(vals)):
    Y[i][:] = solve_ivp(dualchamber(time,[vals[i][18],vals[i][19],vals[i][20],vals[i][21], vals[i][22],
      vals[i][23]],vals[i][0],vals[i][1],vals[i][2],vals[i][3],vals[i][4],vals[i][5], vals[i][6],
      vals[i][7],vals[i][8],vals[i][9],vals[i][10],vals[i][11], vals[i][12],vals[i][13],vals[i][14],vals[i][15],
      vals[i][16],vals[i][17]),[time[0],time[-1]], Zinitial)[len(sol)-1]

print('nn====X Sobol output====nn')
Si_X = sobol.analyze(problem, Y[:,0], print_to_console=True)
print('nn====Y Sobol output====nn')
Si_Y = sobol.analyze(problem, Y[:,1], print_to_console=True)
print('nn====C1 Sobol output====nn')
Si_C1 = sobol.analyze(problem, Y[:,2], print_to_console=True)
print('nn====C2 Sobol output====nn')
Si_C2 = sobol.analyze(problem, Y[:,3], print_to_console=True)
print('nn====A1 Sobol output====nn')
Si_A1 = sobol.analyze(problem, Y[:,4], print_to_console=True)
print('nn====A2 Sobol output====nn')
Si_A2 = sobol.analyze(problem, Y[:,5], print_to_console=True)
  

Любая помощь была бы отличной. Я действительно ценю любую обратную связь.
Спасибо.