Используйте Gekko и Python, чтобы согласовать численное решение ODE с данными - PullRequest
Новая
       4

Используйте Gekko и Python, чтобы согласовать численное решение ODE с данными

4 голосов
/ 09 марта 2020

Используйте Gekko для подбора числового решения ОДУ для данных.

Привет всем! Мне было интересно, если это возможно, чтобы соответствовать коэффициентам ODE с помощью GEKKO. Я безуспешно пытался повторить приведенный здесь пример .

. Это то, что я придумал (но он ошибочен - и я должен упомянуть, что мои математические навыки, к сожалению, довольно плохие): 

import numpy as np
from gekko import GEKKO

tspan = [0, 0.1, 0.2, 0.4, 0.8, 1]
Ca_data = [2.0081,  1.5512,  1.1903,  0.7160,  0.2562,  0.1495]

m = GEKKO(remote=False)

t = m.Param(value=tspan)
m.time = t
Ca_m = m.Param(value=Ca_data)

Ca = m.Var()

k = m.FV(value=1.3)
k.STATUS = 1

m.Equation( Ca.dt() == -k * Ca)

m.Obj( ((Ca-Ca_m)**2)/Ca_m )

m.options.IMODE = 2
m.solve(disp=True)
print(k.value[0]) #2.58893455 is the solution

Может кто-нибудь помочь мне здесь? Большое спасибо, Мартин

(Это мой первый пост здесь - пожалуйста, будьте осторожны, если я сделал что-то неподходящее.)

1 Ответ

2 голосов
/ 10 марта 2020

Ваше решение было близко, но вам нужно было:

  • Больше узлов (по умолчанию = 2) для повышения точности. Гекко добавляет только те точки, которые вы определяете. См. дополнительную информацию о расположении .
  • . Определите Ca как m.CV() для использования встроенной модели ошибок вместо m.Var() и m.Obj с NODES> = 3. В противном случае внутренние узлы каждого интервала размещения также сопоставляются с измерениями, и это дает немного неправильный ответ.
  • Установите EV_TYPE=2 для использования квадратичной ошибки. Целевое абсолютное значение EV_TYPE=1 (по умолчанию) дает правильный, но немного другой ответ. 
import numpy as np
from gekko import GEKKO

m = GEKKO(remote=False)
m.time = [0, 0.1, 0.2, 0.4, 0.8, 1]
Ca_data = [2.0081,  1.5512,  1.1903, 0.7160,  0.2562,  0.1495]
Ca = m.CV(value=Ca_data); Ca.FSTATUS = 1 # fit to measurement
k = m.FV(value=1.3); k.STATUS = 1        # adjustable parameter
m.Equation(Ca.dt()== -k * Ca)            # differential equation

m.options.IMODE = 5   # dynamic estimation
m.options.NODES = 5   # collocation nodes
m.options.EV_TYPE = 2 # squared error
m.solve(disp=True)    # display solver output
print(k.value[0])     # 2.58893455 is the curve_fit solution

Решение - k=2.5889717102. График показывает соответствие измеренным значениям. 

import matplotlib.pyplot as plt  # plot solution
plt.plot(m.time,Ca_data,'ro')
plt.plot(m.time,Ca.value,'bx')
plt.show()

Parameter Estimation Solution

Есть дополнительные учебные пособия и материал курса при оценке параметров с помощью дифференциальных и алгебраических моделей c.

...