Benders decomposition 1解释+pythoncode（不用自行推导dual版本）

已于 2023-02-06 11:14:20 修改
阅读量1.4k 收藏 18
点赞数 4
文章标签： python 线性代数 交通物流 Powered by 金山文档
于 2023-02-06 10:38:35 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_44149045/article/details/128896662
版权

什么是Benders decompostion？

  1. 数学模型解释：

  • MIP 原问题

  • 限制性主问题（主）

  • 子问题（子）

  • 算法流程图：

  1. 算法解释

  • 适用于MIP模型（有整数和连续变量）

  • 将MIP模型进行分解，分解为连续部分（连续变量）和离散部分（整数/01变量）

  • 离散部分作为限制性主问题（主），主问题中的q为变量

  • 连续部分作为子问题（子），子问题中的y_bar由主问题output的y值决定

  • 检测子问题（子）的值infeasible还是optimal，从而找到合适的cut

  • 将cut加入master problem继续求解，直到算法结束

  • 算法结束条件：当主问题的q_bar，和子问题中算出的q_star=cx相同时，算法停止

  1. 代码python code+模块化解释

  • 算列说明：

  • f_i 是fix_cost , c_i是cost，I 是one_matrix 

class bender:
    def __init__(self):
        self.y_num=5
        self.x_num=5
        self.fix_cost=[7 for i in range(self.y_num)]
        self.cost=[1,1,1,1,1]
        self.M = [[1,0,0,1,1],[0,1,0,0,1],[0,0,1,1,0]]
        self.b=[8,3,5]
        self.one_Matrix=list(np.identity(5))
        self.y_rhs=[8,3,5,5,3]
        #########关于master problem的参数#######
        self.q_LB=np.inf
        self.iter=0

  • 整理算法流程main函数部分：

##########################initialization############
mp=master_pro()
mp.optimize()
mp.display()
print('_++++++++++++++++++++++++_')
y_bar,q_bar=update_y_q(mp)
print('_++++++++++++++++++++++++_')
u,q_star,status=dual(y_bar)
mp._sp_status = status
mp._q_star = q_star
mp._q_bar = q_bar
mp=bender_cut_dual(mp,u)
mp.optimize()
mp.display()
y_bar,q_bar=update_y_q(mp)
# # #############start to loop################
while mp._iter<10:
    u, q_star, status = dual(y_bar)
    print('***********sub_prob_q ',q_star)
    print('subdual_{}_\n_{}'.format(u,q_star))
    mp._sp_status = status
    mp._q_star = q_star
    mp._q_bar = q_bar  # q_LB
    if mp._q_star != mp._q_bar:
        new_mp = bender_cut_dual(mp, u)
        new_mp.optimize()
        new_mp.display()
        y_bar,q_bar = update_y_q(new_mp)
        mp=new_mp
    else:
        break

  • 主问题

###############master model,但未添加cut
def master_pro():
    data = bender()
    # construct master problem
    mp = gp.Model('Master')
    mp._iter=0
    y = {}
    for i in range(data.y_num):
            y[i] = mp.addVar(lb=0,ub=1,vtype=GRB.BINARY, name="y" + str(i+1))
    q = mp.addVar(lb=0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name='q')  # q based on different value

    # set objective function
    obj = LinExpr(0)
    for i in range(data.y_num):
            obj.addTerms(data.fix_cost[i], y[i])
    obj+=q
    mp.setObjective(obj, GRB.MINIMIZE)
    #constraint
    # mp.addConstr(q>=0)
    mp.setParam('LazyConstraints', 1)
    mp._q=q
    mp._y=y
    mp.setParam('OutputFlag', 0)
    return mp

  • 拿到y_bar 和q_bar值

#get y_bar
def update_y_q(model):
    y_bar = []
    q_bar=None
    for v in model.getVars():
        if 'y' in v.varName:
            y_bar.append(v.x)

        else:
            q_bar = v.x
    print('q_bar', q_bar)
    print('y_bar', y_bar)
    return y_bar,q_bar

  • 将值带入sub问题中，并拿到extreme ray /extreme point

  • extreme ray：sub infeasible情况，用farkasdual可以拿到值（设置好参数，参数在代码里）

  • extreme point：sub optimal情况，用pi可以拿到值

