1

来回浏览文档，我能够在 Gekko中设置 动态参数估计。

这是代码，测量值如下所示（该文件MeasuredAlgebrProductionRate_30min_18h.csv在我的系统上命名，并;用作分隔符）：

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

#%% Read measurement data from CSV file
t_x_q_obs = np.genfromtxt('MeasuredAlgebrProductionRate_30min_18h.csv', delimiter=';')
#t_obs, x_obs, q_obs = t_xq_obs[:,0:3]

#%% Initialize Model
m = GEKKO(remote=False)
m.time = t_x_q_obs[:,0] #np.arange(0, 18/24+1e-6, 1/2*1/24)

# Declare parameter
V_liq   = m.Param(value = 159.0)

# Declare FVs
k_1     = m.FV(value = 0.80)
k_1.STATUS = 1
f_1     = m.FV(value = 10.0)
f_1.STATUS = 1

# Diff. Variables
X       = m.Var(value = 80.0) # at t=0
Y       = m.Var(value = 80.0*0.2)

rho_1   = m.Intermediate(k_1*X)
#q_prod  = m.Intermediate(0.52*f_1*X/24)
#X       = m.CV(value = t_x_q_obs[:,1])
q_prod  = m.CV(value = t_x_q_obs[:,2])

#%% Equations
m.Equations([X.dt() == -rho_1, Y.dt() == 0, q_prod == 0.52*f_1*X/24])

m.options.IMODE = 5
m.solve(disp=False)

#%% Plot some results
plt.plot(m.time, np.array(X.value)/10, label='X')
plt.plot(t_x_q_obs[:,0], t_x_q_obs[:,2], label='q_prod Meas.')
plt.plot(m.time, q_prod.value, label='q_prod Sim.')
plt.xlabel('time')
plt.ylabel('X / q_prod')
plt.grid()
plt.legend(loc='best')
plt.show()

0.0208333333 NaN 30.8306036
0.0416666667 NaN 29.1200832
0.0625 74.866 28.7700549
0.0833333333 NaN 29.2318865
0.104166667 NaN 30.7727362
0.125 NaN 29.8743804
0.145833333 NaN 29.9923447
0.166666667 NaN 30.9169679
0.1875 NaN 28.5956184
0.208333333 NaN 27.7361632
0.229166667 NaN 26.6669496
0.25 NaN 27.17477
0.270833333 75.751 23.6270346
0.291666667 NaN 23.0646928
0.3125 NaN 23.6442113
0.333333333 NaN 23.089118
0.354166667 NaN 22.9101616
0.375 NaN 22.7453854
0.395833333 NaN 23.2182759
0.416666667 NaN 21.4901903
0.4375 NaN 21.1449899
0.458333333 NaN 20.7093537
0.479166667 NaN 20.3109086
0.5 NaN 20.6825141
0.520833333 NaN 19.199583
0.541666667 NaN 19.6173416
0.5625 NaN 19.5543139
0.583333333 NaN 20.4501879
0.604166667 NaN 18.7678061
0.625 NaN 18.4629262
0.645833333 NaN 18.3730322
0.666666667 NaN 19.5375442
0.6875 NaN 18.1975297
0.708333333 NaN 18.0370627
0.729166667 NaN 17.5734727
0.75 NaN 18.8632046

到目前为止，一切都很好。假设我在某些时间点（第一列）也有 X（第二列）的测量值，其余的不可用（因此NaN）。我想调整k_1and f_1，以便模拟和观察到的变量尽可能地匹配X 。 q_prod

这对 Gekko 可行吗？如果是这样，怎么做？

m.time另一个问题：如果元素多于观察变量的时间点，Gekko 会抛出错误。X但是，我的和的初始值Y指的是t=0，而不是t=0.0208333333。因此，后面的注释掉的部分m.time =，见上文。（在 处t=0的测量不可用。）Gekko 中的初始条件是指 的第一个元素m.time，就像在 Matlab 中所做的那样，还是指t=0？

4

1 回答 1

3

如果您缺少测量值，则可以包含一个非数字值，例如 NaN，Gekko 会忽略目标函数中的该条目。这是一个测试用例，其中有一个 NaN 值ym：

具有 NaN 数据值的非线性回归

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([0.1,0.2,np.nan,0.5,0.8,2.0])
m = GEKKO(remote=False)
x = m.Param(value=xm,name='x')
a = m.FV()
a.STATUS=1
y = m.CV(value=ym,name='y')
y.FSTATUS=1
m.Equation(y==0.1*m.exp(a*x))
m.options.IMODE = 2
m.options.SOLVER = 1
m.solve(disp=True)
print('Optimized, a = ' + str(a.value[0]))
plt.plot(xm,ym,'bo')
plt.plot(xm,y.value,'r-')
m.open_folder()
plt.show()

当您打开运行文件夹m.open_folder()并查看数据文件时，值列gk_model0.csv中有 NaN 。y

y,x
0.1,0
0.2,1
nan,2
0.5,3
0.8,4
2.0,5

这IMODE=2是一个稳态回归问题，但显示了与动态估计问题相同的情况。有更多关于使用m.options.EV_TYPE=1（默认）或m.options.EV_TYPE=2估计的目标函数以及如何在数据文件中处理不良值的信息。当测量值是非数值时，该坏值将从目标函数求和中删除。这是一个带有动态模型的版本：

具有固定初始条件的动态回归

固定初始条件

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正如您所观察到的，您的测量值需要具有相同的长度。如果您缺少值，那么您可以将 a 添加np.nan到数据范围的开头。默认情况下，Gekko 使用value属性中指定的第一个值来设置初始条件。如果您不希望 Gekko 使用该值fixed_initial=False，请为您的CV.

具有自由初始条件的动态回归

自由初始条件

y = m.CV(value=ym,name='y',fixed_initial=False)
于 2019-08-23T19:04:39.737 回答