如何定义时间相关的离散参数？

Question

最近用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，在 time = 2 时变为 2.5。但是我收到以下错误：

  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 是否需要更改，并相应地覆盖当前值。

Answer 1

您可以从 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 函数将问题转换为混合整数非线性规划形式，可能需要更长的时间才能解决。这是一个示例，其中 q_in=0 当 X>8 和 q_in=2.5 当 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()

还有一个 few other examples here 关于使用 Gekko 求解具有时变输入的 ODE。

如何定义时间相关的离散参数？

How to define time-dependent, discrete parameter?

python-3.x

differential-equations

gekko