使用Gurobi求解复杂装箱问题时遇MLinExpr幂运算TypeError求助

尝试使用Gurobi v10.0.1求解复杂装箱问题，构建模型时遭遇TypeError，无法对变量相减得到的MLinExpr对象执行幂运算。

代码实现 #

from gurobipy import GRB, Model
import numpy as np


"""
Create data
"""
items = np.arange(50)
radii = 0.2 * np.ones_like(items)
areas = np.pi * radii ** 2                                          # Areas of items in m2
min_dist = np.add.outer(radii, radii)

pallet_length = 1.2                                                 # x-dimension of pallet in m
pallet_width = 0.8                                                  # y-dimension of pallet in m
pallet_area = pallet_length * pallet_width                          # pallet area in m2


def master(patterns, vtype):
    m = Model()
    x = m.addMVar(patterns.shape[1], lb=0, obj=1, vtype=vtype)
    contraints = m.addConstr(patterns @ x >= 1)

    m.optimize()

    if vtype == GRB.CONTINUOUS:
        return m.objVal, contraints.pi
    else:
        return m.objVal, x.X


def sub(duals):
    m = Model()
    m.setParam("NonConvex", 2)

    x = m.addMVar((len(items)), lb=0, obj=-duals, vtype=GRB.BINARY, name="x")
    cx = m.addMVar((len(items)), lb=radii, ub=pallet_length - radii, name="cx")
    cy = m.addMVar((len(items)), lb=radii, ub=pallet_width - radii, name="cy")

    # Pallet size constraint, the area of packed stacks cannot exceed the pallet area.
    m.addConstr(areas @ x <= pallet_area, name=f"area")

    for i in range(len(radii) - 1):
        for j in range(i + 1, len(radii)):
            lhs = (cx[i] - cx[j])** 2 + (cy[i] - cy[j]) ** 2
            rhs = (x[i] + x[j] - 1) * min_dist[i, j] ** 2
            m.addConstr(lhs >= rhs, name=f"overlap[{i},{j}]")

    def callback(model, where):
        if where == GRB.Callback.MIPNODE:
            if model.cbGet(GRB.Callback.MIPNODE_OBJBST) < -1:
                model.terminate()

    m.optimize(callback)

    return x.X if m.objVal < -1 else None


# Initial pattern: put every item on its own pallet.
patterns = np.eye(len(items))
num_pallets, duals = master(patterns, GRB.CONTINUOUS)

while (new_pattern := sub(duals)) is not None:
    print("Current number of pallets needed (LP):", num_pallets)

    # Check if we already have this pattern. If so, we are done as well.
    if any(np.allclose(new_pattern, column) for column in patterns.T):
        break

    patterns = np.column_stack((patterns, new_pattern))
    num_pallets, duals = master(patterns, GRB.CONTINUOUS)

min_pallets, used_patterns = master(patterns, GRB.BINARY)
unused = set(items)

print("Number of pallets needed (IP):", min_pallets)

for idx, used in enumerate(used_patterns):
    if used:
        on_pallet = items[patterns[:, idx].astype(bool)]
        on_pallet = {item for item in on_pallet if item in unused}
        unused -= on_pallet
        print(f"Put items {on_pallet} on same pallet.")

错误信息 #

lhs = (cx[i] - cx[j])** 2 + (cy[i] - cy[j]) ** 2
          ~~~~~~~~~~~~~~^^~~~
TypeError: unsupported operand type(s) for ** or pow(): 'MLinExpr' and 'int'

解决方案 #

Gurobi支持二次约束，但不能直接对MLinExpr对象使用**2运算符，需要通过变量相乘的方式构建二次表达式，以下是两种可行的修改方式：

方式1：展开二次项 #

将平方项手动展开，直接构建线性与二次项组合的约束：

for i in range(len(radii) - 1):
    for j in range(i + 1, len(radii)):
        # 展开(xi - xj)²和(yi - yj)²
        lhs_x = cx[i] * cx[i] - 2 * cx[i] * cx[j] + cx[j] * cx[j]
        lhs_y = cy[i] * cy[i] - 2 * cy[i] * cy[j] + cy[j] * cy[j]
        lhs = lhs_x + lhs_y
        rhs = (x[i] + x[j] - 1) * min_dist[i, j] ** 2
        m.addConstr(lhs >= rhs, name=f"overlap[{i},{j}]")

方式2：使用QuadExpr构建二次表达式 #

借助Gurobi的QuadExpr对象组装二次项，代码可读性更强：

from gurobipy import QuadExpr  # 需导入QuadExpr

for i in range(len(radii) - 1):
    for j in range(i + 1, len(radii)):
        lhs = QuadExpr()
        # 添加x方向的二次项
        lhs.add(cx[i], cx[i])
        lhs.add(cx[j], cx[j])
        lhs.add(-2, cx[i], cx[j])
        # 添加y方向的二次项
        lhs.add(cy[i], cy[i])
        lhs.add(cy[j], cy[j])
        lhs.add(-2, cy[i], cy[j])
        rhs = (x[i] + x[j] - 1) * min_dist[i, j] ** 2
        m.addConstr(lhs >= rhs, name=f"overlap[{i},{j}]")

额外注意事项 #

• 你已经设置m.setParam("NonConvex", 2)，这是正确的，因为这类二次约束属于非凸问题，必须开启Gurobi的非凸求解模式才能处理。
• 你的Gurobi版本v10.0.1完全支持非凸二次规划，无需版本升级。

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