什么是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来求解

  • 最终结果展示

  • 主函数

  • 子函数

# python # 线性代数 # 交通物流
Logo
腾讯云开发者社区

腾讯云面向开发者汇聚海量精品云计算使用和开发经验，营造开放的云计算技术生态圈。

更多推荐

cover

动物识别系统 python实现+UI

avatar 腾讯云开发者社区

Linux下安装MySQL8.0(超详细)

Linux下安装MySQL8.0(超详细)

avatar 腾讯云开发者社区
cover

H3C多WAN口负载均衡实验

avatar 腾讯云开发者社区