📢作者： 小小明-代碼實體
倉庫選址問題
倉庫選址是運輸行業最常見的問題，目标是使整個配送系統的總成本最低。
例如有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])