0

我正在使用 PyOMO 来模拟半批量反应。

考虑一个描述半间歇式反应器的 ODE 系统，其中一种反应物在t1单位时间内以给定的体积流量进料，反应一直持续到t end，并且显然t1 < t end。

要指定流程中的停止点，我可以使用条件规则（假设t1 = 3.5*60）：

def _vol_flow_in_schedule(mod,t):
 if t<=3.5*60:
  return mod.vol_flow_in[t] == (12.3/1000)/(3.5*60)
 else:
  return mod.vol_flow_in[t] == 0
m1.vol_flow_in_schedule = Constraint(m1.time,rule=_vol_flow_in_schedule)

这将产生不连续性（然后我的模型不会收敛）。我想要做的是使用一个 sigmoidal 函数，它将流转换为零而没有不连续性。

尽管我需要参考时间变量本身，但要实现 sigmoidal。

下面的 MATLAB 代码给了我想要的结果：

t=[0:1:500];
acc=2;  %Acceleration parameter, higher values yields sharper change.
time_of_step=3.5*60;
init_value = (12.3/1000)/(3.5*60);
end_value = 0;
sigmoidal=(init_value+(end_value-init_value)/2)...
          +((end_value-init_value)/2)*atan((t-time_of_step)*acc)/atan(max(t));

但是，此实现需要在函数中显式地使用时间变量。如何访问 PyOMO 规则中的时间变量？我尝试了以下方法，但出现“无法将标量组件 't_of_step' 视为数组”错误：

m1.init_value = Param(initialize = (12.3/1000)/(3.5*60))
m1.end_value  = Param(initialize = 0)
m1.t_of_step  = Param(initialize = 210)
m1.acc        = Param(initialize = 5)
.
.

def _vol_flow_sigmoidal (mod,t):
 return mod.vol_flow_in[t] == (mod.init_value+(mod.end_value-mod.init_value)/2)+((mod.end_value-mod.init_value)/2)*atan((t-mod.t_of_step)*mod.acc)/atan(1500)
m1.vol_flow_sigmoidal = Constraint(m1.time,rule=_vol_flow_sigmoidal)

希望我已经清楚地描述了我所追求的。欢迎任何提示，

谢谢！ 萨尔

4

1 回答 1

1

你如何声明m1.time索引？

我的猜测是您正在使用 NumPy 数组来初始化m1.time索引。Pyomo 中存在一个已知问题（参见问题 #31），其中 NumPy 运算符重载和 Pyomo 运算符重载最终会相互冲突（基本上，NumPy 被愚弄认为 Pyomo 标量实际上已被索引并试图将它们视为数组） .

我能够通过以下完整示例重现该错误：

# pyomo 4.4.1
from pyomo.environ import *                                                     
import numpy as np                                                              

m1 = ConcreteModel()                                                            
m1.time = Set(initialize=np.array([0,100,200,300,400,500]))                     
m1.vol_flow_in = Var(m1.time)                                                   

m1.init_value = Param(initialize = (12.3/1000)/(3.5*60))                        
m1.end_value  = Param(initialize = 0)                                           
m1.t_of_step  = Param(initialize = 210)                                         
m1.acc        = Param(initialize = 5)                                           

def _vol_flow_sigmoidal (mod,t):                                                
 return mod.vol_flow_in[t] == (mod.init_value+(mod.end_value-mod.init_value)/2)\
+((mod.end_value-mod.init_value)/2)*atan((t-mod.t_of_step)*mod.acc)/atan(1500)  
m1.vol_flow_sigmoidal = Constraint(m1.time,rule=_vol_flow_sigmoidal)

有两种可行的方法，都基于避免使用 NumPy 数组来初始化 Pyomo 集。您可以完全避免 Numpy：

m1.time = Set(initialize=[0,100,200,300,400,500])

或将 NumPy 数组显式转换为列表：

timeArray = np.array([0,100,200,300,400,500])
m1.time = Set(initialize=timeArray.tolist())

最后，为了完整起见，还有两个注意事项：

  1. 这也适用于初始化ContinuousSet对象pyomo.dae

  2. 即使您避免显式 PyomoSet声明，您也会看到相同的行为。也就是说，以下也会产生错误：

m1.time = np.array([0,100,200,300,400,500])
# ...
m1.vol_flow_sigmoidal = Constraint(m1.time,rule=_vol_flow_sigmoidal)

这是因为 Pyomo 会在幕后悄悄地为您创建 Set 对象，m1.vol_flow_sibmodial_index然后使用该 Set 来索引约束。

于 2016-09-01T03:24:56.803 回答