Python 使用GEKKO将两个总体拟合到测量值如何优化第一个数据点

默认情况下，将gekko设置为初值问题，其中t=0处的值是固定的，并且方程不会在该初始点处求解.利用所计算的初始值仍然适合于该框架的一种方式是设置fixed_initial=False并且计算微分方程初始条件.然而，这并不能解决时间范围中第一个 node 的方程被停用的问题.

要克服这个问题，请try 在时间范围中设置另一个非常小的时间增量，时间步长非常小，例如1e-5.将分数f创建为FV类型，并设置STATUS=1以使用优化器计算该单个值.我在0.2和0.8之间加入了f的边界，只是为了让问题更有趣，它们可以改变.

# 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