PYOMO抽象建模中的嵌套分离

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.d1和m.d2 mt的每个元素，但它们覆盖了每个元素每次都有相同的名字以来。您应该看到为此记录的警告消息。因此，如果您将模型的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)

但是，通过mt索引划分的清洁度（可能更容易），因为这基本上是基本上是独特的名称在做什么。

Answer 2

block（因此，ed规则）被传递给块（或分离）作为第一个参数。因此，您的具体模型也可能看起来像这样：

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：'nontype'对象没有属性'component'，

如果您可以看一下它，那将是很棒的。

非常感谢

Answer 4

这是与BIGM相当良好的代码，但与MBIGM或Hull Transformation相当不错：

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))