在混合整数规划中，如果两个二元决策变量相乘，例如 x[i,j] 和 y[j,k]，模型表达式就变成了非线性函数。很多实际业务场景只需要使用 CBC、GLPK、HiGHS 这类线性求解器，此时不处理乘积项就无法继续。其实二元变量之间的乘积完全可以精确转化成一组线性约束，关键是把乘积定义为一个新的辅助变量。
下面从原理推导到 Pyomo 代码实现，再讨论连续变量扩展和常见错误。
一、乘积项为什么会让线性求解器失效
线性规划以及混合整数线性规划要求目标函数和所有约束都是决策变量的线性函数。简单理解，每一项只能是常数乘以单个变量再相加，不允许两个变量相乘、相除或使用指数、对数等非线性函数。二元变量 x 和 y 取值虽然是零或一，但 x * y 在表达式中仍然是二次项。Pyomo 在构造表达式时会保留乘法结构，模型因此被标记为非线性。
一旦模型包含非线性项，CBC、GLPK 等 MILP 求解器无法读取。即便安装了 Ipopt、Bonmin 等非线性求解器，二元变量乘积也会逼着求解器走 MINLP 路线，求解速度通常远慢于等价线性化之后的 MILP。此外，MINLP 求解器对规模更敏感，可能难以得到全局最优解。正确做法是引入辅助变量 z 表示乘积，再添加线性不等式，使 z 的可行域恰好在乘积值为 1 时取 1，其余情况取 0。
线性化之后模型仍然保持 MILP 结构，求解器可以继续使用分支定界和割平面，性能更稳定。尤其当模型中存在大量配对变量时，例如 i、j、k 三个维度的分配关系，预先做好线性化可以显著减少调试时间。
二、二元变量乘积的等价线性化推导
设 x 和 y 都是二元变量，定义辅助变量 z = x * y。由于 x、y 只能取 0 或 1，乘积也只有 0 和 1 两种结果。我们可以用以下三个约束来刻画这个关系：
z小于等于xz小于等于yz大于等于x + y - 1
逐一验证四种组合：当 x=0、y=0 时，前两条约束使 z 小于等于 0，第三条约束使 z 大于等于 -1，因此 z=0。当 x=0、y=1 时，第一条约束使 z 小于等于 0，第三条约束使 z 大于等于 0，同样得到 0。当 x=1、y=0 时也类似。只有在 x=1、y=1 时，第三条约束变成 z 大于等于 1，而前两条只限制 z 小于等于 1，于是 z=1。因此三组不等式与乘积定义完全等价。
从模型结构看，辅助变量 z 可以被声明为连续变量并设置上下界为 0 到 1。因为约束会强制 z 取到 0 或 1，没有必要再增加一个二元变量。把 z 保持为连续变量可以减少整数变量的数量，对分支定界更友好。如果后续约束要求 z 必须为二元，也可以将 z 声明为 Binary，两种情况在最优解上是等价的。
当需要表示多个二元变量同时为 1 的乘积时，例如 z = x1 * x2 * ... * xn，可以引入边界约束 z 小于等于每个 xi，再加上 z 大于等于所有 xi 之和减去 n-1。这样 z 只有在全部 xi 为 1 时才等于 1，其他情况都为 0。不过变量数量较多时，逐对线性化再组合通常更容易维护。
三、在 Pyomo 中落地：完整模型与代码
下面用一个三层分配示例来演示。假设有若干任务 i、中间节点 j 和资源 k。二元变量 x[i,j] 表示任务 i 是否分配给节点 j，y[j,k] 表示节点 j 是否使用资源 k。当且仅当两个变量同时为 1 时，任务 i 和资源 k 之间才产生配对。为了把配对成本放入目标函数，可以定义 z[i,j,k] 表示 x[i,j] 与 y[j,k] 的乘积。
from pyomo.environ import (
ConcreteModel, Var, Constraint, Objective, Binary, NonNegativeReals,
minimize, SolverFactory, Set
)
model = ConcreteModel(name="bilinear_product_linearization")
# 三个简单集合
model.I = Set(initialize=[1, 2, 3])
model.J = Set(initialize=[1, 2, 3])
model.K = Set(initialize=[1, 2, 3])
# 原始二元变量
model.x = Var(model.I, model.J, within=Binary)
model.y = Var(model.J, model.K, within=Binary)
# 辅助变量 z[i,j,k] 表示 x[i,j] * y[j,k]
model.z = Var(model.I, model.J, model.K, within=NonNegativeReals, bounds=(0, 1))
# 每个任务只选择一个中间节点
def row_rule(m, i):
return sum(m.x[i, j] for j in m.J) == 1
model.row_cons = Constraint(model.I, rule=row_rule)
# 每个中间节点最多使用一个资源
def col_rule(m, j):
return sum(m.y[j, k] for k in m.K) <= 1
model.col_cons = Constraint(model.J, rule=col_rule)
# 线性化约束一：z 小于等于 x
def lin1(m, i, j, k):
return m.z[i, j, k] <= m.x[i, j]
model.lin1_cons = Constraint(model.I, model.J, model.K, rule=lin1)
# 线性化约束二：z 小于等于 y
def lin2(m, i, j, k):
return m.z[i, j, k] <= m.y[j, k]
model.lin2_cons = Constraint(model.I, model.J, model.K, rule=lin2)
# 线性化约束三：z 大于等于 x + y - 1
def lin3(m, i, j, k):
return m.z[i, j, k] >= m.x[i, j] + m.y[j, k] - 1
model.lin3_cons = Constraint(model.I, model.J, model.K, rule=lin3)
# 目标函数：最小化任务与资源的配对成本
cost = {
(1, 1): 2, (1, 2): 4, (1, 3): 5,
(2, 1): 3, (2, 2): 2, (2, 3): 6,
(3, 1): 4, (3, 2): 1, (3, 3): 3,
}
def obj_rule(m):
return sum(cost[i, k] * m.z[i, j, k]
for i in m.I for j in m.J for k in m.K)
model.obj = Objective(rule=obj_rule, sense=minimize)
这段代码中，z 被声明为连续变量并限制在 0 到 1 之间，这是最常用的做法。随后三条线性化约束分别对每个 i、j、k 索引建立，数量等于三层索引组合数。虽然约束数量增加，但模型仍然是 MILP，结构非常清晰。
求解时只需要将模型交给线性求解器即可。例如可以使用 CBC：
solver = SolverFactory("cbc")
solution = solver.solve(model)
for i in model.I:
for j in model.J:
for k in model.K:
if model.z[i, j, k].value > 0.5:
print(f"i={i}, j={j}, k={k}, z={model.z[i,j,k].value:.2f}")
如果本机没有安装 CBC，也可以把 "cbc" 替换为 "glpk" 或 "appsi_highs"。判断辅助变量取值时使用 0.5 作为阈值，而不是严格等于 1，能够避免求解器数值容差带来的误判。求解完成后，z 的值就对应原始乘积项，目标函数和约束都无需再出现 x[i,j] * y[j,k]。
四、扩展到二元变量与连续变量的乘积
另一种常见场景是二元变量 x 与连续变量 y 相乘，比如 x 表示某条产线是否启用，y 表示产线产量，z = x * y 表示实际产量。这种形式只需要一个 Big-M 约束组，其中 U 是 y 的上界，L 是下界。
当 x=0 时，约束 z 大于等于 L*x 且小于等于 U*x 会强制 z=0。当 x=1 时，另外两条约束 z 小于等于 y - L*(1-x) 和 z 大于等于 y - U*(1-x) 则强制 z=y。整体写出来如下：
from pyomo.environ import Var, Binary, NonNegativeReals
# y 在 0 到 100 之间连续变化
model.x = Var(within=Binary)
model.y = Var(bounds=(0, 100))
model.z = Var(within=NonNegativeReals)
U = 100
L = 0
def c1(m):
return m.z <= U * m.x
model.c1 = Constraint(rule=c1)
def c2(m):
return m.z >= L * m.x
model.c2 = Constraint(rule=c2)
def c3(m):
return m.z <= m.y - L * (1 - m.x)
model.c3 = Constraint(rule=c3)
def c4(m):
return m.z >= m.y - U * (1 - m.x)
model.c4 = Constraint(rule=c4)
这里的 U 和 L 必须紧贴 y 的真实取值范围。上界取得过大会造成 Big-M 约束松弛，分支定界时需要更多节点才能收敛；上界取得过小则可能切掉最优解。工程上应优先根据业务含义给出尽量紧的上界，而不是随意填一个很大的数。
与二元乘积不同，二元与连续变量的乘积线性化中，当 x 是二进制时这些约束是精确的，并非近似。它们经常用于固定成本、启停决策、批量生产等模型。如果 y 也有整数限制，只需把 y 声明为 Integer 或相关类型，线性化结构完全相同。
五、调试方法与常见误区
在 Pyomo 中处理乘积约束时，最常见的错误是只在数学层面知道线性化，却没有把原始乘积从模型中移除。比如目标函数里仍然写着 cost * model.x[i,j] * model.y[j,k]，同时又加上了辅助变量约束。这样模型依然是非线性，求解器一样报错。正确做法是所有需要乘积的地方都使用辅助变量 z。
如果求解时 Pyomo 提示无法将模型写入 LP 文件，或者求解器返回失败状态，可以先用 model.pprint() 检查约束。重点查看目标函数和约束中是否还残留两个变量相乘的项。另一个常见问题是辅助变量没有设置上下界，导致某些约束没有把 z 限制住，模型出现无界或不可行。对于二元乘积，设置 bounds=(0, 1) 通常会立刻解决。
数值容差也需要留意。当求解器返回 z=0.999999 或 1.000001 时，业务上应该按 0.5 的阈值判断，而不是严格比较是否等于 1。如果模型规模很大，可以考虑把辅助变量保持为连续类型，并通过减少 Big-M 值或引入 SOS1 约束来改善 LP 松弛。总的来说，二元变量乘积的线性化是一项基础但非常实用的建模技术，掌握之后可以放心地在 Pyomo 中使用线性求解器处理大量配对关系。