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

Posted 2023-02-18

【中文标题】如何定义时间相关的离散参数？

【英文标题】：How to define time-dependent, discrete parameter?

【发布时间】：2019-12-18 00:17:27

【问题描述】：

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

【问题讨论】：

你能找到错误的确切位置吗？从文档中我得到的印象是if3 在这里被错误应用，因为t 不是模型的静态参数。使用q_in = lambda t: 0 if t<2 else 2.5 会稍微好一些，如果您可以拆分积分区间，则更好，这样第二段的 IC 是第一段的最后一个值。

gekko.readthedocs.io/en/latest/imode.html#dynamics ODE 求解器不仅希望 ODE 中的右侧是连续的，而且是平滑的。您确实在使用超出规格的包。可行的是使用 ODE à la gekko.readthedocs.io/en/latest/… 的多个实例，其中第一个实例具有 q_in=0，第二个实例具有 q_in=2.5，第二个实例的状态通过 m.Connection 链接到时间索引处的第一个实例时间t=2.

【参考方案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=0when X>8 和 q_in=2.5 when 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。

【讨论】：

实际上，我可以让if3 对我的代码进行以下小改动：m.time = np.arange(0,5,1/12); t = m.Param(value = m.time); ... q_in = m.Var(); q_in = m.if3(t - 2, 0, 2.5);

在过多地编辑 OP 之前，我会将其标记为已解决（可能在 t 期间进行 一个 跳转）并发布一个处理多个跳转的新问题，因为我还不能让它工作。