Docplex自定义回调实现子问题求解与割平面添加问题求助

我现在正尝试实现一个自定义回调函数，在主模型的优化过程中求解一个子问题，具体想完成这几个步骤：

• 从主模型的当前解中获取特定变量的取值
• 基于这些值求解一个独立的子优化问题
• 如果子问题的解提示需要添加割平面，就把这个割平面作为新约束加入主模型

我想用Python实现，但现在完全摸不着头绪，也试过Xpress之类的工具，但文档对我帮助不大。下面是我写的代码，真心希望能得到大家的帮助！

from docplex.mp.model import Model
from cplex.callbacks import UserCutCallback
from docplex.mp.callbacks.cb_mixin import ConstraintCallbackMixin
import random

class CustomCutCallback(ConstraintCallbackMixin, UserCutCallback):
    def __init__(self, env):
        UserCutCallback.__init__(self, env)
        ConstraintCallbackMixin.__init__(self)
        self.eps = 1e-6
        self.nb_cuts = 0
        self.cts = []

    def add_cut_constraint(self, cut):
        self.register_constraint(cut)

    def __call__(self, context):
        """
        This method is invoked by CPLEX during optimization.
        It receives the 'context' to interact with the model and variables.
        """
        print("Callback called with context")
        m = context.model  # The main model from the context
        m_sol = m.solution

        x_values = {var: m_sol.get_value(var) for var in m.variables if var.name.startswith('x')}
        w_values = {var: m_sol.get_value(var) for var in m.variables if var.name.startswith('w')}
        m2_cuts = self.solve_subproblem(x_values, w_values, m, context)

    def solve_subproblem(self, x_values, w_values, m, context):
        """
        Solves the subproblem and generates cuts if necessary.
        """
        m2 = Model(name='subproblem')

        client_range = range(len(x_values))
        deposit_range = range(len(w_values))
        plant_range = range(len(w_values[0]))

        alpha = m2.continuous_var_matrix(client_range, deposit_range, name='alpha', lb=0)
        beta = m2.continuous_var_matrix(deposit_range, plant_range, name='beta', lb=0)

        m2.add_constraints(alpha[i, j] + (x_values.get((i, j), 0) * beta[j, k]) <= x_values.get((i, j), 0) * w_values.get((j, k), 0)
                           for i in client_range for j in deposit_range for k in plant_range)

        m2.maximize(m2.sum(alpha[i, j] * x_values.get((i, j), 0) for i in client_range for j in deposit_range) + 
                    m2.sum(beta[j, k] * w_values.get((j, k), 0) for j in deposit_range for k in plant_range))

        m2.solve()
        print(m2.solution)

        # Here, you perform an evaluation of the cut values
        for i in client_range:
            for j in deposit_range:
                for k in plant_range:
                    if  sum(m2.solution.get_value(alpha[i, j]) * x_values.get((i, j), 0) for i in client_range for j in deposit_range) + \
                    sum(m2.solution.get_value(beta[j, k]) * w_values.get((j, k), 0) for j in deposit_range for k in plant_range) > \
                    sum(transport_cost_deposit_client[j][i] * d[i] * x[i, j] for i in client_range for j in deposit_range) + \
                    sum(transport_cost_plant_deposit[j][k] * w[j, k] for j in deposit_range for k in plant_range):
                        m.add_constraint(sum(m2.solution.get_value(alpha[i, j]) * x[i, j] for i in client_range for j in deposit_range) + \
                                         sum(m2.solution.get_value(beta[j, k]) * w[j, k] for j in deposit_range for k in plant_range))


