Pyomo中GDP半连续变量混合问题报错及解决咨询

刚接触Pyomo的GDP建模，尝试求解书中第163页“9.4 含半连续变量的混合问题”示例时出现报错。示例定义了三个半连续变量x₁、x₂、x₃，目标为最小化非零变量数量，同时满足体积约束，但运行时抛出PyomoException，提示无法将非恒定表达式转换为布尔值。

完整代码 #

# scont.py
from pyomo.environ import*
from pyomo.gdp import*

L = [1,2,3]
U = [2,4,6]
index = [0,1,2]

model = ConcreteModel()
model.x = Var(index, within=Reals, bounds=(0,20))

# Each disjunct is a semi-continuous variable
# x[k] == 0 or L[k] <= x[k] <= U[k]


def d_rule(block, k, i):
    m = block.model()
    if i==0:
        block.c = Constraint(expr=m.x[k] == 0)
    else:
        block.c = Constraint(expr=L[k] <= m.x[k] <= U[k])
    
model.d = Disjunct(index, [0,1], rule=d_rule)

#There are three disjunctions

def D_rule(block, k):
    model = block.model()
    return[model.d[k,0], model.d[k,1]]

model.D = Disjunction(index, rule=D_rule)

# Minimize the number of x variables that are non zero 

model.o = Objective(expr=sum(model.d[k,1].indicator_var for k in index))


# Satisfy a demand that is met by these variables

model.c = Constraint(expr=sum(model.x[k]for k in index)>= 7)


xfrm = TransformationFactory('gdp.bigm')
xfrm.apply_to(model)

solver = SolverFactory('glpk')
status = solver.solve(model)

报错信息 #

ERROR: Constructing component 'd' from data=None failed: PyomoException:
Cannot convert non-constant Pyomo expression (1  <=  x[0]) to bool. This
error is usually caused by using a Var, unit, or mutable Param in a
Boolean context such as an "if" statement, or when checking container
membership or equality. For example,
m.x = Var() >>> if m.x >= 1: ...     pass
and
m.y = Var() >>> if m.y in [m.x, m.y]: ...     pass
would both cause this exception.

报错核心原因是：Pyomo不支持在约束中直接使用L[k] <= m.x[k] <= U[k]这种链式比较写法——它会尝试将整个表达式转换为布尔值，但表达式包含变量，无法直接完成转换。需要将链式比较拆分为两个独立的约束。

修改后的完整代码如下：

# scont.py
from pyomo.environ import*
from pyomo.gdp import*

L = [1,2,3]
U = [2,4,6]
index = [0,1,2]

model = ConcreteModel()
model.x = Var(index, within=Reals, bounds=(0,20))

# 每个析取式对应半连续变量的两种状态：0 或 处于[L[k], U[k]]区间
def d_rule(block, k, i):
    m = block.model()
    if i == 0:
        # x[k]等于0的情况
        block.c_eq = Constraint(expr=m.x[k] == 0)
    else:
        # 将链式约束拆分为上下界两个独立约束
        block.c_lb = Constraint(expr=m.x[k] >= L[k])
        block.c_ub = Constraint(expr=m.x[k] <= U[k])
    
model.d = Disjunct(index, [0,1], rule=d_rule)

# 定义析取关系：每个x[k]必须满足两种状态之一
def D_rule(block, k):
    return [model.d[k,0], model.d[k,1]]

model.D = Disjunction(index, rule=D_rule)

# 目标：最小化非零变量的数量（显式指定最小化方向）
model.o = Objective(expr=sum(model.d[k,1].indicator_var for k in index), sense=minimize)

# 体积约束：变量总和不小于7
model.c = Constraint(expr=sum(model.x[k] for k in index) >= 7)

# 应用GDP转MIP的big-M转换
xfrm = TransformationFactory('gdp.bigm')
xfrm.apply_to(model)

# 调用求解器并输出结果
solver = SolverFactory('glpk')
status = solver.solve(model)

print("求解状态：", status)
print("非零变量数量：", value(model.o))
print("变量取值：")
for k in index:
    print(f"x[{k}] = {round(value(model.x[k]), 2)}")

关键修改说明 #

1. 拆分链式约束：将L[k] <= m.x[k] <= U[k]拆分为m.x[k] >= L[k]和m.x[k] <= U[k]两个独立约束，避免Pyomo尝试将表达式转换为布尔值。
2. 显式指定目标方向：虽然Pyomo默认目标为最小化，但显式写出sense=minimize让代码逻辑更清晰。
3. 添加结果输出：增加打印代码，方便直接查看求解后的变量取值和目标值。

内容的提问来源于stack exchange，提问作者Samuel M