Python Gekko ODE扁平结果

Question

我一直在和Gekko(https://gekko.readthedocs.io/en/latest/)玩。为了测试它，我从这里的https://julia.quantecon.org/continuous_time/seir_model.html实现了Covid模型。但结果只是平面图。它看起来一点也不像链接的结果。有没有人看到我做错了什么？

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

m = GEKKO(remote=False)    # create GEKKO model

# constants
y = 1/18
o = 1/5.2
R = 3

# create variables
i_0 = 1E-7
e_0 = 4.0 * i_0
s_0 = 1.0 - i_0 - e_0
r_0 = 0.0

s = m.Var(s_0)
e = m.Var(e_0)
i = m.Var(i_0)
r = m.Var(r_0)

# create equations
m.Equation(s.dt()==-y*R*s*i)
m.Equation(e.dt()==-y*R*s*i - o*e)
m.Equation(i.dt()==o*e - y*i)
m.Equation(r.dt()==y*i)

m.time = np.linspace(0,350)

# solve ODE
m.options.IMODE = 4
m.solve()

# plot results
plt.plot(m.time, s, label="s")
plt.plot(m.time, e, label="e")
plt.plot(m.time, i, label="i")
plt.plot(m.time, r, label="r")
plt.legend()
plt.show()

Answer 1

这是一个错误的等式：

m.Equation(e.dt()==-y*R*s*i - o*e)

它应该是：

m.Equation(e.dt()==y*R*s*i - o*e)

你还需要这些选项来提高精度，因为问题的值很小，Gekko默认使用1e-6容差：

# solve ODE
m.options.IMODE = 7
m.options.OTOL  = 1e-8
m.options.RTOL  = 1e-8
m.options.NODES = 3

如果你不想细化求解器容差，那么你也可以缩放方程。

sc = 1000
m.Equation(sc*s.dt()==-y*R*s*i*sc)
m.Equation(sc*e.dt()==(y*R*s*i - o*e)*sc)
m.Equation(sc*i.dt()==(o*e - y*i)*sc)
m.Equation(sc*r.dt()==y*i*sc)

以下是该脚本的工作版本，并进行了一些其他修改，以提高准确性和求解速度：

​

​

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

m = GEKKO(remote=False)    # create GEKKO model

# constants
y = 1/18
o = 1/5.2
R = 3

# create variables
i_0 = 1E-7
e_0 = 4.0 * i_0
s_0 = 1.0 - i_0 - e_0
r_0 = 0.0

s = m.Var(s_0)
e = m.Var(e_0)
i = m.Var(i_0)
r = m.Var(r_0)

# create equations
m.Equation(s.dt()==-y*R*s*i)
m.Equation(e.dt()==y*R*s*i - o*e)
m.Equation(i.dt()==o*e - y*i)
m.Equation(r.dt()==y*i)

m.time = np.linspace(0,350)
m.time = np.insert(m.time,1,np.logspace(-8,-1,10))
    
# solve ODE
m.options.IMODE = 7
m.options.OTOL  = 1e-8
m.options.RTOL  = 1e-8
m.options.NODES = 3
m.solve()

# plot results
plt.plot(m.time, s, label="s")
plt.plot(m.time, e, label="e")
plt.plot(m.time, i, label="i")
plt.plot(m.time, r, label="r")
plt.legend()
plt.show()

这里是另一个SEIR model for COVID-19 in Gekko，以及一个关于医疗保健资源优化的简单案例研究。

​

​

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

t_incubation = 5.1 
t_infective = 3.3 
R0 = 2.4 
N = 100000 

# initial number of infected and recovered individuals
e_initial = 1/N
i_initial = 0.00
r_initial = 0.00
s_initial = 1 - e_initial - i_initial - r_initial

alpha = 1/t_incubation
gamma = 1/t_infective
beta = R0*gamma

m = GEKKO()
u = m.FV(0)
s,e,i,r = m.Array(m.Var,4)
s.value = s_initial
e.value = e_initial
i.value = i_initial
s.value = s_initial
m.Equations([s.dt()==-(1-u)*beta * s * i,\
             e.dt()== (1-u)*beta * s * i - alpha * e,\
             i.dt()==alpha * e - gamma * i,\
             r.dt()==gamma*i])

t = np.linspace(0, 200, 101)
t = np.insert(t,1,[0.001,0.002,0.004,0.008,0.02,0.04,0.08,\
                   0.2,0.4,0.8])
m.time = t
m.options.IMODE=7
m.solve(disp=False)

# plot the data
plt.figure(figsize=(8,5))
plt.subplot(2,1,1)
plt.plot(m.time, s.value, color='blue', lw=3, label='Susceptible')
plt.plot(m.time, r.value, color='red',  lw=3, label='Recovered')
plt.ylabel('Fraction')
plt.legend()

plt.subplot(2,1,2)
plt.plot(m.time, i.value, color='orange', lw=3, label='Infective')
plt.plot(m.time, e.value, color='purple', lw=3, label='Exposed')
plt.ylim(0, 0.2)
plt.xlabel('Time (days)')
plt.ylabel('Fraction')
plt.legend()

plt.show()

该案例研究展示了如何计算最佳交互参数以满足医疗保健约束。

​

​