我正在尝试在 python 中解决模型并将模型的未知参数拟合到实验数据。该模型由 2 个 ODE 组成,我使用 scipy.integrate.solve_ivp 对其进行求解。模型的参数是未知的,因此,我想使用多种方法来拟合它们。我的第一选择是差分进化,因为它可以为更复杂的模型提供非常好的结果(之前将它用作另一个包的一部分)。然而,问题是,当我给出我的模型(我希望它找到计算点和实验点之间的最小二乘的全局最小值)时,它会找到模型变得不稳定的参数(LSODA 无法整合它,因为内部步长为 0)。我试图捕捉 LSODA 在这些情况下抛出的运行时警告,但这并没有帮助。
我的模型和代码如下:
def aggregation_model(time, z, k1, k2, k3, k_1, k_2):
GPVI_0 = 55.5
GPVI, GPVI_Clust = z
dGPVIdt = - k1 * GPVI_Clust * GPVI + k_1 * (GPVI_0 - GPVI) - 2 * k2 * GPVI * GPVI
dGPVI_Clustdt = - (k_2 * GPVI_Clust + k3) * GPVI_Clust + k2 * GPVI * GPVI
return [dGPVIdt, dGPVI_Clustdt]
def res_squares(parameters):
time = np.linspace(0, 300, 1000)
timepoints = [0, 25, 50, 75, 100, 125, 150, 300]
val = [0, 2, 5, 10, 15, 20, 20, 18]
model_calc = solve_ivp(aggregation_model, [0, 300], [55.5, 0], args=parameters, max_step=100000,
dense_output=True, method='LSODA', rtol=1e-12, atol=1e-6)
solution = model_calc.sol(time)
transposed = list(map(list, zip(*solution.T)))
size = [0, ]
for i in range(1, len(transposed[0])):
size.append((55.5 - transposed[0][i]) / transposed[1][i])
diff = []
for i in range(len(timepoints)):
indexval = min(range(len(time)), key=lambda j: abs(time[j] - timepoints[i]))
diff.append((val[i] - size[indexval])**2)
# print(np.sqrt(np.sum(diff)))
output = np.sum(diff)
return output
if __name__ == '__main__':
from scipy.integrate import solve_ivp
from scipy.optimize import differential_evolution
import numpy as np
parameterBounds = []
parameterBounds.append([-1000000, 1000000]) # parameter bounds for a
parameterBounds.append([-1000000, 1000000]) # parameter bounds for a
parameterBounds.append([-1000000, 1000000]) # parameter bounds for a
parameterBounds.append([-1000000, 1000000]) # parameter bounds for a
parameterBounds.append([-1000000, 1000000]) # parameter bounds for a
result = differential_evolution(res_squares, parameterBounds, seed=3, disp=True, init='random')
print(result.x)
sol = solve_ivp(aggregation_model, [0, 300], [55.5, 0], args=parvalues,
dense_output=True, method='LSODA', rtol=1e-6, atol=1e-12)