📢作者： 小小明-代码实体
📢博客主页：https://blog.csdn.net/as604049322
📢欢迎点赞 👍 收藏 ⭐留言 📝 欢迎讨论！
本文链接：https://blog.csdn.net/as604049322/article/details/126324174
仓库选址问题
仓库选址是运输行业最常见的问题，目标是使整个配送系统的总成本最低。
例如有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))
cost
每个选址点的运营成本和最大容量：
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)
data
可以看到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%}")
result
最终我们较为完美的解决了钢材切割问题。
基于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()]
dual
这样我们计算出对偶变量值如下：
[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，说明该切割方案为更优的切割方案。于是我们将该切割方案添加到主问题中：
模式第一种第二种第三种第四种第五种36272292161
代码：
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()]
dual
对偶变量值：
[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%}")
result
下面我们将整个过程的代码整理一下：
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')
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')
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])
result
参考：https://blog.csdn.net/weixin_38442390/article/details/121864422
Original: https://blog.csdn.net/as604049322/article/details/126324174
Author: 小小明-代码实体
Title: OR-Tools求解仓库选址和钢材取料问题
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/694983/
转载文章受原作者版权保护。转载请注明原作者出处！