def sub(y_bar):
    data=bender()
    sp_dual = Model('Sub')
    x = {}
    for i in range(data.x_num):
        x[i] = sp_dual.addVar(vtype=GRB.CONTINUOUS, name="x" + str(i) )
    # set objective
    obj = LinExpr(0)
    for i in range(data.x_num):
        obj.addTerms(data.cost[i], x[i])
    sp_dual.setObjective(obj, GRB.MINIMIZE)

    # set constraints
    # Constraints
    for i in range(3):
        Cons = gp.LinExpr(0)
        for j in range(data.y_num):
            Cons.addTerms(data.M[i][j], x[j])
        sp_dual.addConstr(Cons == data.b[i], name='x' + str(i))
        Cons.clear()

    for i in range(data.x_num):
        Cons = gp.LinExpr(0)
        for j in range(data.y_num):
            Cons.addTerms(data.one_Matrix[i][j], x[j])
        sp_dual.addConstr(Cons <= data.y_rhs[i]*y_bar[i], name='y' + str(i))
        Cons.clear()

    sp_dual.setParam('InfUnbdInfo', 1)
    sp_dual.setParam('Method', 1)
    sp_dual.setParam('Presolve', 0)

    sp_dual.optimize()
    sp_dual.display()
    #check dual varibale

    print("========u=============")
    if sp_dual.status==GRB.INFEASIBLE:
        u = {}
        count = 0
        for c in sp_dual.getConstrs():
            u[count]=c.farkasdual
            count+=1
        print(u)
        q_star = None
        return u, q_star, sp_dual.status

    elif sp_dual.status==GRB.OPTIMAL:
        u = {}
        count = 0
        for c in sp_dual.getConstrs():
            u[count] = c.pi
            count += 1
        print(u)
        q_star = sp_dual.objval
        return u,q_star,sp_dual.status

  • 拿到bender cut，函数bender_cut_dual

def bender_cut_dual(model,u):
    print('______________')

    if model._iter >=0:
        data=bender()
        if model._sp_status==GRB.OPTIMAL:
            if abs(model._q_star-model._q_bar ) > 1e-6:
                print('~~~~~~~~~~~optimal~~~~~~~~~~~~~~~~~~')
                print('iteration:', model._iter)
                print('model._q_star',model._q_star)
                print('model._q_bar',model._q_bar)
                #add feasibility cut:
                lazy=LinExpr(0)
                sum=0
                for j in range(len(u)):
                    if j<3:
                        sum+=u[j]*data.b[j]
                    else:
                        coef=u[j]*data.y_rhs[j-3]
                        lazy.addTerms(coef,model._y[j-3])
                model.addConstr(lazy+sum<=model._q)###not sure q how to iter
                model._iter += 1
            else:
                model.terminate()

        elif model._sp_status == GRB.INFEASIBLE:
            print('~~~~~~~~~~~infeasible~~~~~~~~~~~~~~~~~~')
            print('iteration:', model._iter)
            lazy = LinExpr(0)
            sum = 0
            for j in range(len(u)):
                if j < 3:
                    sum +=np.dot(u[j],data.b[j])
                else:
                    coef = u[j] * data.y_rhs[j - 3]
                    lazy.addTerms(coef, model._y[j - 3])
            model.addConstr(lazy + sum >= 0)  ###not sure q how to iter
            model._iter += 1
        else:
            model.terminate()

        return model

  • 整体代码

import numpy as np
from gurobipy import *
import gurobipy as gp

class bender:
    def __init__(self):
        self.y_num=5
        self.x_num=5
        self.fix_cost=[7 for i in range(self.y_num)]
        self.cost=[1,1,1,1,1]
        self.M = [[1,0,0,1,1],[0,1,0,0,1],[0,0,1,1,0]]
        self.b=[8,3,5]
        self.one_Matrix=list(np.identity(5))
        self.y_rhs=[8,3,5,5,3]
        #########关于master problem的参数#######
        self.q_LB=np.inf
        self.iter=0
        self.dual_Matrix=[[1,0,0,1,0,0,0,0],[0,1,0,0,1,0,0,0],[0,0,1,0,0,1,0,0],[1,0,1,0,0,0,1,0],[1,1,0,0,0,0,0,1]]

