从列生成，分支定界到分支定价算法解决大规模一维下料问题（一）

最近在调研大规模的cutting stock问题，发现了分支定价(price and bound) 方法。但是网上关于分支定价(price and bound)的资料很少，对自己的学习过程做个记录，希望能帮助到有需要的小伙伴。

第一部分主要讲cutting stock问题的定义，以及如何用列生成算法解决cutting stock问题。

一. cutting stock问题的定义

问题描述：有一批钢管数量为 N，长度都为 L；假设m个顾客，来钢厂下订单，每位顾客需要长度为 Li的钢管 Di根。

优化目标: 如何切割钢管使得既能满足客户的需求又使得的消耗的钢管数量最少？

整数规划模型如下

二. 列生成算法

列生成算法的主要流程如下：

step1：给出模型的初始解，其解代表可行列

step2: 求解主问题（RMP）并获取乘子π

step3: 求解子问题 ，并计算检验数，如果小于0，向限制主问题中添加一列求解，返回step2

step4: 若检验数大于等于零，停止算法，输出最优解

初始原问题时一个集合覆盖模型，一开始列出所有的列是不现实的，因此我们一开始只能从有限列出发，不断添加列，直到包含有限列的模型（限制主问题）与原问题模型具有相同的线性最优解；为了使用列生成，必须将原问题松弛为线性模型，因为列生成只能对线性模型使用。

python代码如下：

import gurobipy as grb
from gurobipy import GRB
import numpy as np

# column就是pattern
def master_problem(column, vtype, demand_number_array):
    m = grb.Model()
    x = m.addMVar(shape=column.shape[1], lb=0, vtype=vtype)
    m.addConstr(lhs=column @ x >= demand_number_array)
    m.setObjective(x.sum(), GRB.MINIMIZE)
    m.optimize()

    if vtype == GRB.CONTINUOUS:
        return np.array(m.getAttr('Pi', m.getConstrs()))
    else:
        return m.objVal, np.array(m.getAttr('X'))

# 将主问题松弛成线性问题
def restricted_lp_master_problem(column,demand_number_array):
    return master_problem(column, GRB.CONTINUOUS,demand_number_array)

# 找到所有的cut pattern后，求解整数优化主问题。
def restricted_ip_master_problem(column,demand_number_array):
    return master_problem(column, GRB.INTEGER,demand_number_array)

# 子問題求解，寻找pattern
def knapsack_subproblem(kk,demand_width_array,demand_number_array,roll_width):
    m = grb.Model()
    x = m.addMVar(shape=kk.shape[0], lb=0, vtype=GRB.INTEGER)  # 整数规划
    m.addConstr(lhs=demand_width_array @ x <= roll_width)
    m.setObjective(1 - kk @ x, GRB.MINIMIZE)
    m.optimize()

    flag_new_column = m.objVal < 0   # 判别数
    if flag_new_column:
        new_column = m.getAttr('X')  # 子问题求解，找到新的列
    else:
        new_column = None
    return flag_new_column, new_column

def test(roll_width,demand_width_array,demand_number_array):
    # 初始定义最简单的pattern
    initial_cut_pattern = np.diag(np.floor(roll_width / demand_width_array))
    flag_new_cut_pattern = True
    new_cut_pattern = None
    cut_pattern = initial_cut_pattern  # 最近的pattern
    while flag_new_cut_pattern:
        if new_cut_pattern:
            cut_pattern = np.column_stack((cut_pattern, new_cut_pattern))    # 最新的pattern。按列合并
        kk = restricted_lp_master_problem(cut_pattern, demand_number_array)  # 将主问题松弛成线性问题
        flag_new_cut_pattern, new_cut_pattern = knapsack_subproblem(kk, demand_width_array, demand_number_array, roll_width)
    # 找到所有的cut pattern后，求解整数优化主问题。
    minimal_stock, optimal_number = restricted_ip_master_problem(cut_pattern, demand_number_array)
    print(f'demand_width_array: {demand_width_array}')
    print(f'demand_number_array: {demand_number_array}')
    print('result:')
    print(f'minimal_stock: {minimal_stock}')

# 生成测试数据
def generate_data(roll_width,customer_num):
    roll_width = np.array(roll_width)
    demand_width_array = np.random.randint(1, 3, size=(customer_num))
    demand_number_array = np.random.randint(50, 200, size=(customer_num))
    return roll_width, demand_width_array, demand_number_array

# print(generate_data(120,100))
roll_width, demand_width_array, demand_number_array = generate_data(roll_width=120,customer_num=2)
test(roll_width, demand_width_array, demand_number_array)