如何在 ODE 右侧实现矩形脉冲（不连续）？

Question

除了之外，我还想在模拟时间内进行几次跳跃。正如@LutzL 所建议的，我尝试在一个循环中对每个阶段执行一个模拟，并使用 Connection 方法将它们拼接在一起，最终状态为 Phase 1 == initial states of Phase 2 等。但是我得到一个错误，说

Exception: @error: Model Expression *** Error in syntax of function string: Missing operator

Position: 6
11.55,11.55,11.55,11.55,11.55 ?

喂食次数和频率理解为：[t_start, t_end]，所以有两次喂食事件（分别从t = 2和t = 3.1开始）。一个需要 0.7，另一个需要 0.2（天）。所以有五个阶段(numOfPhases)：第一次喂食前，第一次喂食，第一次和第二次喂食之间，第二次喂食，第二次喂食后。

这是我的代码：

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
from gekko import GEKKO
# Construct input schedule
feedingTimes = np.array([2, 2.7, 3.1, 3.3])
feedingRates = np.array([60, 0, 430, 0])
tf      = 5
# Number of phases
numOfPhases = len(feedingTimes) + 1
# Construct time vectors (unequally spaced for now)
timeVector = np.zeros((numOfPhases, 10))
timeVector[0,:] = np.linspace(0, feedingTimes[0], 10)
for tI in np.arange(1,numOfPhases-1):
    timeVector[tI,:] = np.linspace(feedingTimes[tI-1], feedingTimes[tI], 10)
timeVector[-1,:] = np.linspace(feedingTimes[-1], tf, 10) 

#Initialize Model
m = GEKKO(remote=False)

q_in    = [m.Var(0.0, lb=-2, ub=500, fixed_initial=False) for i in range(numOfPhases)]

# (Example of) same parameter for each phase
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)
# Variables (one version of x for each phase)
X       = [m.Var(value = 11.55) for i in range(numOfPhases)]
Y       = [m.Var(value = 11.55*0.2) for i in range(numOfPhases)]

rho_1   = [m.Intermediate(k_1*X) for i in range(numOfPhases)]
q_prod  = [m.Intermediate(0.52*f_1*X) for i in range(numOfPhases)]

# Equations (different for each phase)
for pI in range(numOfPhases):
    m.time = timeVector[pI,:]
    m.Equations([X[pI].dt() == q_in[pI]/V_liq*(X_in - X[pI]) - rho_1[pI], \
                 Y[pI].dt() == q_in[pI]/V_liq*(Y_in - Y[pI])])

# Connect phases together at endpoints
for pI in range(numOfPhases-1):
    m.Connection(X[pI+1], X[pI], 1, len(m.time)-1, 1, 3)
    m.Connection(Y[pI+1], Y[pI], 1, len(m.time)-1, 1, 3)

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

plt.plot(m.time, q_in.value, label='q_in')
plt.plot(m.time, X.value, label='X')
plt.xlabel('time')
plt.ylabel('X')
plt.grid()
plt.legend(loc='best')
plt.show()

Answer 1

中间体 rho_1 和 q_prod 定义不正确，因为 X 是一个变量列表，而不是单个变量。

rho_1   = [m.Intermediate(k_1*X[i]) for i in range(numOfPhases)]
q_prod  = [m.Intermediate(0.52*f_1*X[i]) for i in range(numOfPhases)]

但是，我认为您的问题不需要 X 值的数组，因为您希望一段的最终条件等于下一段的初始条件并且您的方程式不会改变从一个部分到下一个部分。不允许 m.time 的多个定义。您只需为每个解决方案定义一次时间，并且每个段必须具有相同的时间。有一种方法可以对每个片段进行时间缩放，使其可变，但这是一个更高级的主题。此外，timeVector 不能有两个相等的连续值（必须递增）。我在每个喂食阶段的开始添加了 1e-5。

timeVector = np.zeros((numOfPhases, 10))
timeVector[0,:] = np.linspace(0, feedingTimes[0], 10)
for tI in np.arange(1,numOfPhases-1):
    timeVector[tI,:] = np.linspace(feedingTimes[tI-1]+1e-5, feedingTimes[tI], 10)
timeVector[-1,:] = np.linspace(feedingTimes[-1]+1e-5, tf, 10)

timeVector 是一个 (5,10) numpy 数组。您需要将其展平为 Gekko 的一维数组。我使用 numpy.reshape 函数进行了此修改。

m.time = np.reshape(timeVector, -1)

如果您需要针对不同时间段使用不同的方程式，那么您可能需要使用不同的方程式定义一个 X 值数组。因为您对每个部分使用相同的方程式，所以我建议使用顺序方法，如下所示。

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
from gekko import GEKKO
# Construct input schedule
feedingTimes = np.array([2, 2.7, 3.1, 3.3])
feedingRates = np.array([60, 0, 430, 0])
tf      = 5
# Number of phases
numOfPhases = len(feedingTimes) + 1
# Construct time vectors (unequally spaced for now)
timeVector = np.zeros((numOfPhases, 10))
timeVector[0,:] = np.linspace(0, feedingTimes[0], 10)
for tI in np.arange(1,numOfPhases-1):
    timeVector[tI,:] = np.linspace(feedingTimes[tI-1]+1e-5, feedingTimes[tI], 10)
timeVector[-1,:] = np.linspace(feedingTimes[-1]+1e-5, tf, 10)

feedVector = np.zeros(tf*10)
for i in range(len(feedingRates)):
    feedVector[(i+1)*10:(i+2)*10] = feedingRates[i]

#Initialize Model
m = GEKKO(remote=False)
m.time = np.reshape(timeVector, -1)

q_in    = m.Param(value=feedVector)

# (Example of) same parameter for each phase
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)
# Variables (one version of x for each phase)
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)

# Equations (different for each phase)
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=True)

plt.plot(m.time, q_in.value, label='q_in')
plt.plot(m.time, X.value, label='X')
plt.xlabel('time')
plt.grid()
plt.legend(loc='best')
plt.show()

如果您最终需要切换到优化/控制，并且每个段只需要一个进给率值，那么我建议参数 MV_STEP_HOR 为 described in the documentation。

如何在 ODE 右侧实现矩形脉冲（不连续）？

How to implement rectangular pulses (discontinuities) on ODE right-hand side?

python-3.x

differential-equations

gekko