###############master model,但未添加cut
def master_pro():
    data = bender()
    # construct master problem
    mp = gp.Model('Master')
    mp._iter=0
    y = {}
    for i in range(data.y_num):
            y[i] = mp.addVar(lb=0,ub=1,vtype=GRB.BINARY, name="y" + str(i+1))
    q = mp.addVar(lb=0,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name='q')  # q based on different value

    # set objective function
    obj = LinExpr(0)
    for i in range(data.y_num):
            obj.addTerms(data.fix_cost[i], y[i])
    obj+=q
    mp.setObjective(obj, GRB.MINIMIZE)
    #constraint
    # mp.addConstr(q>=0)
    mp.setParam('LazyConstraints', 1)
    mp._q=q
    mp._y=y
    mp.setParam('OutputFlag', 0)
    return mp

def bender_cut_dual(model,u):
    print('______________')

    if model._iter >=0:
        data=bender()
        if model._sp_status==GRB.OPTIMAL:
            if abs(model._q_star-model._q_bar ) > 1e-6:
                print('~~~~~~~~~~~optimal~~~~~~~~~~~~~~~~~~')
                print('iteration:', model._iter)
                print('model._q_star',model._q_star)
                print('model._q_bar',model._q_bar)
                #add feasibility cut:
                lazy=LinExpr(0)
                sum=0
                for j in range(len(u)):
                    if j<3:
                        sum+=u[j]*data.b[j]
                    else:
                        coef=u[j]*data.y_rhs[j-3]
                        lazy.addTerms(coef,model._y[j-3])
                model.addConstr(lazy+sum<=model._q)###not sure q how to iter
                model._iter += 1
            else:
                model.terminate()

        elif model._sp_status == GRB.INFEASIBLE:
            print('~~~~~~~~~~~infeasible~~~~~~~~~~~~~~~~~~')
            print('iteration:', model._iter)
            lazy = LinExpr(0)
            sum = 0
            for j in range(len(u)):
                if j < 3:
                    sum +=np.dot(u[j],data.b[j])
                else:
                    coef = u[j] * data.y_rhs[j - 3]
                    lazy.addTerms(coef, model._y[j - 3])
            model.addConstr(lazy + sum >= 0)  ###not sure q how to iter
            model._iter += 1
        else:
            model.terminate()

        return model
#get y_bar
def update_y_q(model):
    y_bar = []
    q_bar=None
    for v in model.getVars():
        if 'y' in v.varName:
            y_bar.append(v.x)

        else:
            q_bar = v.x
    print('q_bar', q_bar)
    print('y_bar', y_bar)
        # print('jjjjjjjjjjjjjjjjj',y_bar)
    return y_bar,q_bar


def dual_varible(model):
    if model.status==GRB.INFEASIBLE:
        u={}
        count=0
        for v in model.getVals():
            u[count]=v.unbdRay
        print(u)
    elif model.status==GRB.OPTIMAL:
        u={}
        count=0
        for v in model.getVars():
            u[count]=v.pi
        return u

def sub(y_bar):
    data=bender()
    sp_dual = Model('Sub')
    x = {}
    for i in range(data.x_num):
        x[i] = sp_dual.addVar(vtype=GRB.CONTINUOUS, name="x" + str(i) )
    # set objective
    obj = LinExpr(0)
    for i in range(data.x_num):
        obj.addTerms(data.cost[i], x[i])
    sp_dual.setObjective(obj, GRB.MINIMIZE)

    # set constraints
    # Constraints
    for i in range(3):
        Cons = gp.LinExpr(0)
        for j in range(data.y_num):
            Cons.addTerms(data.M[i][j], x[j])
        sp_dual.addConstr(Cons == data.b[i], name='x' + str(i))
        Cons.clear()

    for i in range(data.x_num):
        Cons = gp.LinExpr(0)
        for j in range(data.y_num):
            Cons.addTerms(data.one_Matrix[i][j], x[j])
        sp_dual.addConstr(Cons <= data.y_rhs[i]*y_bar[i], name='y' + str(i))
        Cons.clear()

    sp_dual.setParam('InfUnbdInfo', 1)
    sp_dual.setParam('Method', 1)
    sp_dual.setParam('Presolve', 0)

    sp_dual.optimize()
    sp_dual.display()
    #check dual varibale

    print("========u=============")
    if sp_dual.status==GRB.INFEASIBLE:
        u = {}
        count = 0
        for c in sp_dual.getConstrs():
            u[count]=c.farkasdual
            count+=1
        print(u)
        q_star = None
        return u, q_star, sp_dual.status

    elif sp_dual.status==GRB.OPTIMAL:
        u = {}
        count = 0
        for c in sp_dual.getConstrs():
            u[count] = c.pi
            count += 1
        print(u)
        q_star = sp_dual.objval
        return u,q_star,sp_dual.status






