1

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

异常：@error：模型表达式 *** 函数字符串语法错误：缺少运算符

位置：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()
4

1 回答 1

2

中间体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中描述的参数。

于 2019-08-17T13:40:33.197 回答