关于Bender decomposition分解算法解析附python代码

1.Bender decomposition分解算法的原理

2.Bender decomposition算法的python实现

最近在学习鲁棒优化的时候看到了一个新算法也就是本文的Bender 分解算法 本着求知的态度，看了网上的众多解释可教程，大概算是有点了解了，顺手写个blog记录一下吧，附上复现的python代码和matlab代码

一、Bender decomposition算法的原理

Bender分解算法是用来求解MIP问题的，特别是对于大规模混合整数规划问题有奇效，注：当然这也是我道听途说，gpt也是这样说的，实际上的求解效率笔者没有和其他精确算法进行对比过
其本质思想就是将复杂的变量给置换出来，考虑如下模型(MIP)：
min ⁡ x , y c x + f y s . t . A x + B y = b x ⩾ 0 y ∈ Y ⊆ R n \begin{aligned}&\min_{x,y}&&cx+fy\\&\mathrm{s.t.}&&Ax+By=b\\&&&x\geqslant0\\&&&y\in Y\subseteq\mathbb{R}^n\end{aligned} ​x,ymin​s.t.​​cx+fyAx+By=bx⩾0y∈Y⊆Rn​
其中x是连续决策变量，y是整数变量
其中 c ∈ R 1 × m , f ∈ R 1 × n c\in\mathbb{R}^{1\times m},f\in\mathbb{R}^{1\times n} c∈R1×m,f∈R1×n ,是行向量，其维度分别为 m × 1 m\times1 m×1和 n × 1 n\times1 n×1。 A , B A,B A,B是约束矩阵，其维度分别为 A ∈ R r × m , B ∈ R r × n A\in\mathbb{R}^{r\times m},B\in\mathbb{R}^{r\times n} A∈Rr×m,B∈Rr×n 。 b b b为右端常数项。
这个模型进行求解的时候是比较困难的，于是Jacques F. Benders 在1962年想到将y这个复杂变量给提前确定好，也就是固定一个值，于是就得了一个LP问题：
SP: min ⁡ x c x s.t. A x = b − B y ˉ x ⩾ 0 \begin{aligned}&\text{SP:}\quad\min_{x}\quad cx\\&\text{s.t.}\quad Ax=b-B\bar{y}\\&x\geqslant0\end{aligned} ​SP:xmin​cxs.t.Ax=b−Byˉ​x⩾0​
这个问题里面只有一个线性变量 x x x不是np-hard问题，自然可以很轻松的求解，学术上也就称这个问题为子问题(subproblem)
同时对比初始问题余下的也就叫主问题(master problem),如下所示：
M P : min ⁡ y f y s . t . y ∈ Y ⊆ R n \begin{array}{rcl}\mathrm{MP}:&\min_y&fy\\&\mathrm{s.t.}&y\in Y\subseteq\mathbb{R}^n\end{array} MP:​miny​s.t.​fyy∈Y⊆Rn​
这里可以看出sp问题是无界的时候，也就是存在一个y的值可以让子问题的最优解是无界的，所以显然原问题也是无界的
显然求出来两个问题的最优值之和并不是原问题的最优解，主问题显然没有考虑到子问题的约束而直接求解，子问题是在主问题得出的任意一个解上做出的优化，显然没有考虑到两个约束问题之间的关联。所以这里涉及到一个专业名词recourse，就是以一个切平面的形式追加到主问题中。这两个切平面是: ( α p i ) T ( b − B y ) ⩽ q , ∀ i ∈ P (\alpha_p^i)^T(b-By)\leqslant q,\quad\forall i\in\mathcal{P} (αpi​)T(b−By)⩽q,∀i∈P ( α p i ) T ( b − B y ) ⩽ 0 , ∀ i ∈ P (\alpha_p^i)^T(b-By)\leqslant 0,\quad\forall i\in\mathcal{P} (αpi​)T(b−By)⩽0,∀i∈P
前者是叫做最优切平面，后者是可行切平面 （这里涉及到了极点和极线的知识，不懂的可以看看关于这个分解算法的经典文献或者看看晚上的一些讲解，不求甚解）来源于子问题的对偶问题(看sp问题可以看到其duality问题的约束条件里面并没有复杂变量y）。这里的这个p是一个辅助变量， α p i \alpha_p^i αpi​在前者是极点，后者是极线
所以这个时候的主问题mp就变成了：
M P : min ⁡ y f y + q \mathrm M\mathrm P:\quad\min_y\quad fy+q MP:ymin​fy+q
s . t . y ∈ Y ⊆ R n s.t. y\in Y\subseteq \mathbb{R} ^n s.t.y∈Y⊆Rn

B e n d e r s   o p t i m a l i t y   c u t s Benders\ optimality \ cuts Benders optimality cuts

B e n d e r s   f e a s i b i l i t y   c u t s Benders \ feasibility \ cuts Benders feasibility cuts

q f r e e qfree qfree
这里的q就是一个替代cx的一个变量，也就是添加子问题影响的变量，令子问题的解是 f ( x ) f(x) f(x),这个q是他的一个下界,所以主要的步骤如下所示：
1.求解主问题的解，给出一个y来，
2.将y给入到子问题中，得到子问题的最优解(两种方法得到这个问题的对偶问题)，一种写出对偶问题来，另一种是不需要写出，来利用所用求解器中特定的函数来进行获取，这个具体可以查看官方文档。首先第一种方式中 当对偶问题是有界可行解的时候 加入第一个切平面，是无界的时候则加入二个切平面，对于第二种方式，当子问题是有界可行的时候使用第一个切平面，是无可行解的时候使用第二个切平面，切平面加入到MP问题中。这里的这个技巧好像和原问题和对偶问题之间的对应关系差不多，具体可见：运小筹有关Bender分解算法的讲解，这个才是大佬
3.更新这个问题的上界和下界，同时循环，直到上界和下届相等，这里的上界显然是 f y + f ( x ) fy+f(x) fy+f(x),因为这里的 f ( x ) f(x) f(x)是任意从主问题中给定一个y求出来的界，所以显然是一个可行解，由于是min问题，所以是MIP问题的一个上界，下界是 f y + p fy+p fy+p这里是因为这里的切平面每次只是添加一部分切平面所以可以看作原问题的一个松弛问题，自认可行域要比原问题MIP的范围要广，也就比其最优值要小。
4.对mp问题进行求解这里的变量有两个也就是y和q，所以可以将p进行更新，也就可以得到新的p

二、Bender decomposition代码实现

具体的问题形式我就不写出来了，读者可以在旁边用笔写，比我用这该死的latex写起来快一点

python版本：

import gurobipy as gp
from gurobipy import GRB

def benders_decomposition():
    # 创建主问题模型
    master = gp.Model("master_problem")
    
    # 添加二进制变量 y 和连续变量 theta
    y = master.addVars(5, vtype=GRB.BINARY, name="y")
    theta = master.addVar(lb=-GRB.INFINITY, name="theta")
    
    # 设置主问题的目标函数
    master.setObjective(theta + 7 * sum(y[i] for i in range(5)), GRB.MINIMIZE)
    
    # 初始约束，保证theta有初始值
    master.addConstr(theta >= -GRB.INFINITY)
    
    # 初始求解主问题
    master.optimize()
    
    # 初始化Benders割的循环
    iteration = 0
    optimal = False
    while not optimal:
        iteration += 1
        print(f"Iteration {iteration}")
        
        # 获取当前y的值
        y_vals = [y[i].x for i in range(5)]
        
        # 创建子问题模型
        subproblem = gp.Model("subproblem")
        x = subproblem.addVars(5, lb=0, name="x")
        
        # 设置子问题的目标函数
        subproblem.setObjective(gp.quicksum(x[i] for i in range(5)), GRB.MINIMIZE)
        
        # 添加子问题的约束
        subproblem.addConstr(x[0] + x[3] + x[4] == 8, name="C1")
        subproblem.addConstr(x[1] + x[4] == 3, name="C2")
        subproblem.addConstr(x[2] + x[3] == 5, name="C3")
        subproblem.addConstr(x[0] <= 8 * y_vals[0], name="R1")
        subproblem.addConstr(x[1] <= 3 * y_vals[1], name="R2")
        subproblem.addConstr(x[2] <= 5 * y_vals[2], name="R3")
        subproblem.addConstr(x[3] <= 5 * y_vals[3], name="R4")
        subproblem.addConstr(x[4] <= 3 * y_vals[4], name="R5")
        
        # 求解子问题
        subproblem.optimize()
        
        # 检查子问题的可行性
        if subproblem.status == GRB.OPTIMAL:
            sub_obj = subproblem.ObjVal
            
            # 获取对偶变量
            pi1 = subproblem.getConstrByName("C1").Pi
            pi2 = subproblem.getConstrByName("C2").Pi
            pi3 = subproblem.getConstrByName("C3").Pi
            duals = [subproblem.getConstrByName(f"R{i+1}").Pi for i in range(5)]
            
            # 添加Benders割
            master.addConstr(
                theta >= sub_obj - 
                (pi1 * 8 + pi2 * 3 + pi3 * 5) + 
                sum(duals[i] * 8 * y[i] for i in range(5))
            )
        elif subproblem.status == GRB.INFEASIBLE:
            print("Subproblem is infeasible.")
            subproblem.computeIIS()
            subproblem.write("subproblem.ilp")

            # 获取导致不可行性的约束
            infeasible_constraints = [c for c in subproblem.getConstrs() if c.IISConstr]
            
            # 添加可行性割
            feasibility_cut = gp.LinExpr()
            for constr in infeasible_constraints:
                if "R" in constr.ConstrName:
                    index = int(constr.ConstrName[1:]) - 1
                    feasibility_cut += y[index]
            master.addConstr(feasibility_cut >= 1)
        else:
            print("Unexpected subproblem status:", subproblem.status)
            break
        
        # 求解主问题
        master.optimize()
        
        # 检查收敛
        if master.status == GRB.OPTIMAL:
            # 检查主问题目标值和子问题目标值是否收敛
            if theta.x - sub_obj < 1e-6:
                optimal = True
    
    # 打印最优解
    print("Optimal objective value:", master.ObjVal)
    for i in range(5):
        print(f"y_{i} =", y[i].x)

benders_decomposition()

matlab版本：

% 方程案例来自于：哔哩哔哩------课程四 基于Gurobi 的列生成和 Benders 分解方法
%% Benders Decomposition 示例程序

clear variables
close all %关闭所有figure
warning off  %关闭所有警告消息的显示
clc

%% 构建系数矩阵

% 1. 问题的标准形式
% min   cx + dy
% s.t.  Ax + By >= f
%       x >= 0
%       x is continious
%       y is binary

c = [ 1, 1, 1, 1, 1];
d = [ 7, 7, 7, 7, 7];
A = [-1, 0, 0,-1,-1;
      1, 0, 0, 1, 1;
      0,-1, 0, 0,-1;
      0, 1, 0, 0, 1;
      0, 0,-1,-1, 0;
      0, 0, 1, 1, 0;
     -1, 0, 0, 0, 0;
      0,-1, 0, 0, 0;
      0, 0,-1, 0, 0;
      0, 0, 0,-1, 0;
      0, 0, 0, 0,-1];  % 系数矩阵A
B = [ 0, 0, 0, 0, 0;
      0, 0, 0, 0, 0;
      0, 0, 0, 0, 0;
      0, 0, 0, 0, 0;
      0, 0, 0, 0, 0;
      0, 0, 0, 0, 0;
      8, 0, 0, 0, 0;
      0, 3, 0, 0, 0;
      0, 0, 5, 0, 0;
      0, 0, 0, 5, 0;
      0, 0, 0, 0, 3];  % 系数矩阵A
f = [-8; 8;-3; 3;-5; 5; 0; 0; 0; 0; 0];  % 系数矩阵f
%% 采用 yalmip + cplex 直接建模求解

% 1. 变量设置
x = sdpvar(5,1,'full');  % 变量x
y = binvar(5,1,'full');  % 变量y

% 2. 模型构建
obj = c*x + d*y;  % 目标函数
constr = [];  % 约束条件
constr = [constr, A*x + B*y >= f];
constr = [constr, x >= 0];

% 3. 模型求解
opts = sdpsettings('solver','cplex','verbose',2);  % yalmip参数设置
diag = optimize(constr,obj,opts);  % 模型求解
if diag.problem == 0  % 有可行解
    res_cplex = value(obj);  % 记录可行解
    sol_cplex = [value(x),value(y)];
else
    disp('The original problem is infeasible!')
    pause();
end

%% 采用Benders Decomposition求解

% 1. 算法初始化
lb = [];  % 下界
ub = [];  % 上界
sub_coef = [];  % 记录对偶系数
sub_sol  = [];  % 记录子问题中主问题传递的解
sub_obj  = [];  % 记录子问题的目标函数

% 2. 主函数
while true
    % 1) 变量设置
    x = sdpvar(5,1);
    y = binvar(5,1);
    z = sdpvar(size(A,1),1);  % 不可行松弛新增变量
    t = sdpvar(1,1);  % 主问题新增变量
    
    % 2) 求解主问题
    obj = d*y + t;
    constr = [];
    constr = [constr, t >= 0];
    if ~isempty(sub_coef)  % 添加Benders Cut
        for i = 1:size(sub_coef,2)
            constr = [constr, t >= sub_obj(i) - sub_coef(:,i)' * (y - sub_sol(:,i))];
        end
    end
    opts = sdpsettings('solver','cplex','verbose',2);
    optimize(constr,obj,opts);
    y_star = value(y);  % 记录当前y的值
    lb = [lb, value(obj)];  % 更新下界
    
    % 3) 求解子问题
    obj = c*x;
    constr = [];
    constr = [constr, A*x + B*y >= f];
    constr = [constr, x >= 0];
    constr = [constr, y == y_star];  % 将主问题结果应用于子问题
    opts = sdpsettings('solver','cplex','verbose',2,'relax',1);  % 将所有0-1变量松弛
    diag = optimize(constr,obj,opts);
    
    sub_coef = [sub_coef, dual(constr(end))];  % 更新对偶系数
    sub_sol  = [sub_sol , y_star];  % 更新子问题中主问题传递的解
    sub_obj  = [sub_obj , value(obj)];  % 更新子问题的目标函数
    ub = [ub, d * y_star + value(obj)];  % 更新上界
    
    % 4) 收敛判断
    gap = (ub(end) - lb(end)) / ub(end);
    if gap <= 1e-5
        res_benders = ub(end);
        sol_benders = [value(x),value(y)];
        break;
    end
end

%% 输出结果

disp('');
disp('cplex得到的问题最优解为：');
sol_cplex
disp(['最优值：', num2str(res_cplex)]);
disp('');
disp('Benders Decomposition得到的问题最优解为：');
sol_benders
disp(['最优值：', num2str(res_benders)]);

这里matlab代码中的disp(“”)可以用来打印一个空行，因为disp函数执行一次就要换行，但是还有一个东西他执行起来是不换行的，那就是fprintf（相信学过c语言的对这个打印函数是有印象的，因为matlab中的很多源文件的是c写的，所以可以使用，使用时一般要加上换行符\n）

注：赶时间 关于原理很多细节没有写清楚，这里有视频详解：https://www.youtube.com/watch?app=desktop&v=vQzpydNOWDY