##########################initialization############
mp=master_pro()
mp.optimize()
mp.display()
print('_++++++++++++++++++++++++_')
y_bar,q_bar=update_y_q(mp)
print('_++++++++++++++++++++++++_')
# u,q_star,status=sub_dual(y_bar)
u,q_star,status=sub(y_bar)
mp._sp_status = status
mp._q_star = q_star
mp._q_bar = q_bar
# mp=bender_cut(mp,u)
mp=bender_cut_dual(mp,u)
mp.optimize()
mp.display()
y_bar,q_bar=update_y_q(mp)
# # #############start to loop################
while mp._iter<10:
    u, q_star, status = sub(y_bar)
    # u,q_star,status=sub_dual(y_bar)
    print('***********sub_prob_q ',q_star)
    print('subdual_{}_\n_{}'.format(u,q_star))
    mp._sp_status = status
    mp._q_star = q_star
    mp._q_bar = q_bar  # q_LB
    if mp._q_star != mp._q_bar:
        # new_mp = bender_cut(mp,u)
        new_mp = bender_cut_dual(mp, u)
        new_mp.optimize()
        new_mp.display()
        y_bar,q_bar = update_y_q(new_mp)
        mp=new_mp
    else:
        break

  • bender decompostion2，根据推导并构建的sub-dual来求解

  • 最终结果展示

  • 主函数

  • 子函数

