3

最近用GEKKO搭建了一个小模型。它包含一个实际随时间变化的参数。我该如何实施？我尝试使用if3，但它给出了一个错误。

这是MWE：

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Started on 10-08-2019

@author: winkmal
"""
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
#Initialize Model
m = GEKKO(remote=False)

# Parameters
k_1     = m.Param(value = 0.19)
f_1     = m.Param(value = 29.0)
V_liq   = m.Param(value = 159.0)

q_in    = m.Param(value = 2.5)
X_in    = m.Param(value = 271.77)
Y_in    = m.Param(value = 164.34)

X       = m.Var(value = 11.55)
Y       = m.Var(value = 11.55*0.2)
rho_1   = m.Intermediate(k_1*X)
q_prod  = m.Intermediate(0.52*f_1*X)

m.time  = np.arange(0,5,1/12)

m.Equations([X.dt() == q_in/V_liq*(X_in - X) - rho_1, \
             Y.dt() == q_in/V_liq*(Y_in - Y)]) 
#Dynamic simulation
m.options.IMODE = 4
m.solve(disp=False)

plt.plot(m.time, X.value)
plt.xlabel('time')
plt.ylabel('X')
plt.show()

我尝试了以下方法：

q_in    = m.if3(m.time - 2, 0, 2.5)

所以q_in最初是 0，然后在 2.5 时变为 2.5 time = 2。但我收到以下错误：

  File "/usr/local/lib/python3.7/site-packages/gekko/gekko.py", line 1838, in solve
    raise Exception(apm_error)

Exception:  @error: Equation Definition
 Equation without an equality (=) or inequality (>,<)
 (((1-int_v5))*([-2.-1.91666667-1.83333333-1.75-1.66666667-1.58333333
 STOPPING...

你知道我怎么能做到这一点吗？实际上，这个变量在 0 到 60 之间跳跃了几次，我在 CSV 文件中有可用的时间点。理想情况下，我可以创建一个循环，在每次迭代时检查是否q_in需要更改，并相应地覆盖当前值。

4

1 回答 1

2

您可以从 CSV 读取输入，并将随时间变化的值分配给q_in.value参数初始化期间（参见示例 #1）或在每个时间积分间隔值更改的循环中（参见示例 #2）。示例 1 和 2 都产生以下结果，但示例 1 更快。

积分结果

m.options.IMODE=7如果您的时间跨度很长，示例 1 的选项也可能会更快。IMODE=7使用顺序求解方法而不是同时求解方法。

示例 1

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
t = np.arange(0,5,1/12)
step = [0 if z<2 else 2.5 for z in t]
m = GEKKO(remote=False)
k_1     = m.Param(value = 0.19)
f_1     = m.Param(value = 29.0)
V_liq   = m.Param(value = 159.0)
q_in    = m.Param(value = step)
X_in    = m.Param(value = 271.77)
Y_in    = m.Param(value = 164.34)
X       = m.Var(value = 11.55)
Y       = m.Var(value = 11.55*0.2)
rho_1   = m.Intermediate(k_1*X)
q_prod  = m.Intermediate(0.52*f_1*X)
m.time  = t
m.Equations([X.dt() == q_in/V_liq*(X_in - X) - rho_1, \
             Y.dt() == q_in/V_liq*(Y_in - Y)]) 
m.options.IMODE = 4
m.solve(disp=False)
plt.plot(m.time,q_in.value,label=r'$q_{in}$')
plt.plot(m.time, X.value,label='X')
plt.plot(m.time, Y.value,label='Y')
plt.legend()
plt.xlabel('time')
plt.show()

示例 2

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
t = np.arange(0,5,1/12)
m = GEKKO(remote=False)
k_1     = m.Param(value = 0.19)
f_1     = m.Param(value = 29.0)
V_liq   = m.Param(value = 159.0)
q_in    = m.Param()
X_in    = m.Param(value = 271.77)
Y_in    = m.Param(value = 164.34)
X       = m.Var(value = 11.55)
Y       = m.Var(value = 11.55*0.2)
rho_1   = m.Intermediate(k_1*X)
q_prod  = m.Intermediate(0.52*f_1*X)
m.time  = [t[0],t[1]]
m.Equations([X.dt() == q_in/V_liq*(X_in - X) - rho_1, \
             Y.dt() == q_in/V_liq*(Y_in - Y)]) 
m.options.IMODE = 4
# store Xs and Ys for plotting
for i in range (1,len(t)):
    q_in.value = 0 if t[i]<2 else 2.5
    m.solve(disp=False)
    if i==1:
        Xs = [X.value[0]]
        Ys = [Y.value[0]]
    Xs.append(X.value[1])
    Ys.append(Y.value[1])
step = [0 if z<2 else 2.5 for z in t]
plt.plot(t,step,label=r'$q_{in}$')
plt.plot(t, Xs,label='X')
plt.plot(t, Ys,label='Y')
plt.legend()
plt.xlabel('time')
plt.show()

如果您需要q_in依赖某些变量的值，那么您可以使用该m.if3函数。然而，这是一个更难解决的问题，因为该m.if3函数将问题转换为可能需要更长时间才能解决的混合整数非线性规划形式。这是一个示例 where q_in=0whenX>8和q_in=2.5when X<=8。但是，它并没有对我收敛。我不知道为什么，我需要做一些额外的挖掘，但我虽然你想拥有它，以防它对你有用。

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=False)
k_1     = m.Param(value = 0.19)
f_1     = m.Param(value = 29.0)
V_liq   = m.Param(value = 159.0)
X_in    = m.Param(value = 271.77)
Y_in    = m.Param(value = 164.34)
X       = m.Var(value = 11.55,name='X')
Y       = m.Var(value = 11.55*0.2,name='Y')
rho_1   = m.Intermediate(k_1*X)
q_prod  = m.Intermediate(0.52*f_1*X)
q_in    = m.if3(8-X, 0.0, 2.5)
m.time  = np.arange(0,5,1/12)
m.Equations([X.dt() == q_in/V_liq*(X_in - X) - rho_1, \
             Y.dt() == q_in/V_liq*(Y_in - Y)]) 
m.options.IMODE = 6
m.options.SOLVER = 1
m.solve(disp=True)
plt.plot(m.time,q_in.value,label=r'$q_{in}$')
plt.plot(m.time, X.value,label='X')
plt.plot(m.time, Y.value,label='Y')
plt.legend()
plt.xlabel('time')
plt.show()

这里还有一些关于使用 Gekko 求解具有时变输入的 ODE的其他示例。

于 2019-08-12T14:16:15.513 回答