使用 GEKKO 将两个群体拟合到测量结果

Question

我需要使用 GEKKO 将取决于 3 个参数的两个不同偏方程定义的两个总体的总和拟合到类型（测量、时间）的测量列表中。 微分方程为：

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

（a,b,c为参数）

我需要将 X 拟合到测量值

Answer 1

以下是样本测量的

X

值：

# 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]

Gekko 调整参数

a

、

b

和

c

。

# 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

X

的值通过

S

和

R

的微分方程以及与

X=S+R

的代数方程计算。

# 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

测量的

X

值与预测的

X

值之间的平方误差最小化。

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

需要一些配置参数才能在动态估计模式下运行 (

IMODE=5

)，并且每个时间步长有 3 个配置节点 (

NODES=3

)。

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

这会产生一个优化的解决方案，使平方误差最小化，但不是唯一的，因为

a

和

b

是共线参数。求解器可以增加

a

和

b

并得到相同的答案。

这是完整的脚本：

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
   
# 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()

其他示例问题可从动态优化课程获得。