python - 通过拟合实验数据估计复微分方程组参数时如何增加预测数据？

Question

我正在尝试在我的程序中估计我的 ODE 集的参数，始终最小化我的实验数据和预测数据之间的误差。

问题是，我可以获得很好的预测和非常好的拟合，但我只能估计与我的实验数据相同数量的点，这会返回一个非常奇怪的输出。

能否请您提供更多有关如何获得更准确的预测点集的信息？

代码可以在下面找到。

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
xm = np.array([0,1,2,3,4,5])
ym = np.array([2.0,1.5,np.nan,2.2,3.0,5.0])
m = GEKKO(remote=False)
m.time = xm
a = m.FV(lb=0.1,ub=2.0)
a.STATUS=1
y = m.CV(value=ym,name='y',fixed_initial=False)
y.FSTATUS=1
m.Equation(y.dt()==a*y)
m.options.IMODE = 5
m.options.SOLVER = 1
m.solve(disp=True)
print('Optimized, a = ' + str(a.value[0]))
plt.figure(figsize=(6,2))
plt.plot(xm,ym,'bo',label='Meas')
plt.plot(xm,y.value,'r-',label='Pred')
plt.ylabel('y')
plt.ylim([0,6])
plt.legend()
plt.show()

如果我替换变量 m.time 以获得更多预测数据:
m.time = np.linspace(0,5,30)我收到错误:引发异常('数据数组必须具有相同的长度，并匹配动态问题中的时间离散化')异常(exception):数据数组必须具有相同的长度，并匹配动态问题中的时间离散化

Answer 1

下面显示了两个选项(方法 1 和 2)。您可以 plot the interpolating nodes为您提供更高的分辨率或创建新的模拟模型。

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
xm = np.array([0,1,2,3,4,5])
ym = np.array([2.0,1.5,np.nan,2.2,3.0,5.0])
m = GEKKO(remote=False)
m.time = xm
a = m.FV(lb=0.1,ub=2.0)
a.STATUS=1
y = m.CV(value=ym,name='y',fixed_initial=False)
y.FSTATUS=1
m.Equation(y.dt()==a*y)
m.options.IMODE = 5
m.options.SOLVER = 1
m.options.CSV_WRITE = 2  # For Method 1
m.options.NODES = 3      # For Method 1 (options 3-6)
m.solve(disp=True)
print('Optimized, a = ' + str(a.value[0]))

# Method 1: Plot interpolating nodes
import json
with open(m.path+'//results_all.json') as f:
    results = json.load(f)

# Method 2: Re-simulate with more points
sim = GEKKO(remote=False)
ap = a.value[0]
xp = np.linspace(0,7); sim.time=xp
yp = sim.Var(y.value[0])
sim.Equation(yp.dt()==ap*yp)
sim.options.NODES = 3
sim.options.IMODE=4; sim.solve()

plt.figure(figsize=(6,2))
plt.plot(xm,ym,'bo',label='Meas')
plt.plot(xm,y.value,'gs-.',label='Pred Original')
plt.plot(results['time'],results['y'],'kx-',\
         MarkerSize=10,label='Pred Method 1')
plt.plot(xp,yp,'r.--',label='Pred Method 2')
plt.ylabel('y')
plt.ylim([0,10])
plt.legend()
plt.show()

第三种选择是重置原始模型的 .value，但这可能很乏味。相反，您也可以像使用 example of Moving Horizon Estimation and Model Predictive Control 一样在循环中创建估计模型和仿真模型。使用相同模型但在它们之间传递参数:

# use remote=True for MacOS
mhe = GEKKO(name='tclab-mhe',remote=False)
mpc = GEKKO(name='tclab-mpc',remote=False)

# create 2 models (MHE and MPC) in one loop
for m in [mhe,mpc]:
    # Parameters with bounds
    m.K1 = m.FV(value=0.607,lb=0.1,ub=1.0)
    m.K2 = m.FV(value=0.293,lb=0.1,ub=1.0)
    m.K3 = m.FV(value=0.24,lb=0.1,ub=1.0)
    m.tau12 = m.FV(value=192,lb=100,ub=200)
    m.tau3 = m.FV(value=15,lb=10,ub=20)
    m.Ta = m.Param(value=23.0) # degC

    m.Q1 = m.MV(value=0,lb=0,ub=100,name='q1')
    m.Q2 = m.MV(value=0,lb=0,ub=100,name='q2')

    # Heater temperatures
    m.TH1 = m.SV(value=T1m[0])
    m.TH2 = m.SV(value=T2m[0])
    # Sensor temperatures
    m.TC1 = m.CV(value=T1m[0],name='tc1')
    m.TC2 = m.CV(value=T2m[0],name='tc2')

    # Temperature difference between two heaters
    m.DT = m.Intermediate(m.TH2-m.TH1)

    # Equations
    m.Equation(m.tau12*m.TH1.dt()+(m.TH1-m.Ta)==m.K1*m.Q1+m.K3*m.DT)
    m.Equation(m.tau12*m.TH2.dt()+(m.TH2-m.Ta)==m.K2*m.Q2-m.K3*m.DT)
    m.Equation(m.tau3*m.TC1.dt()+m.TC1==m.TH1)
    m.Equation(m.tau3*m.TC2.dt()+m.TC2==m.TH2)