# Main model
def build_location_model(transport_cost_deposit_client, transport_cost_plant_deposit, p, q, d, **kwargs):
    m = Model(name='location', **kwargs)

    num_deposits = len(transport_cost_deposit_client)
    num_plants = len(transport_cost_plant_deposit[0])
    num_clients = len(d)
    
    deposit_range = range(num_deposits)
    plant_range = range(num_plants)
    client_range = range(num_clients)

    x = m.binary_var_matrix(client_range, deposit_range, name='x')
    w = m.integer_var_matrix(deposit_range, plant_range, name='w', lb=0)
    y = m.binary_var_list(deposit_range, name='y')
    h = m.binary_var_list(plant_range, name='h')

    m.add_constraints(m.sum(x[i, j] for j in deposit_range) == 1 for i in client_range)
    m.add_constraints(m.sum(w[j, k] for k in plant_range) == y[j] for j in deposit_range)
    m.add_constraints(x[i, j] <= y[j] for i in client_range for j in deposit_range)
    m.add_constraint(m.sum(h[k] for k in plant_range) == q)
    m.add_constraint(m.sum(y[j] for j in deposit_range) == p)
    m.add_constraints(m.sum(d[i] * x[i,j] for i in client_range) == m.sum(w[j, k] for k in plant_range) for j in deposit_range)
    m.add_constraints(w[j, k] <= (m.sum(d[i] for i in client_range) * h[k]) for j in deposit_range for k in plant_range)

    transport_cost = m.sum(transport_cost_deposit_client[j][i] * d[i] * x[i, j] for i in client_range for j in deposit_range) + \
                        m.sum(transport_cost_plant_deposit[j][k] * w[j, k] for j in deposit_range for k in plant_range)
    m.minimize(transport_cost)
    m.parameters.preprocessing.presolve = 0

    # Register the callback
    cut_cb = m.register_callback(CustomCutCallback)

    # Configure CPLEX parameters for cuts
    params = m.parameters
    params.mip.cuts.mircut = -1

    m.solve()

    return m



# Test function
def solve_model():
    num_deposits = 10
    num_plants = 4
    num_clients = 20

    TRANSPORT_COST_DEPOSITS_CLIENTS = [
        [random.randint(20, 100) for _ in range(num_clients)] for _ in range(num_deposits)
    ]

    TRANSPORT_COST_PLANTS_DEPOSITS = [
        [random.randint(30, 80) for _ in range(num_plants)] for _ in range(num_deposits)
    ]

    p = 5
    q = 3
    d = [random.randint(5, 20) for _ in range(num_clients)]

    m = build_location_model(TRANSPORT_COST_DEPOSITS_CLIENTS, TRANSPORT_COST_PLANTS_DEPOSITS, p, q, d)
    if m.solution is None:
        print("No valid solution found.")
        return None
    print(m.solution)
    return m


if __name__ == "__main__":
    solve_model()

针对你代码的问题分析与修复建议 #

1. 回调中获取变量值的正确方式 #

你当前用m.solution.get_value(var)获取的是全局最优解（如果存在的话），但回调触发时我们需要的是当前MIP节点的松弛解，应该用context.get_values(var)来获取：

# 更高效的方式：提前把主模型的x、w变量传给回调，避免遍历所有变量
x_vals = {(i,j): context.get_values(self.x_vars[i,j]) for i,j in self.x_vars}

2. 子问题维度匹配错误 #

你用len(x_values)来推导client_range是错误的——x_values是变量字典，长度是x变量的总数，不是客户数量。建议在回调初始化时，直接传入主模型的client_range、deposit_range、plant_range以及x、w变量矩阵，避免动态推导出错。

3. 割平面添加的错误逻辑 #

在回调中不能直接调用m.add_constraint()修改主模型，必须用ConstraintCallbackMixin提供的register_constraint方法。同时你代码中transport_cost_deposit_client等参数在回调中未定义，需要提前传入回调。

4. 子问题性能优化 #

每次回调都重建子模型会非常低效，建议提前初始化子模型结构，每次仅更新约束系数和目标函数即可。

调整后的核心回调代码示例 #