奶伦呀❤️
关注
  • 4
    点赞
  • 踩
  • 18
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
matlab2016代码-NIHBA:网络遮断概念的应变设计，并通过混合Benders分解算法求解
05-24
matlab2016代码美国国立卫生研究院 NIHBA（带有混合Benders算法的网络拦截）是一个软件平台，用于识别代谢工程中尽可能多的遗传修饰策略。 它基于基因组规模的代谢模型（GSMM），并将宿主细胞和代谢工程师视为玩游戏的两个不同年龄网络：代谢工程师试图最大程度地违反宿主细胞的主要目标。 这种游戏在博弈论中被称为网络拦截。 NIHBA采用混合Bender分解算法来解决此网络拦截问题，从而在没有数字问题的情况下实现最佳设计策略，从而实现了高效率。 入门 这些说明将为您提供在本地计算机上运行并运行的项目的副本，以进行开发和测试。 先决条件 MATLAB2016和更高版本的cobratoolbox（最新版本：3.0）Gurobi7.5（或包含LP和ILP求解器的其他软件包） 安装 1，直接从该存储库下载整个NIHBA软件包或使用$ git clone --depth = 1 NIHBA进行克隆 2，启动MATLAB并将工作目录更改为NIHBA作为当前文件夹。 然后运行nihba_setup nihba_setup 注意：nihba_setup将帮助您将所有NIHBA文件添加到matl
线性规划技巧: Benders Decomposition
胡拉哥
03-19 1万+
Benders分解由Jacques F. Benders在1962年提出1. 它是一种把线性规划问题分解为小规模子问题的技巧. 通过迭代求解主问题和子问题, 从而逼近原问题的最优解. 与列生成(Column Generation)相比, Benders分解是一种行生成(Row Generation)技巧: 主问题的约束来自子问题的解. 本文介绍的内容基于 Georrion和Graves2. 问题描...
关于Bender decomposition分解算法解析附python代码
最新发布
2301_79911052的博客
06-10 970
最近在学习鲁棒优化的时候看到了一个新算法也就是本文的Bender 分解算法 本着求知的态度，看了网上的众多解释可教程，大概算是有点了解了，顺手写个blog记录一下吧，附上复现的python代码和matlab代码。
运筹系列7：大规模线性规划的Benders/DW分解
kittyzc的博客
04-04 3万+
1. 行生成算法 行生成就是指的不断添加约束的算法。 因为在求解矩阵中，一个约束条件对应一行，因此添加约束条件的方法自然叫做行生成算法。相对应的，添加变量的方法就叫做列生成算法。 这一节先看行生成算法，用在求解变量不多，但是约束条件特别多的情况下。 2. Benders分解 Benders分解（Benders Decomposition，BD）的基本思路是：使用子问题（primal ...
Benders Decomposition初认识
bkave的博客
11-24 2109
最近在阅读文献时，注意到一种求解混合整数规划问题（MIP）的方法——Benders Decomposition，先写下一点初步的学习理解。
Benders分解算法详解
热门推荐
qq_42756072的博客
08-01 3万+
前言 最近在学习Benders Decomposition，看了不少的文章博客，但我感觉写的都不是很清楚，细节地方解释的也不够到位，于是就去读了英文论文原文，发现原文写的还是相当清楚的。本文也基于那篇论文，对Benders Decomposition进行总结。 背景 Benders Decomposition由Jacques F. Benders提出，用来求解混合整数规划问题(MIP:Mixed Integer Programming problem)，即同时包含整数和连续变量的极值问题，换句话说就是线性规
详解 Benders 分解与一个算例的 python 代码
天天向上的专栏
11-15 6839
听说过 benders 分解几年了，一直没有细看，最近论文里面也见到，还是有必要了解一下它的基本思想与用法的。
Generalized Benders Decomposition 学习笔记论文对照版）
05-07
针对Generalized Benders Decomposition原始论文逐步进行推导，直到推导出Generalized Benders Decomposition的一般解法
Benders_Decomposition.rar_Benders算法_ROBUST benders_bende_benders
07-14
一个好的最短路径算法参考文献——A benders decomposition approach
Benders-Decomposition
03-16
**Benders分解法详解** Benders分解法，也称为Benders策略或Benders分解，是运筹学领域中解决大规模线性规划问题的一种有效方法。该方法由J.F. Benders于1962年提出，主要用于将一个复杂的优化问题分解为两个相对较...
基于MATLAB实现的benders分解算法+使用说明文档.zip
05-24
基于MATLAB实现的benders分解算法+使用说明文档.zip 1、代码压缩包内容 主函数：main.m； 调用函数：其他m文件；无需运行 运行结果效果图； 2、代码运行版本 Matlab 2020b；若运行有误，根据提示GPT修改；若不会，...
BendersDecomposition.rar_BendersDecomposition_benders_benders分解_
07-15
**Benders分解法详解** Benders分解法，也称为Benders策略或Benders分解，是一种在运筹学和优化领域广泛使用的求解方法，尤其适用于大规模的混合整数线性规划（MILP）问题。这种方法由荷兰数学家J.F. Benders于1962...
benders分解算法 逻辑思路整理（加星）
sinat_41348401的博客
05-25 4067
一般我们经常可以看到的是经典的benders，此处我也结合经典的benders梳理下思路。
Gurobi运筹学开发教程03：python+gurobi实现bender decomposition 算法示例代码
m1m2m3mmm的博客
05-10 1995
一个简单的bender分解示例
Lecture09:Benders分解--理论部分
AmieeXue的博客
01-05 5332
Benders分解的理论推导
Benders分解求固定花费的运输问题
xiaoyao_de的博客
07-27 798
参考数据魔术师的一篇Benders分解算法(附代码及注释) 添加链接描述 该算法在求解模型时子问题用的是对偶模型。如果子问题限制特别多，可以在不用对偶模型的条件下完成该算法。用cplex中的getDual语句可得到对偶模型的极点与极射线，同时修改需要子问题的限制(上下界用setUB,LB)与目标函数，也可以用来求解。 ...
运筹学读书笔记---Benders Decomposition
IELBHJY的博客
04-12 3487
1、前言 最近在读这本厚厚的运筹学书籍，发现这本书写的真的好。大规模优化方法系列，有个叫benders分解的算法，之前只是听过名字，不懂具体原理。仔细啃了一下，还是很有意思的。 2、分解策略 如果y值固定，上述问题就变成一个线性规划问题。因此，在第l次迭代时。选取部分整数决策变量y(l)并固定相应的值，从而将整数规划转化为只包含连续决策...
Benders分解(Benders Decomposition)算法：算法原理+具体算例
weixin_42991799的博客
09-12 9595
Lecture10a：Benders分解 -- 应用部分
AmieeXue的博客
01-06 4021
benders分解方法、算法和流程
benders decomposition
03-16
Benders分解是一种用于解决大规模线性规划问题的分解技术。它将原问题分解成一个主问题和若干个子问题，通过求解子问题来逐步优化主问题的解。Benders分解的优点是可以有效地处理大规模问题，并且可以利用问题的特殊结构来提高求解效率。
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值