Pyomo抽象建模中的嵌套Disjunctions

Question

我正在研究一个带有一些析取的小型优化模型。我在具体模型中的做法效果很好：

from pyomo.environ import *
m = ConcreteModel()
m.d1 = Disjunct()
m.d2 = Disjunct()
m.d1.sub1 = Disjunct()
m.d1.sub2 = Disjunct()
m.d1.disj = Disjunction(expr=[m.d1.sub1, m.d1.sub2])
m.disj = Disjunction(expr=[m.d1, m.d2])

但现在我把具体的模型转换成了抽象的公式。我能够修复所有内容，而不是嵌套析取。我的做法是：

#Disjunct 1        
def _op_mode1(self, op_mode, t):
m = op_mode.model()
op_mode.c1 = po.Constraint(expr=m.x[t] == True)
#Disjunct 2       
def _op_mode2(self, op_mode, t):
m = op_mode.model()
op_mode.c1 = po.Constraint(expr=m.x[t] == False)
#Disjunction 1
def _op_modes(self,m, t):
return [m.mode1[t], m.mode2[t]]
#Adding Components
self.model.del_component("mode1")
self.model.del_component("mode1_index")
self.model.add_component("mode1", pogdp.Disjunct(self.model.T, rule=self._op_mode1))
self.model.del_component("mode2")
self.model.del_component("mode2_index")
self.model.add_component("mode2", pogdp.Disjunct(self.model.T, rule=self._op_mode1))
self.model.del_component("modes")
self.model.del_component("modes_index")
self.model.add_component("modes", pogdp.Disjunction(self.model.T, rule=self._op_modes))`

正如我之前提到的，这很好用。但我还没有找到任何方法来嵌套这些脱节。Pyomo alsways抱怨析取的第二层；sub1"；。

有人能给我一个提示吗？

许多问候

Joerg

Answer 1

上述最新模型的问题是，您为m.T的每个元素声明m.d1和m.d2，但它们每次都会相互覆盖，因为它们具有相同的名称。您应该看到为此记录的警告消息。因此，如果取消对模型的pprint的注释，您将看到只有最后一个声明(对x[10]有约束(。因此，m.disjunction_中的前9个不连接是不存在的不连接的析取。最简单的修复方法是在声明析取时为它们提供唯一的名称：

import pyomo.environ as pyo
import pyomo.gdp as pogdp
model = pyo.ConcreteModel()
model.T = pyo.RangeSet(0, 10)
model.x=pyo.Var(model.T,bounds=(-2, 10))
model.y=pyo.Var(model.T,bounds=(20, 30))
# This was also a duplicate declaration:
#model.disjunction_ = pogdp.Disjunction(model.T)
def d1(m, t):
disj = pogdp.Disjunct()
disj.c1= pyo.Constraint(expr=m.x[t] <= 10)
m.add_component('d1_%s' % t, disj)
return disj
def d2(m, t):
disj = pogdp.Disjunct()
disj.c1= pyo.Constraint(expr=m.x[t] >= 10)
m.add_component('d2_%s' % t, disj)
return disj
# sum x,y
def obj_rule(m):
return pyo.quicksum(pyo.quicksum([m.x[t] + m.y[t]], linear=False) for t in
m.T)
model.obj = pyo.Objective(rule=obj_rule)
def _op_mode_test(m, t):
disj1 = d1(m, t)
disj2 = d2(m, t)
return [disj1, disj2]
model.disjunction_ = pogdp.Disjunction(model.T, rule=_op_mode_test)

然而，通过m.T对Disjuncts进行索引也会更干净(而且可能更容易(，因为这基本上就是唯一名称所做的。

Answer 2

Block(以及Disjunct规则(被传递给要填充的块(或析取(作为第一个参数。因此，一个"；摘要"；你的具体模型可能看起来像这样：

model = AbstractModel()
@model.Disjunct()
def d1(d):
# populate the `d` disjunct (i.e., `model.d1`) here
pass
@model.Disjunct()
def d2(d):
@d.Disjunct()
def sub1(sd):
# populate the 'sub1' disjunct here
pass
@d.Disjunct()
def sub2(sd):
# populate the 'sub2' disjunct here
pass
d.disj = Disjunction(expr=[d.sub1, d.sub2])
model.disj = Disjunction(expr=[model.d1, model.d2])

还有一个更基本的问题是为什么您要将您的模型转换为"；摘要"；Pyomo抽象模型大多是为来自AMPL建模的人所熟悉而设计的。虽然它们将与块结构模型一起工作，因为AMPL从未真正考虑到块，类似的面向块的抽象模型往往不必要地繁琐。

Answer 3

这是我们的新模型：

import pyomo.environ as pyo
import pyomo.gdp as pogdp
model = pyo.ConcreteModel()
model.T = pyo.RangeSet(0,10)
model.x=pyo.Var(model.T,bounds=(-2, 10))
model.y=pyo.Var(model.T,bounds=(20, 30))
model.disjunction_=pogdp.Disjunction(model.T)
def d1(m,t):
m.d1 = pogdp.Disjunct()
m.d1.c1= pyo.Constraint(expr=m.x[t] <=10)
def d2(m,t):
m.d2 = pogdp.Disjunct()
m.d2.c1= pyo.Constraint(expr=m.x[t] >=10)
# sum x,y
def obj_rule(m):
return pyo.quicksum(pyo.quicksum([m.x[t] + m.y[t]], linear=False) for t in m.T)
model.obj = pyo.Objective(rule=obj_rule)
def _op_mode_test(m,t):
d1(m,t)
d2(m,t)
return [m.d1,m.d2]
model.disjunction_=pogdp.Disjunction(model.T,rule=_op_mode_test)
#model.pprint()
pyo.TransformationFactory('gdp.bigm').apply_to(model)
solver = pyo.SolverFactory('baron')
solver.solve(model)
print(pyo.value(model.obj))

我认为这与RangeSet有关。对于单个步骤，它是有效的，但对于多个步骤，它会抛出错误：AttributeError:"NoneType"对象没有属性"component">

如果你能看一下就太好了。

非常感谢

Answer 4

以下是适用于bigm但不适用于mbigm或外壳转换的代码：

import pyomo.environ as pyo
import pyomo.gdp as pogdp
model = pyo.ConcreteModel()
model.T = pyo.RangeSet(2)
model.x=pyo.Var(model.T,bounds=(1, 10))
model.y=pyo.Var(model.T,bounds=(1, 100))

def _op_mode_sub(m, t):
m.disj1[t].sub1 = pogdp.Disjunct()
m.disj1[t].sub1.c1= pyo.Constraint(expr=m.y[t] == 60)
m.disj1[t].sub2 = pogdp.Disjunct()
m.disj1[t].sub2.c1= pyo.Constraint(expr=m.y[t] == 100)
return [m.disj1[t].sub1, m.disj1[t].sub2]
def _op_mode(m, t):
m.disj2[t].c1= pyo.Constraint(expr=m.y[t] >= 3)
m.disj2[t].c2= pyo.Constraint(expr=m.y[t] <= 5)
return [m.disj1[t], m.disj2[t]]
model.disj1 = pogdp.Disjunct(model.T)
model.disj2 = pogdp.Disjunct(model.T)
model.disjunction1sub = pogdp.Disjunction(model.T, rule=_op_mode_sub)
model.disjunction1 = pogdp.Disjunction(model.T, rule=_op_mode)

def obj_rule(m, t):
return pyo.quicksum(pyo.quicksum([m.x[t] + m.y[t]], linear=False) for t in m.T)
model.obj = pyo.Objective(rule=obj_rule)
model.pprint()
gdp_relax=pyo.TransformationFactory('gdp.bigm')
gdp_relax.apply_to(model)
solver = pyo.SolverFactory('glpk')
solver.solve(model)
print(pyo.value(model.obj))