Pyomo调用Gurobi求解含阶梯分段函数的二次规划未得到最优解

• 你使用Pyomo的Piecewise组件搭配SOS2表示法时，组件默认对相邻定义域断点做线性插值，而非你预期的阶梯常数映射：在你给出的3个定义域点[5,10,11]之间，Z会被模型允许取10~20.5、9~10之间的插值，此时X*Z的最大值会出现在插值区间内（你得到的7.37*15.53≈114远大于你预期的5*20.5=102.5），因此结果不符合需求。
• 当前的断点配置没有显式定义阶梯的常数段，无法强制Z只能取20.5/10/9三个固定值。

方案1：调整Piecewise断点适配阶梯逻辑 #

将每个阶梯区间的左右端点都显式加入定义域列表，对应区间内Z值相同，此时线性插值结果为常数，符合阶梯要求：

from pyomo.core import *

# 调整断点：每个阶梯区间左右端点都列示，对应Z值不变
# 若需要支持X<5的场景，可修改为DOMAIN_PTS = [0.,5.,5.,10.,10.,11.]，RANGE_PTS = [20.5,20.5,10.,10.,9.,9.]，同时调整X的下界为0
DOMAIN_PTS = [5., 5., 10., 10., 11.]
RANGE_PTS = [20.5, 10., 10., 9., 9.]

# Define model and variables
model = ConcreteModel()
model.X = Var(bounds=(5,11))
model.Z = Var()

# Set piecewise constraint
model.con = Piecewise(model.Z,model.X, 
                      pw_pts=DOMAIN_PTS ,
                      pw_constr_type='EQ',
                      f_rule=RANGE_PTS,
                      force_pw=True,
                      pw_repn='SOS2') 

model.obj = Objective(expr= model.Z * model.X, sense=maximize)
opt = SolverFactory('gurobi')
opt.options['NonConvex'] = 2
opt.solve(model)

print(value(model.X)) # 输出5.0
print(value(model.Z)) # 输出20.5
print(value(model.obj)) # 输出102.5

方案2：手动写大M约束实现（更灵活，支持严格区间边界） #

如果需要实现区间不重叠的严格边界（比如5 < X <=10），可以用二进制变量+大M约束直接定义阶梯逻辑，可控性更强：

from pyomo.core import *

model = ConcreteModel()
model.X = Var(bounds=(0,11)) # 可根据需求调整X的上下界
model.Z = Var()
# 定义三个二进制变量分别对应三个区间
model.y1 = Var(domain=Binary) # X<=5的标识
model.y2 = Var(domain=Binary) # 5<X<=10的标识
model.y3 = Var(domain=Binary) # 10<X<=11的标识
M = 1e5 # 大M值取大于变量上下界差的数值即可

# 约束：只能有一个区间激活
model.sum_y = Constraint(expr = model.y1 + model.y2 + model.y3 == 1)
# 区间1约束：y1=1时 X<=5，Z=20.5
model.x1_con = Constraint(expr = model.X <= 5 + M*(1-model.y1))
model.z1_con = Constraint(expr = model.Z == 20.5 - M*(1-model.y1))
# 区间2约束：y2=1时 5<X<=10，Z=10
model.x2_con1 = Constraint(expr = model.X >=5 + 1e-6 - M*(1-model.y2)) # 1e-6处理严格大于边界，不需要可去掉
model.x2_con2 = Constraint(expr = model.X <=10 + M*(1-model.y2))
model.z2_con = Constraint(expr = model.Z == 10 - M*(1-model.y2))
# 区间3约束：y3=1时 10<X<=11，Z=9
model.x3_con1 = Constraint(expr = model.X >=10 + 1e-6 - M*(1-model.y3))
model.x3_con2 = Constraint(expr = model.X <=11 + M*(1-model.y3))
model.z3_con = Constraint(expr = model.Z ==9 - M*(1-model.y3))

model.obj = Objective(expr= model.Z * model.X, sense=maximize)
opt = SolverFactory('gurobi')
opt.options['NonConvex'] = 2
opt.solve(model)

print(value(model.X))
print(value(model.Z))
print(value(model.obj))

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