使用 GEKKO 将两个群体拟合到测量结果如何优化第一个数据点

Question

我有这段代码来尝试将两个总体的总和拟合到测量数据系列。

dS/dt = (a - b) * S
dR/dt = (a - b - c)* R
X(t) = S(t) + R(t)

我还想优化第一点，它现在是固定的，我曾想过使用一个名为 f 的新变量，它位于

0

和

1

之间，代表人口

1

和人口

2

。 像这样的东西：

S0 = x_meas[0] * f
R0 = x_meas[0] * (1-f)

这是代码：

from gekko import GEKKO
m = GEKKO(remote=False)

# data
m.time = [0,0.1,0.2,0.3,0.4,0.5,0.8,0.85,0.9,0.95,1]
x_meas = [1.1,1.2,1.56,1.77,1.89,2.1,2.3,2.5,2.7,2.2,2.1]

# parameters, single value over time horizon
p = m.Array(m.FV,3,lb=0.5,ub=10); a,b,c = p
# turn STATUS On for adjustment by solver
for fv in p:
   fv.STATUS=1

f0 = m.FV(value=0.01, lb=0, ub=1)
f0.STATUS=1
# variables
S,R,X = m.Array(m.Var,3,lb=0)
# change initial conditions
S0,R0 = [0.35,0.65]; X0 = S0+R0
S.value=S0; R.value=R0; X.value=X0

# measured values
Xm = m.Param(x_meas)

# equations
m.Equations([S.dt()==(a-b)*S,
             R.dt()==(a-b-c)*R,
             X==S+R])

# minimize difference to measurements
m.Minimize((X-Xm)**2)

m.options.IMODE = 5   # dynamic estimation
m.options.NODES = 3   # collocation nodes
m.solve(disp=False)   # display solver output

print(f'a: {a.value[0]}, b: {b.value[0]}, c: {c.value[0]}')

import matplotlib.pyplot as plt
plt.figure(figsize=(7,3))
plt.plot(m.time,S.value,'b:',linewidth=3,label='S')
plt.plot(m.time,R.value,'r--',linewidth=2,label='R')
plt.plot(m.time,X.value,'k--',label='X')
plt.plot(m.time,Xm.value,'kx',label='Xm')
plt.legend(); plt.grid(); plt.xlabel('Time')
plt.tight_layout()
plt.savefig('fit.png',dpi=300)
plt.show()

Answer 1

默认情况下，

gekko

设置为初始值问题，其中

t=0

处的值是固定的，并且方程在该初始点未求解。使用计算的初始值仍然适合该框架的一种方法是设置

fixed_initial=False

并计算微分方程初始条件。但是，这并不能解决时间范围内第一个节点的方程被停用的问题。

要克服这个问题，请尝试在时间范围内以非常小的时间步长设置另一个非常小的时间增量，例如

1e-5

。创建分数

f

作为

FV

类型，并设置

STATUS=1

以使用优化器计算该单个值。我在

f

和

0.2

之间添加了

0.8

的界限，只是为了让问题更有趣，并且它们可以更改。

# initial condition variables
f = m.FV(value=0.5,lb=0.2,ub=0.8)
f.STATUS = 1

仅在该小时间步长上激活

R

、

S

和

f

之间关系的方程。

# activate only for initial condition
init = np.zeros(len(m.time)); init[1]=1
init = m.Param(init)
m.Equations([init*(S-x_meas[0]*f)==0,
             init*(R-x_meas[0]*(1-f))==0])

问题的其余部分保持不变。

from gekko import GEKKO
import numpy as np
m = GEKKO(remote=True)

# data
m.time = [0,1e-5,0.1,0.2,0.3,0.4,0.5,0.8,0.85,0.9,0.95,1]
x_meas = [1.1,1.1,1.2,1.56,1.77,1.89,2.1,2.3,2.5,2.7,2.2,2.1]

# parameters, single value over time horizon
p = m.Array(m.FV,3,lb=0.5,ub=10); a,b,c = p
# turn STATUS On for adjustment by solver
for fv in p:
   fv.STATUS=1

# variables
S,R,X = m.Array(m.Var,3,value=0.5,lb=0,fixed_initial=False)

# initial condition variables
f = m.FV(value=0.5,lb=0.2,ub=0.8)
f.STATUS = 1

# activate only for initial condition
init = np.zeros(len(m.time)); init[1]=1
init = m.Param(init)
m.Equations([init*(S-x_meas[0]*f)==0,
             init*(R-x_meas[0]*(1-f))==0])

# change initial condition for X
X.value=x_meas[0]

# measured values
Xm = m.Param(x_meas)

# equations
m.Equations([S.dt()==(a-b)*S,
             R.dt()==(a-b-c)*R,
             X==S+R])

# minimize difference to measurements
m.Minimize((X-Xm)**2)

m.options.IMODE = 5   # dynamic estimation
m.options.NODES = 2   # collocation nodes
m.options.SOLVER = 1  # APOPT solver
m.solve(disp=True)    # display solver output

print(f'a: {a.value[0]}, b: {b.value[0]}, c: {c.value[0]}')
print(f'S0: {S.value[0]}, R0: {R.value[0]}, f:{f.value[0]}')

import matplotlib.pyplot as plt
plt.figure(figsize=(7,3))
plt.plot(m.time,S.value,'b:',linewidth=3,label='S')
plt.plot(m.time,R.value,'r--',linewidth=2,label='R')
plt.plot(m.time,X.value,'k--',label='X')
plt.plot(m.time,Xm.value,'kx',label='Xm')
plt.legend(); plt.grid(); plt.xlabel('Time')
plt.tight_layout()
plt.savefig('fit.png',dpi=300)
plt.show()

它成功解决了：

 Iter    Objective  Convergence
   20  1.01871E+00  1.15106E-04
   21  1.01863E+00  1.02840E-04
   22  1.01863E+00  3.60679E-09
   23  1.01863E+00  5.74752E-10
   24  1.01863E+00  5.74752E-10
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   6.410000000323635E-002 sec
 Objective      :    1.01863407387195     
 Successful solution
 ---------------------------------------------------
 
a: 1.6769964544, b: 0.50197865393, c: 0.5
S0: 0.21999741496, R0: 0.87999405984, f:0.2