Оптимизация трех переменных в дифференциальном уравнении на основе доступных данных из решения дифференциального уравнен ⇐ Python

[Anonymous]

У меня есть следующие точки данных (я называю их фактическими точками данных):

Код: Выделить всё

y_data = np.array([0, 32.1463583, 33.1915926, 37.9100309, 39.2501778, 40.8225707, 48])
t_data = np.array([0, 26.75, 72.25, 163.4166667, 209.25, 525, 1250])

и у меня есть дифференциальное уравнение, включающее y и t:

Код: Выделить всё

dy/dt=(1-y)/((a+b*t)*exp(-E/3060.8))

моя цель — оптимизировать a, b и E так, чтобы решение приведенного выше дифференциального уравнения лучше всего соответствовало моим фактическим данным (y_data).
для этого. У меня в голове были следующие шаги:

Код: Выделить всё

1- set initial guess for a=1, b=20, E=10000
2- create a for loop ( for 1000 iterations)
3- solve differential equation using ODEint with initial guess
4- Find difference between calculated value of y from solution of differential equation  at time equal to "t" with actual "y" value
5- print error
6- update a to a , b, E and repeat from step 3

и вот мой код:

Код: Выделить всё

import numpy as np
from scipy.integrate import odeint

# Given data
y_data = np.array([0, 32.1463583, 33.1915926, 37.9100309, 39.2501778, 40.8225707, 48])
t_data = np.array([0, 26.75, 72.25, 163.4166667, 209.25, 525, 1250])

# Initial guess for parameters
a = 1
b = 20
E = 10000

# Number of iterations
iterations = 1000

# Tolerance for convergence
tolerance = 1e-6

# Step size for numerical gradient approximation
h = 1e-6

# Perform optimization
for i in range(iterations):
# Calculate gradients using numerical gradient approximation
def calculate_error(a, b, E):
def model(y, t):
return (1 - y) / ((a + b * t) * np.exp(-E / 3060.8))

y_solution = odeint(model, y_data[0], t_data)
error = np.mean(np.abs(y_solution[:, 0] - y_data))
return error
print (error)

grad_a = (calculate_error(a + h, b, E) - calculate_error(a, b, E)) / h
grad_b = (calculate_error(a, b + h, E) - calculate_error(a, b, E)) / h
grad_E = (calculate_error(a, b, E + h) - calculate_error(a, b, E)) / h

# Update parameters
a -= 0.1 * grad_a
b -= 0.1 * grad_b
E -= 0.1 * grad_E

# Check for convergence
if abs(grad_a) < tolerance and abs(grad_b) < tolerance and abs(grad_E) < tolerance:
break

# Output the optimized parameters
print("Optimized Parameters (a, b, E):", a, b, E)

Проблема в том, что этот код не работает и постоянно выдает ошибку «32.223947830786685». Это каким-то образом указывает на то, что функция обновления для обновления a, b и E работает неправильно.
Для справки, я использовал функцию для

Код: Выделить всё

(calculate_error(a + h, b, E) - calculate_error(a, b, E))

чтобы вручную рассчитать градиент a, b, E.
Есть какие-нибудь предложения по устранению моей проблемы? или какое-либо альтернативное решение для поиска параметров ставок a, b, E?