📢作者： 小小明-代码实体
仓库选址问题
仓库选址是运输行业最常见的问题，目标是使整个配送系统的总成本最低。
例如有100个网点，7个备选仓库位置，要选从7个备选点选出使总成本最低的位置作为仓库。
每个备选点所需的日均运营成本和每天能够处理的最大需求量如下表：
每个网点每天有固定的需求量，各网点到指定备选点的成本也固定，上图展示了部分网点的日需求量和其到达各备选点的成本，下面我们使用OR-Tools建模并求解这个问题。
首先我们准备数据，各网点运输至各选址点的成本数据如下：
import pandas as pd
cost = [[7, 7, 4, 3, 8, 2, 6, 4, 10, 2, 10, 7, 3, 4, 8, 10, 6, 5, 6, 3, 10, 9, 4, 2, 6, 1, 7, 2, 9, 3, 8, 3, 5, 8, 9, 1, 7, 9, 7, 5, 1, 9, 1, 5, 5, 10, 8, 7, 5, 9, 9, 1, 4, 5, 1, 7, 7, 10, 5, 8, 2, 8, 4, 7, 1, 5, 10, 7, 5, 6, 5, 6, 6, 1, 2, 2, 9, 8, 2, 10, 4, 1, 9, 3, 3, 4, 3, 3, 10, 9, 7, 10, 1, 9, 2, 6, 2, 5, 3, 5],
[5, 8, 1, 5, 1, 10, 8, 1, 3, 3, 6, 3, 1, 10, 1, 6, 3, 5, 5, 9, 4, 9, 1, 7, 7, 8, 8, 5, 9, 10, 9, 6, 5, 4, 3, 4, 3, 2, 4, 10, 6, 7, 5, 4, 6, 10, 5, 7, 3, 7, 8,
9, 9, 8, 10, 5, 10, 9, 5, 1, 7, 10, 2, 3, 9, 2, 5, 8, 10, 3, 9, 5, 8, 10, 2, 9, 3, 6, 4, 5, 9, 1, 6, 9, 1, 5, 1, 8, 10, 3, 10, 8, 3, 6, 9, 2, 3, 8, 6, 6],
[9, 4, 9, 5, 4, 6, 4, 4, 9, 3, 10, 2, 6, 10, 1, 2, 3, 2, 8, 10, 4, 6, 10, 10, 9, 6, 6, 2, 9, 3, 7, 7, 6, 10, 4, 4, 8, 7, 5, 5, 8, 3, 1, 3, 7, 3, 6, 3, 3, 2,
8, 1, 6, 6, 5, 9, 4, 9, 4, 9, 1, 10, 10, 2, 8, 10, 4, 2, 1, 4, 3, 7, 7, 8, 9, 2, 3, 6, 8, 1, 10, 9, 6, 5, 10, 1, 1, 10, 2, 7, 3, 1, 5, 3, 10, 9, 8, 4, 6, 9],
[8, 7, 10, 4, 3, 5, 8, 2, 3, 2, 2, 7, 1, 4, 8, 4, 8, 3, 2, 1, 3, 2, 6, 5, 5, 9, 8, 3, 1, 2, 9, 5, 9, 2, 9, 1, 10, 2, 2, 3, 1, 5, 8, 3, 3, 4, 7, 2, 2, 6, 6,
6, 3, 6, 1, 7, 7, 8, 2, 4, 5, 6, 6, 3, 6, 6, 3, 5, 6, 4, 3, 2, 3, 3, 3, 3, 10, 10, 3, 10, 6, 9, 2, 7, 10, 8, 3, 5, 2, 2, 2, 8, 10, 1, 5, 10, 2, 3, 4, 4],
[4, 4, 9, 7, 8, 7, 1, 8, 10, 5, 10, 1, 1, 6, 8, 5, 5, 9, 7, 9, 1, 8, 2, 1, 1, 9, 6, 6, 8, 2, 5, 3, 6, 10, 6, 1, 4, 1, 9, 5, 9, 3, 2, 10, 6, 5, 1, 10, 8, 8,
9, 6, 9, 8, 7, 1, 3, 8, 4, 6, 7, 1, 3, 1, 3, 5, 3, 9, 7, 6, 2, 8, 7, 8, 2, 10, 5, 10, 7, 4, 8, 7, 7, 6, 9, 4, 1, 5, 1, 7, 1, 8, 1, 6, 7, 4, 1, 1, 3, 6],
[4, 3, 8, 10, 3, 5, 1, 7, 5, 7, 10, 5, 8, 2, 8, 2, 5, 2, 6, 5, 3, 10, 10, 6, 9, 5, 10, 8, 2, 2, 4, 10, 1, 3, 7, 1, 10, 3, 8, 7, 1, 4, 1, 3, 8, 2, 10, 5, 7, 2,
5, 4, 1, 2, 6, 4, 9, 3, 4, 10, 7, 9, 7, 4, 1, 9, 5, 9, 10, 7, 9, 4, 8, 5, 8, 3, 10, 6, 2, 3, 1, 2, 10, 1, 6, 7, 5, 6, 4, 2, 7, 4, 4, 6, 6, 1, 9, 6, 5, 8],
[6, 9, 1, 7, 5, 7, 5, 5, 4, 8, 1, 8, 6, 7, 6, 1, 4, 5, 7, 8, 2, 5, 10, 5, 8, 7, 1, 9, 5, 8, 4, 1, 8, 9, 10, 7, 4, 3, 10, 7, 1, 1, 2, 2, 6, 5, 2, 4, 5, 10, 3, 2, 1, 1, 4, 9, 10, 5, 6, 7, 10, 9, 4, 10, 8, 7, 10, 8, 1, 9, 4, 7, 7, 6, 4, 9, 4, 4, 3, 10, 8, 3, 1, 2, 7, 7, 8, 10, 1, 2, 7, 6, 7, 1, 4, 7, 10, 2, 9, 8]]
cost = pd.DataFrame(cost, index=list("ABCDEFG"), columns=range(1, 101)) 
每个选址点的运营成本和最大容量：
data = pd.DataFrame({'运营成本': [90, 40, 50, 60, 55, 123, 70],
'最大容量': [100, 50, 100, 100, 150, 100, 100]},
index=list("ABCDEFG"))
# 每个网点的需求量
need = [5, 3, 2, 5, 5, 2, 4, 3, 4, 4, 2, 3, 2, 4, 3, 2, 2, 4, 3, 4, 2, 2, 2, 2, 3, 2, 4,
5, 5, 3, 5, 4, 5, 2, 4, 5, 2, 2, 4, 3, 2, 4, 3, 4, 3, 2, 4, 3, 5, 4, 5, 3, 3, 4,
3, 4, 3, 3, 3, 3, 4, 4, 2, 3, 3, 4, 3, 5, 2, 4, 4, 3, 4, 3, 3, 5, 3, 2, 5, 3, 3,
5, 3, 3, 3, 4, 3, 5, 3, 3, 2, 3, 2, 3, 3, 3, 2, 2, 4, 5] 下面是完整的处理代码：
from ortools.sat.python import cp_model
model = cp_model.CpModel()
# 选择的仓库
use_facility = cost.index.to_series().apply(model.NewBoolVar)
# 每个网点选择的仓库
ser_customer = cost.applymap(str).applymap(model.NewBoolVar)
# 每个网点选一个地址
for c in ser_customer.columns:
model.Add(ser_customer[c].sum() == 1)
# 保证选择的仓库才有数据
for row, facility in zip(ser_customer.values, use_facility.values):
for e in row:
model.Add(e <= facility)
# 每个地址的最大容量
tmp = (ser_customer*need).sum(axis=1).values
for r, c in zip(tmp, data.最大容量.values):
model.Add(r <= c)
# 最小化运输成本 + 仓库运营成本
model.Minimize((ser_customer*cost).sum().sum() +
(use_facility*data.运营成本.values).sum())
solver = cp_model.CpSolver()
status = solver.Solve(model)
print(solver.StatusName(status), "总成本：", solver.ObjectiveValue())
mask = use_facility.apply(solver.Value) == 1
print("选择的仓库：", use_facility[mask].index.to_list())
result = ser_customer.applymap(solver.Value)[mask]
for c, s in result.T.iteritems():
print("仓库", c, "网点:", s[s != 0].index.to_list())
print("各仓库日均处理的需求：")
needs = (result*need).sum(axis=1)
print(needs)
print("总需求:", needs.sum()) OPTIMAL 总成本： 447.0
选择的仓库： ['B', 'C', 'D', 'E']
仓库 B 网点: [3, 5, 8, 23, 33, 35, 37, 60, 63, 66, 70, 82, 85, 90, 96]
仓库 C 网点: [2, 4, 15, 16, 17, 18, 26, 27, 28, 42, 43, 44, 46, 50, 52, 54, 61, 68, 69, 76, 77, 78, 80, 84, 86, 87, 92]
仓库 D 网点: [6, 9, 10, 11, 14, 19, 20, 22, 29, 34, 39, 40, 41, 45, 48, 49, 51, 53, 55, 59, 72, 73, 74, 79, 81, 83, 94, 95, 100]
仓库 E 网点: [1, 7, 12, 13, 21, 24, 25, 30, 31, 32, 36, 38, 47, 56, 57, 58, 62, 64, 65, 67, 71, 75, 88, 89, 91, 93, 97, 98, 99]
各仓库日均处理的需求：
B 50
C 91
D 98
E 94
dtype: int64
总需求: 333 从结果可以看到选择BCDE这四个仓库时，总成本最低。以上代码的逻辑都有完整的注释，唯一的建模难点在于选择仓库，即如下几行代码：
for row, facility in zip(ser_customer.values, use_facility.values):
for e in row:
model.Add(e <= facility) 上面的最优解选择了四个仓库，我们还可以增加约束来获取选择3个以下仓库和5个以上仓库时的最优解，例如获取选择3个以下仓库时的最优解，增加如下约束：
model.Add(use_facility.sum() <= 3) 重新求解得到：
OPTIMAL 总成本： 448.0
选择的仓库： ['C', 'D', 'E']
仓库 C 网点: [4, 5, 6, 15, 16, 17, 18, 26, 28, 35, 43, 44, 46, 49, 50, 52, 54, 61, 68, 69, 70, 76, 77, 78, 80, 84, 86, 87, 92]
仓库 D 网点: [8, 9, 10, 11, 14, 19, 20, 22, 29, 34, 39, 40, 41, 45, 48, 51, 53, 55, 59, 60, 72, 73, 74, 79, 81, 83, 90, 94, 95, 100]
仓库 E 网点: [1, 2, 3, 7, 12, 13, 21, 23, 24, 25, 27, 30, 31, 32, 33, 36, 37, 38, 42, 47, 56, 57, 58, 62, 63, 64, 65, 66, 67, 71, 75, 82, 85, 88, 89, 91, 93, 96, 97, 98, 99]
各仓库日均处理的需求：
C 100
D 100
E 133
dtype: int64
总需求: 333 选择5个以上仓库时的最优解，增加如下约束：
model.Add(use_facility.sum() >= 5) 重新求解：
OPTIMAL 总成本： 482.0
选择的仓库： ['B', 'C', 'D', 'E', 'G']
仓库 B 网点: [5, 8, 17, 23, 33, 35, 37, 60, 63, 66, 70, 75, 82, 85, 96]
仓库 C 网点: [15, 18, 26, 28, 43, 46, 50, 52, 61, 68, 69, 76, 77, 80, 86, 87, 92]
仓库 D 网点: [4, 6, 9, 10, 13, 14, 19, 20, 22, 29, 30, 34, 36, 39, 40, 41, 45, 48, 49, 55, 59, 72, 73, 74, 81, 90, 94, 100]
仓库 E 网点: [1, 2, 7, 12, 21, 24, 25, 38, 47, 56, 57, 62, 64, 65, 67, 71, 88, 89, 91, 93, 97, 98, 99]
仓库 G 网点: [3, 11, 16, 27, 31, 32, 42, 44, 51, 53, 54, 58, 78, 79, 83, 84, 95]
各仓库日均处理的需求：
B 50
C 58
D 95
E 72
G 58
dtype: int64
总需求: 333 可以看到只选择3个仓库时就能得到接近于最优解的结果，所以实际选址时往往会选择CDE这三个仓库。
钢管切割问题
建模思路1
假设需要长度为3米、7米、9米、16米的钢管各25根、30根、14根、8根，目前只有长度为20米的钢管若干，如何切割使消耗钢管的数量最少？
在1939 年前苏联经济学家 Kantorovich 给出了一个建模模型：
下面我们使用OR-tools对这种模型建模：
from ortools.sat.python import cp_model
import numpy as np
import pandas as pd
L = 20
df = pd.DataFrame(
{"len": [3, 7, 9, 16],
"nums": [25, 30, 14, 8]})
model = cp_model.CpModel()
N = int(df.nums.sum())
max_cut = int(L // df.len.min())
print(N, max_cut)
nums = pd.DataFrame(np.empty((N, df.shape[0])), columns=df.len).applymap(
lambda x: model.NewIntVar(0, max_cut, "x"))
chooses = pd.Series([model.NewBoolVar("y") for _ in range(N)])
# 选择最少的根数
model.Minimize(chooses.sum())
# 切割结果大于需要的根数
for x_j, d_j in zip(nums.sum().values, df.nums.values):
model.Add(x_j >= d_j)
# 切割方案不得大于钢材长度
tmp = nums.mul(df.len.values, axis=1).sum(axis=1).values
for x_j, y_j in zip(tmp, (L*chooses).values):
model.Add(x_j <= y_j)
solver = cp_model.CpSolver()
status = solver.Solve(model)
print(solver.StatusName(status))
value = int(solver.ObjectiveValue())
print("总用料：", value) 77 6
OPTIMAL
总用料： 30 在2秒左右的时间得到了结果，可以看到在仅切割30根钢管的情况下可以满足客户需求。
整理一下结果进行展示：
mask = chooses.apply(solver.Value) == 1
result = nums[mask].applymap(solver.Value)
result = result.value_counts(sort=False).rename("根数").reset_index()
result.sort_values([*result.columns], ascending=False, inplace=True)
tmp = result[df.len]
result["每根剩余"] = L-tmp.mul(df.len.values, axis=1).sum(axis=1)
display(result)
count = result.根数.sum()
yl_sum = (result.每根剩余*result.根数).sum()
usage = 1-yl_sum/(count*L)
print(f"原材料长度：{L}，合计用料：{count}，剩余料头：{yl_sum}，利用率：{usage:.2%}")
need = df.set_index("len").nums.convert_dtypes()
need.name = "需求"
report = result.iloc[:, :-2].mul(result.根数, axis=0).sum()
report.name = "切割结果"
report = pd.concat([need, report], axis=1)
display(report) 在上述条件下，思路1还能计算出结果，但是当需求情况稍微复杂一点时，却几小时也无法计算出结果，显然这是无法接受的，于是我们继续研究其他的解决方案。
建模思路2
上述模型对于稍微复杂一点点的模型计算速度都非常慢，我随便测试了一个加刀缝的案例耗时3小时也未计算出结果，所以只能放弃。例如优化目标如下：
下面我们的建模模型如下：
i表示第几种钢管，j表示第几种切割方案，yj表示每种切割方案所采用的根数，mij表示第i种钢管在第j种切割方案下切割根数，ni表示每种需求钢管的需求数。
要建立上述模型，我们首先需求找出所有的切割方案，ortools可以很方便的求出全部可行解。
先准备好数据：
import pandas as pd
L, knife_seam = 6000, 3
df = pd.DataFrame(
{"len": [2500, 1200, 1000, 500, 400],
"nums": [15, 35, 35, 25, 25]})
N = int(df.nums.sum())
print("总需求根数:", N) 总需求根数: 135 然后我们使用ortools找出所有的切割方案，具体建模思路就是，每根钢管切割的余料必须大于等于0并且小于最短的钢管需求长度，因为一根钢管的剩余长度大于最小需求长度就说明还可以继续切分。
最终代码为：
from ortools.sat.python import cp_model
class MyCpSolver(cp_model.CpSolverSolutionCallback):
def __init__(self, nums, result=None):
if result is None:
result = []
cp_model.CpSolverSolutionCallback.__init__(self)
self.nums = nums
self.result = result
def on_solution_callback(self):
self.result.append(self.nums.apply(self.Value).values)
result = []
model = cp_model.CpModel()
nums = df.len.apply(lambda length: model.NewIntVar(
0, L // (length+knife_seam), "x"))
yl = L-((df.len+knife_seam)*nums).sum()
model.Add(yl < df.len.min())
model.Add(yl >= 0)
solver = cp_model.CpSolver()
myCpSolver = MyCpSolver(nums, result)
solver.parameters.enumerate_all_solutions = True
# 限制最大运行时间，防止查找的结果数过多
solver.parameters.max_time_in_seconds = 2
status = solver.Solve(model, myCpSolver)
print(solver.StatusName(status), f"共{len(result)}种切分方案：")
data = pd.DataFrame(result, columns=df.len)
data["每根剩余"] = L-data.mul((df.len+knife_seam).values, axis=1).sum(axis=1)
data.sort_values(by=[*df.len], ascending=False,
inplace=True, ignore_index=True) 可以看到100毫秒左右就找出全部的136种切分方案。
注意：若需求的规格种类非常多，一定要设置solver.parameters.max_time_in_seconds参数，防止找出的切割方案太多，导致后续难以计算出结果。
然后我们继续对这136种切分方案建模，设置每种方案的采用数量，最小化总根数：
model = cp_model.CpModel()
# 每种切割方案切割的根数
y = data.index.to_series().apply(
lambda i: model.NewIntVar(0, N, f'y_{i}'))
# 切出的数量需要满足需求
tmp = data.iloc[:, :-1].mul(y, axis=0).sum()
for row, num in zip(tmp, df.nums.values):
model.Add(row >= num)
# 最小化总根数
model.Minimize(y.sum())
solver = cp_model.CpSolver()
status = solver.Solve(model)
print(solver.StatusName(status))
data["根数"] = y.apply(solver.Value)
result = data[data.根数 != 0].copy()
count = result.根数.sum()
print("总根数：", count, "，采用的切割方案数：", result.shape[0])
display(result)
yl_sum = (result.每根剩余*result.根数).sum()
usage = 1-yl_sum/(count*L)
total_need = (df.nums*df.len).sum()
total = L*result.根数.sum()
print(f"原材料长度：{L}，刀缝：{knife_seam}")
print(f"剩余料头：{yl_sum}，需求利用率：{total_need/total:.2%}，结果利用率：{usage:.2%}")
need = df.set_index("len").nums.convert_dtypes()
need.name = "需求"
report = result.iloc[:, :-2].mul(result.根数, axis=0).sum()
report.name = "切割结果"
report = pd.concat([need, report], axis=1)
report["多出根数"] = report.切割结果-report.需求
display(report) 然后几十毫秒的时间计算出了一种总根数消耗最少的切割方案：
现在我们已经知道了该组合下，最小的总根数为24，我们可以考虑在限制根数为24的情况下，最大化切割结果的总长度，这样我们剩余的余料会更加集中，更容易被后期利用。
编码如下：
model = cp_model.CpModel()
y = data.index.to_series().apply(
lambda i: model.NewIntVar(0, N, f'y_{i}'))
# 每种切割方案切割的根数
tmp = data[df.len].mul(y, axis=0).sum()
# 切出的数量需要满足需求
for row, num in zip(tmp, df.nums.values):
model.Add(row >= num)
# 限制用料必须是最小根数
model.Add(y.sum() == count)
# 切割结果的总长越高越好
model.Maximize((tmp*df.len.values).sum())
solver = cp_model.CpSolver()
status = solver.Solve(model)
print(solver.StatusName(status))
data["根数"] = y.apply(solver.Value)
result = data[data.根数 != 0].copy()
count = result.根数.sum()
print("总根数：", count, "，采用的切割方案数：", result.shape[0])
display(result)
yl_sum = (result.每根剩余*result.根数).sum()
usage = 1-yl_sum/(count*L)
total_need = (df.nums*df.len).sum()
total = L*result.根数.sum()
print(f"原材料长度：{L}，刀缝：{knife_seam}")
print(f"剩余料头：{yl_sum}，需求利用率：{total_need/total:.2%}，结果利用率：{usage:.2%}")
need = df.set_index("len").nums.convert_dtypes()
need.name = "需求"
report = result.iloc[:, :-2].mul(result.根数, axis=0).sum()
report.name = "切割结果"
report = pd.concat([need, report], axis=1).astype("int")
report["多出根数"] = report.切割结果-report.需求
report 最终在仅采用3种切割方案的情况获得了更好的结果，对于多出的根数，我们可以考虑在某个切割方案的某一根上不进行切割。
也可以将这个过程通过程序自动计算：
from collections import Counter
new_plans = Counter()
main = result.copy()
remain = report.多出根数
while remain.sum() != 0:
idx, max_plan, new_plan = None, None, None
max_num = 0
for i, row in enumerate(main[remain.index].values):
plan = np.c_[row, remain.values].min(axis=1)
if plan.sum() > max_num:
max_num = plan.sum()
max_plan = plan
new_plan = row-plan
idx = i
remain -= max_plan
new_plans[tuple(new_plan)] += 1
main.iloc[idx, -1] -= 1
if main.iloc[idx, -1] == 0:
result.drop(index=main.index[idx], inplace=True)
other = pd.DataFrame(new_plans.keys(), columns=remain.index)
other["每根剩余"] = L-other.mul((df.len+knife_seam).values, axis=1).sum(axis=1)
other["根数"] = new_plans.values()
result = pd.concat([main, other])
yl_sum = (result.每根剩余*result.根数).sum()
usage = 1-yl_sum/(count*L)
print(f"原材料长度：{L}，刀缝：{knife_seam}")
print(f"剩余料头：{yl_sum}，利用率：{usage:.2%}") 最终我们较为完美的解决了钢材切割问题。
基于Gurobi的列生成算法
思路2固然可行，但是假如是将几百个客户的需求合并起来一起优化，长度规格种类数就达到几十种，上述方法可能无法在有效时间内找出有效的可行解。
既然并不是所有的切割模式都会用到，我们可以考虑让程序自动寻找出更好的切割模式，不断迭代直到找到最优解，这样并不需要使用全部的切割方式。
该模型要求在求解模型必须获得对应的对偶变量值，目前市面上已知Gurobi求解器能够良好的支持这一需求，暂时没有发现其他能够获取对偶变量值的求解器。
这里我们直接pip安装该库：
pip install gurobipy 安装后，会自动生成一个免费的pip许可证，我这边的有效期到2023-10-25，意味着可以在Python语言中免费试用1年多。
安装后，只要导包并创建Model对象不报错就证明安装成功。若确保pip认证还在有效期内，导包却报错可能是因为以前安装过gurobi的安装包并认证失败，这时我们可以设置相关的环境变量，再重启python程序（若使用jupyter则需要重启jupyter使环境变量生效）。
具体操作是设置GRB_LICENSE_FILE为pip认证文件的位置：
注意：D:\Miniconda3修改为你的python安装路径。
下面我们先以示例1所使用的数据演示列生成算法，假设需要长度为3米、7米、9米、16米的钢管各25根、30根、14根、8根，目前只有长度为20米的钢管若干，如何切割使消耗钢管的数量最少？
列生成算法在一开始可以随意给出切割方案，然后求解模型获得对应的对偶变量，根据对偶值和切割长度约束求解更优的切割方案。
下面我们的初始方案是每种切割方案都只切一种长度，通过gurobi求解器获取第一个对偶变量值：
from gurobipy import *
# 每根钢管长度，所需的长度，所需的根数
L = 20
lengths = [3, 7, 9, 16]
nums = [25, 30, 14, 8]
# 模型
model = Model()
# 创建变量
z = model.addVars(len(nums), vtype=GRB.CONTINUOUS, name='z')
# 目标函数
model.setObjective(quicksum(z), GRB.MINIMIZE)
# 初始切割方案是每种方案只切一种规格的钢管
model.addConstrs(((z[i]*(L//lengths[i])) >= nums[i]
for i in range(len(nums))), name='mainCon')
# 求解
model.Params.OutputFlag = 0
model.optimize()
# 获取对偶变量值
dual = [con.Pi for con in model.getConstrs()] 这样我们计算出对偶变量值如下：
[0.16666666666666666, 0.5, 0.5, 1.0] 下面我们设置变量c表示一个切割方案中每种规格的钢管的切割根数，然后最小化检验数，若校验数小于0，表示该切割方案更优。校验数即1-dual*c，编码如下：
# 定义子模型
subModel = Model()
c = subModel.addVars(len(nums), vtype=GRB.INTEGER, name='c')
# 添加子问题约束
subModel.addConstr(c.prod(lengths) <= L, name='subCon')
# 目标函数
subModel.setObjective(1-c.prod(dual), GRB.MINIMIZE)
subModel.Params.OutputFlag = 0 # 不输出求解信息
subModel.optimize()
check = subModel.objVal
plan = [int(var.X) for var in subModel.getVars()]
print("校验数：", check, "，切割方案：", plan) 校验数： -0.33333333333333326 ，切割方案： [2, 2, 0, 0] 可以看到校验数小于0，说明该切割方案为更优的切割方案。于是我们将该切割方案添加到主问题中：
| 模式 | 第一种 | 第二种 | 第三种 | 第四种 | 第五种 |
| 3 | 6 | 2 | |||
| 7 | 2 | 2 | |||
| 9 | 2 | ||||
| 16 | 1 |
代码：
column = Column(plan, model.getConstrs())
z_new = model.addVar(vtype=GRB.CONTINUOUS, column=column)
z[z.keys()[-1]+1]=z_new
model.setObjective(quicksum(z), GRB.MINIMIZE)
# 求解更新后的主问题
model.Params.OutputFlag = 0
model.optimize()
dual = [con.Pi for con in model.getConstrs()] 对偶变量值：
[0.0, 0.5, 0.5, 1.0] 然后我们再次求解校验数和当前更优的切割方案：
subModel.setObjective(1-c.prod(dual), GRB.MINIMIZE)
subModel.Params.OutputFlag = 0 # 不输出求解信息
subModel.optimize()
check = subModel.objVal
plan = [int(var.X) for var in subModel.getVars()]
print("校验数：", check, "，切割方案：", plan) 校验数： 0.0 ，切割方案： [0, 0, 2, 0] 此时检验数已经等于0，说明未找到更优的切割方案，那么前面的结果就是最佳结果。当然为了防止前面求解的结果是小数，我们可以将变量类型修改为整数后，重新求解一遍：
for var in model.getVars():
var.setAttr('vtype', GRB.INTEGER)
model.setObjective(quicksum(z), GRB.MINIMIZE)
# 求解更新后的主问题
model.Params.OutputFlag = 0
model.optimize() 整理一下结果：
import pandas as pd
# 输出
print('最少使用钢材数为：', int(model.objVal))
plans = []
for i, length in enumerate(lengths):
plan = [0]*len(nums)
plan[i] = L//length
plans.append(plan)
plans.append([2, 2, 0, 0])
result = pd.DataFrame(plans, columns=lengths)
result["每根剩余"] = L-(result*lengths).sum(axis=1)
result["根数"] = [int(var.X) for var in model.getVars()]
result = result.query("根数!=0")
yl_sum = (result.每根剩余*result.根数).sum()
usage = 1-yl_sum/(result.根数.sum()*L)
print(f"原材料长度：{L}")
print(f"剩余料头：{yl_sum}，利用率：{usage:.2%}") 下面我们将整个过程的代码整理一下：
import pandas as pd
from gurobipy import *
# 每根钢管长度，所需的长度，所需的根数
L, knife_seam = 6000, 3
lengths = [2500, 1200, 1000, 500, 400]
nums = [15, 35, 35, 25, 25]
# 定义主模型和子模型
model = Model()
subModel = Model()
# 创建变量每种切割方案所采用的根数
z = model.addVars(len(nums), obj=1, vtype=GRB.CONTINUOUS, name='z')
# 设置obj后，会有默认的目标函数
# 初始切割方案是每种方案只切一种规格的钢管
plans = []
for i, length in enumerate(lengths):
plan = [0]*len(nums)
plan[i] = L//(length+knife_seam)
plans.append(plan)
model.addConstr((z[i]*(L//(length+knife_seam))) >= nums[i], name='mainCon')
# 求解,OutputFlag = 0表示不输出日志
model.Params.OutputFlag = 0
model.optimize()
# 获取对偶变量值
dual = [con.Pi for con in model.getConstrs()]
# 每种规格钢管的切割根数
c = subModel.addVars(len(nums), obj=dual, vtype=GRB.INTEGER, name='c')
# 添加子问题约束
subModel.addConstr(
c.prod([length+knife_seam for length in lengths]) <= L, name='subCon')
# 修改目标函数为取最大值
subModel.setAttr(GRB.Attr.ModelSense, GRB.MAXIMIZE)
subModel.Params.OutputFlag = 0 # 不输出求解信息
subModel.optimize()
while subModel.objVal > 1:
plan = [int(var.X) for var in subModel.getVars()]
plans.append(plan)
print("校验数：", 1-subModel.objVal, "，切割方案：", plan)
column = Column(plan, model.getConstrs())
z_new = model.addVar(obj=1, vtype=GRB.CONTINUOUS, column=column)
z[z.keys()[-1]+1] = z_new
# 求解更新后的主问题
model.Params.OutputFlag = 0
model.optimize()
# 更新对偶值
dual = [con.Pi for con in model.getConstrs()]
for i in range(len(nums)):
c[i].obj = dual[i]
subModel.Params.OutputFlag = 0 # 不输出求解信息
subModel.optimize()
print("校验数：", check)
for var in model.getVars():
var.setAttr('vtype', GRB.INTEGER)
# 求解更新后的主问题
model.Params.OutputFlag = 0
model.optimize()
# 输出
print('最少使用钢材数为：', int(model.objVal))
result = pd.DataFrame(plans, columns=lengths)
result["每根剩余"] = L - \
(result*[length+knife_seam for length in lengths]).sum(axis=1)
result["根数"] = [int(var.X) for var in model.getVars()]
result = result.query("根数!=0")
yl_sum = (result.每根剩余*result.根数).sum()
count = result.根数.sum()
usage = 1-yl_sum/(count*L)
print(f"原材料长度：{L}，刀缝：{knife_seam}")
print(f"剩余料头：{yl_sum}，利用率：{usage:.2%}")
print("总根数：", count, "，采用的切割方案数：", result.shape[0])
display(result)
report = pd.DataFrame({"需求": nums}, index=lengths)
report["切割结果"] = result.iloc[:, :-2].mul(result.根数, axis=0).sum()
report["多出根数"] = report.切割结果-report.需求
display(report) 对于该结果我们也可以自动拆分多出根数到新的方案：
from collections import Counter
new_plans = Counter()
main = result.copy()
remain = report.多出根数
while remain.sum() != 0:
idx, max_plan, new_plan = None, None, None
max_num = 0
for i, row in enumerate(main[remain.index].values):
plan = np.c_[row, remain.values].min(axis=1)
if plan.sum() > max_num:
max_num = plan.sum()
max_plan = plan
new_plan = row-plan
idx = i
remain -= max_plan
new_plans[tuple(new_plan)] += 1
main.iloc[idx, -1] -= 1
if main.iloc[idx, -1] == 0:
result.drop(index=main.index[idx], inplace=True)
other = pd.DataFrame(new_plans.keys(), columns=remain.index)
other["每根剩余"] = L-(other*[length+knife_seam for length in lengths]).sum(axis=1)
other["根数"] = new_plans.values()
result = pd.concat([main, other])
yl_sum = (result.每根剩余*result.根数).sum()
usage = 1-yl_sum/(count*L)
print(f"原材料长度：{L}，刀缝：{knife_seam}")
print(f"剩余料头：{yl_sum}，利用率：{usage:.2%}")
print("总根数：", count, "，采用的切割方案数：", result.shape[0])