class CustomCutCallback(ConstraintCallbackMixin, UserCutCallback):
    def __init__(self, env, main_model, x_vars, w_vars, client_range, deposit_range, plant_range, tc_dc, tc_pd, d):
        UserCutCallback.__init__(self, env)
        ConstraintCallbackMixin.__init__(self)
        self.eps = 1e-6
        self.nb_cuts = 0
        self.main_model = main_model
        self.x_vars = x_vars  # 主模型x[i,j]变量矩阵
        self.w_vars = w_vars  # 主模型w[j,k]变量矩阵
        self.client_range = client_range
        self.deposit_range = deposit_range
        self.plant_range = plant_range
        self.tc_dc = tc_dc
        self.tc_pd = tc_pd
        self.d = d
        # 提前初始化子问题
        self._init_subproblem()

    def _init_subproblem(self):
        self.sub_model = Model(name='subproblem')
        self.alpha = self.sub_model.continuous_var_matrix(self.client_range, self.deposit_range, name='alpha', lb=0)
        self.beta = self.sub_model.continuous_var_matrix(self.deposit_range, self.plant_range, name='beta', lb=0)
        # 初始化空约束，后续更新系数
        self.sub_constraints = []
        for i in self.client_range:
            for j in self.deposit_range:
                for k in self.plant_range:
                    self.sub_constraints.append(
                        self.sub_model.add_constraint(self.alpha[i,j] + 0 <= 0)
                    )

    def _update_subproblem(self, x_vals, w_vals):
        # 更新子问题约束与目标
        idx = 0
        for i in self.client_range:
            for j in self.deposit_range:
                for k in self.plant_range:
                    x_val = x_vals[i,j]
                    w_val = w_vals[j,k]
                    # 更新约束
                    self.sub_constraints[idx].lhs = self.alpha[i,j] + x_val * self.beta[j,k]
                    self.sub_constraints[idx].rhs = x_val * w_val
                    idx += 1
        # 更新目标函数
        obj_expr = self.sub_model.sum(self.alpha[i,j] * x_vals[i,j] for i,j in x_vals) + \
                   self.sub_model.sum(self.beta[j,k] * w_vals[j,k] for j,k in w_vals)
        self.sub_model.set_objective('max', obj_expr)

    def __call__(self, context):
        # 仅在MIP节点松弛解可用时触发
        if context.in_relaxation():
            print(f"回调触发于节点 {context.get_node_id()}")
            # 获取当前节点的变量值
            x_vals = {(i,j): context.get_values(self.x_vars[i,j]) for i in self.client_range for j in self.deposit_range}
            w_vals = {(j,k): context.get_values(self.w_vars[j,k]) for j in self.deposit_range for k in self.plant_range}
            
            # 更新并求解子问题
            self._update_subproblem(x_vals, w_vals)
            sub_sol = self.sub_model.solve()
            if not sub_sol:
                return
            
            sub_obj = sub_sol.objective_value
            # 计算主模型对应表达式的当前解值
            main_expr_val = sum(self.tc_dc[j][i] * self.d[i] * x_vals[i,j] for i,j in x_vals) + \
                           sum(self.tc_pd[j][k] * w_vals[j,k] for j,k in w_vals)
            
            # 判断是否需要添加割平面
            if sub_obj > main_expr_val + self.eps:
                # 构造割平面约束（使用主模型变量）
                cut_expr = self.main_model.sum(
                    sub_sol.get_value(self.alpha[i,j]) * self.x_vars[i,j] for i,j in x_vals
                ) + self.main_model.sum(
                    sub_sol.get_value(self.beta[j,k]) * self.w_vars[j,k] for j,k in w_vals
                ) <= sub_obj
                self.register_constraint(cut_expr)
                self.nb_cuts += 1
                print(f"已添加第 {self.nb_cuts} 个割平面")

注册回调的修改 #

在主模型build_location_model函数中，注册回调时传入必要参数：

# 注册回调
cut_cb = m.register_callback(CustomCutCallback, m, x, w, client_range, deposit_range, plant_range, transport_cost_deposit_client, transport_cost_plant_deposit, d)

备注：内容来源于stack exchange，提问作者Javier