组合优化问题
组合优化（Combinatorial Optimization，CO）数学优化研究的一个分支。主要关注的是从有限的对象集合中寻找最优解的问题。这个词的由来主要是由“组合”和“优化”两部分构成。“组合”指的是从有限的对象集合中选择一部分的过程，而“优化”则指的是在满足一定条件下，使得某个目标函数达到最大或最小。即：组合优化问题的解集合有限个的，要做的事情是在有限的集合里面找到目标最优的集合。
组合优化问题在许多实际情况中都有出现，包括经济管理、工业制造、交通运输等领域。比如：
- 0-1 背包问题：如何把体积有限、价值有限，装入容积有限的背包里，获得最大价值的装包？
- 装箱问题，bin packing，如何装箱使得所用的箱子最少?
- TSP的旅行商问题：如何选择一条道路遍历每个城市，路径最短？
- 车间调度问题：车间的加工顺序如何安排完工时间最小？
- 地图块着色问题：如何采用少量的颜色给地图块着色，两两相邻颜色不一样？
组合优化问题特征：
- 解可数的，可行域是有限的
- 最优解一定存在
- 可行方案非常多，无法用枚举的方案，NP-hard的问题
解组合优化的方式有：
- 数学规划的方法，利用数学规划求解器，寻找精确最优解。缺点：耗时。
- 启发式算法：近似方法、贪婪算法、基于规则等。解的质量有基础保证，但是耗时不确定。
- 元启发式算法：遗传算法、蚁群算法、禁忌搜索、模拟退火等。不同的问题需要调参数，解的质量不可靠。
- 机器学习算法：监督学习、强化学习、图神经网络。此方法还在探索的算法领域，对于同分布的数据，机器学习后能进行快速的求解。
本文将以一个简化的装箱问题为例，来讲解采用数学规划的方法来解决装箱这个组合优化问题。
什么是装箱问题？
装箱问题（Bin Packing Problem）是组合优化领域中的一个经典问题，主要涉及如何将一系列对象高效地装入有限数量的容器（或“箱”）中，同时满足特定的约束条件。这个问题的目标是最小化所需使用的箱子数量或者最大化箱子的装载效率，以减少空间或资源的浪费。
在数学领域，装箱问题通常被建模为整数规划问题。它的变体众多，包括一维装箱问题、二维装箱问题、多维装箱问题等，以及它们各自的约束条件，如容器容量、物品形状、摆放方向限制等。
装箱问题可以用在什么场景？
装箱问题广泛应用于各种实际场景，几个典型的应用示例包括：
- 物流与配送: 确定最小数量的货车来运输一定量的货物，同时考虑每辆货车的载重限制、体积限制以及货物的装载平衡。
- 仓库存储: 如何利用有限的仓储空间存储尽可能多的库存货物，或如何将货物高效地摆放到货架上，以便节省空间并方便存取。
- 生产包装: 在产品的包装过程中，决定最少的包装盒数量来装载不同大小的产品，或者在确保产品安全的前提下，使得各个包装盒中的空隙最小化。
- 计算资源分配: 在计算机科学领域，将不同的任务分配给有限数量的处理器或内存资源，以求最大化资源的使用效率。
- 切割优化和材料利用: 在制造业中，比如钢材、木材、玻璃等材料的切割，如何从原材料中切割出所需部件，以减少材料浪费。
- 行李装载: 确定如何在飞机的货舱内或个人行李箱中高效地安置行李和货物，以保证安全、均衡且有效地利用空间。
- 云计算资源优化: 决定如何在云计算环境中分配虚拟机到物理服务器上，以最小化服务器的数量和能耗。
装箱问题
假设我们有一批货物要运输走，这批货物包装盒的比较整齐，侧面面尺寸一样，仅厚度不同、重量不同。运输货车里可以放两列，如上图示意。
货车的总承重是10吨。为了行驶平稳，希望左右摆的货物质量差异不超过500公斤。货车可装货物的总长度（厚度）是5米。
要装的货物的信息如此表格文件所示： data/items-1010.csv
部分数据如下：
| 序号 | 包装和厚(cm) | 重量(kg) |
|---|---|---|
| 1 | 101 | 505 |
| 2 | 61 | 305 |
| 3 | 126 | 630 |
| …… | …… | …… |
| 100 | 85 | 765 |
| …… | …… | …… |
| 200 | 60 | 600 |
| …… | …… | …… |
| 1010 | 80 | 640 |
当数据维度太高的时候，采用数学规划方式计算会很慢。
我们可以拆解问题，比如分块去算某100个货物的装载方案，切割前100条数据为 data/items-100.csv 。更多分块数据方案用户可以进data文件夹自己修改。
装箱问题的数学规划建模
下面我们采用数学规划的方法来建模这个问题。
首先：定义几何，描述数据，方便后面引用
- 货物集合
I
I
I
- 它的描述信息参数是 厚度 l i l_{i} li 和 重量 w i w_{i} wi
- 假设有车辆集合 B B B
- 问题中有左侧和右侧的问题，定义集合 s = l e f t , r i g h t s = { left, right } s=left,right
- 关于左右侧的差异在500kg以内，在表达式上，会需要一个绝对值，属于非线性。变更成线性更简单。这里我们假设所有的火车上都是左边的重，右边的轻。则左边重量上限 W l e f t m a x = 5250 W_{leftmax} = 5250 Wleftmax=5250 kg ，右边重量上限 W r i g h t m a x = 4750 W_{rightmax} = 4750 Wrightmax=4750 kg。左边比右边重的差异上限是 W l e g t r i g h t = 500 W_{legt_right} = 500 Wlegtright=500 kg
- 两边的装货厚度要求是参数 L m a x = 500 L_{max} = 500 Lmax=500 cm
然后，设置变量来表达要解决的问题。
- 货物如何摆放，设置 0-1 变量 x i , b , s x_{i,b,s} xi,b,s，其中 i ∈ I i \in I i∈I, b ∈ B b \in B b∈B, s ∈ S s \in S s∈S。 当变量取值为1的时候，则代表货物摆放在该位置，为0的时候则不放置。
由此，我们可以用来表达装货的约束。
- 每个货物只能放一个位置,且至少是一个位置： ∑ b ∈ B , s ∈ S x i , b , s = 1 , ∀ i ∈ I \sum_{b \in B, s \in S} x_{i,b,s} = 1, \forall i \in I ∑b∈B,s∈Sxi,b,s=1,∀i∈I
- 每个货车上，装载的货物重量限制：
- 左边限制： ∑ i ∈ I x i , b , l e f t ⋅ W i ≤ W l e f t m a x , ∀ b ∈ B \sum_{i \in I} x_{i,b,left} \cdot W_{i} \leq W_{leftmax} , \forall b \in B ∑i∈Ixi,b,left⋅Wi≤Wleftmax,∀b∈B
- 右边限制： ∑ i ∈ I x i , b , r i g h t ⋅ W i ≤ W r i g h t m a x , ∀ b ∈ B \sum_{i \in I} x_{i,b,right} \cdot W_{i} \leq W_{rightmax} , \forall b \in B ∑i∈Ixi,b,right⋅Wi≤Wrightmax,∀b∈B
- 左右差异限制： 0 ≤ ( ∑ i ∈ I x i , b , l e f t ⋅ W i ) − ( ∑ i ∈ I x i , b , r i g h t ⋅ W i ) ≤ W l e g t r i g h t , ∀ b ∈ B 0 \leq (\sum_{i \in I} x_{i,b,left} \cdot W_{i}) - (\sum_{i \in I} x_{i,b,right} \cdot W_{i}) \leq W_{legt_right} , \forall b \in B 0≤(∑i∈Ixi,b,left⋅Wi)−(∑i∈Ixi,b,right⋅Wi)≤Wlegtright,∀b∈B
- 每辆货车的每一边的装货长度不超过5米
∑ b ∈ B x i , b , s ⋅ L i ≤ L m a x , ∀ b ∈ B , s ∈ S \sum_{b \in B} x_{i,b,s} \cdot L_{i} \leq L_{max}, \forall {b \in B, s \in S} ∑b∈Bxi,b,s⋅Li≤Lmax,∀b∈B,s∈S
此外，我们需要设置目标，用的车辆数最少。
我们可以设置新的变量，增加约束来表述和上面的装载方案 x i , b , s x_{i,b,s} xi,b,s的关系。
- 设置新变量车辆是否有使用 y b ∈ { 0 , 1 } y_{b} \in \{0,1\} yb∈{0,1}
- 由此，目标函数表达为： m i n i m i z e ∑ b ∈ B y b minimize \sum_{b \in B} y_{b} minimize∑b∈Byb
- 增加约束：
- 每个车辆上货物的摆放变量加和大于等于车辆被使用。
- 即如果没有放货物，都是0，等于关系。如果有放1个货物，都是1，等于关系。如果放两个及以上货物，x的加和大于y。
- 数学公式表达为： ∑ i ∈ I , s ∈ S x i , b , s ≥ y b , ∀ b ∈ B \sum_{ i \in I, s \in S} x_{i,b,s} \geq y_{b}, \forall b \in B ∑i∈I,s∈Sxi,b,s≥yb,∀b∈B
- 再构造y大于x的限制。有车辆被使用的y大于等于每一个x，即如果车辆被使用，y是1，和这个车有关系的x是1或者0，
- 数学公式表达为： y b ≥ x i , b , s , ∀ i ∈ I , b ∈ B , s ∈ S y_{b} \geq x_{i,b,s}, \forall i \in I, b \in B, s\in S yb≥xi,b,s,∀i∈I,b∈B,s∈S
- 每个车辆上货物的摆放变量加和大于等于车辆被使用。
汇总得到数学公式如下：
min ∑ b ∈ B y b s.t. ∑ i ∈ I x i , b , l e f t ⋅ W i ≤ W l e f t m a x , ∀ b ∈ B ∑ i ∈ I x i , b , r i g h t ⋅ W i ≤ W r i g h t m a x , ∀ b ∈ B 0 ≤ ( ∑ i ∈ I x i , b , l e f t ⋅ W i ) − ( ∑ i ∈ I x i , b , r i g h t ⋅ W i ) ≤ W l e g t r i g h t , ∀ b ∈ B ∑ b ∈ B x i , b , s ⋅ L i ≤ L m a x , ∀ b ∈ B , s ∈ S ∑ i ∈ I , s ∈ S x i , b , s ≥ y b , ∀ b ∈ B y b ≥ x i , b , s , ∀ i ∈ I , b ∈ B , s ∈ S x i , b , s ∈ { 0 , 1 } y b ∈ { 0 , 1 } \begin{array}{rll} \min & \sum_{b \in B} y_{b} & \\ \text{s.t.} & \sum_{i \in I} x_{i,b,left} \cdot W_{i} \leq W_{leftmax} , & \forall b \in B \\ & \sum_{i \in I} x_{i,b,right} \cdot W_{i} \leq W_{rightmax} , & \forall b \in B \\ & 0 \leq (\sum_{i \in I} x_{i,b,left} \cdot W_{i}) - (\sum_{i \in I} x_{i,b,right} \cdot W_{i}) \leq W_{legt_right} , & \forall b \in B \\ & \sum_{b \in B} x_{i,b,s} \cdot L_{i} \leq L_{max}, & \forall b \in B, s \in S \\ & \sum_{ i \in I, s \in S} x_{i,b,s} \geq y_{b}, & \forall b \in B \\ & y_{b} \geq x_{i,b,s}, & \forall i \in I, b \in B, s\in S \\ & x_{i,b,s} \in\{0,1\} & \\ & y_{b} \in \{0,1\} & \\ \end{array} mins.t.∑b∈Byb∑i∈Ixi,b,left⋅Wi≤Wleftmax,∑i∈Ixi,b,right⋅Wi≤Wrightmax,0≤(∑i∈Ixi,b,left⋅Wi)−(∑i∈Ixi,b,right⋅Wi)≤Wlegtright,∑b∈Bxi,b,s⋅Li≤Lmax,∑i∈I,s∈Sxi,b,s≥yb,yb≥xi,b,s,xi,b,s∈{0,1}yb∈{0,1}∀b∈B∀b∈B∀b∈B∀b∈B,s∈S∀b∈B∀i∈I,b∈B,s∈S
采用建模语言 MindOpt APL 编程和求解
下面的内容包含了code源码，如果您当前在浏览项目内容页，请复制项目后，点击右上角的NoteBook按钮，进入NoteBook环境中查看和运行代码。
In [1]:
clear model;
option modelname binpacking_01; #定义存储文件名
# ----------建模--------Start----
# binpacking_01.mapl
## 导入数据-------
print "导入数据-------";
param fileDir = "data/";
set Items = {read fileDir+"items.csv" as "<1n>" skip 1}; #读取货物序号
print "总共有{}个货物待装载"% card(Items);
param ItemSize[Items] = read fileDir+"items.csv" as "<1n> 2n" skip 1; #读取货物的长度（厚度），单位厘米cm
param ItemWeight[Items] = read fileDir+"items.csv" as "<1n> 3n" skip 1; #读取货物的重量，单位千克kg
param BinLength = 500; #货车的长度,单位厘米cm
param BinLoadCapacity = 10*1000; #货车的载重，单位kg
param BinLeftRightBalanceDiff = 500; #或者左右两边的重量差异千克
# 因为左右无定数，假设货物都是重的放左边，轻的右边，则左边载重不超过 5250，右边不超过4750.
param BinLoadCapacity_Left = BinLoadCapacity * 0.5 + 500 * 0.5;
param BinLoadCapacity_right = BinLoadCapacity * 0.5 - 500 * 0.5;
# 假设总共有x辆货车
param Bin_Num = 20;
set Bins = {1..Bin_Num};
print "先假设有{}辆货车"%Bin_Num;
## 左右两边
set Sides = {"left","right"};
## 数学建模-------
print "数学建模--------------";
# 变量----
var x_Item2Bin[Items*Bins*Sides] >= 0 binary; #货物放在哪个车的哪一边
var x_BinUsed[Bins] >= 0 binary;; #总共被使用的货车数
# 目标----
# 目标函数：最小化使用货车数量
minimize TotalBins: sum { b in Bins} x_BinUsed[b];
# 货车的总承重是10吨。为了行驶平稳，希望左右摆的货物质量差异不超过500公斤。货车可装货物的总长度（厚度）是5米。 要装的货物的信息如下表格所示：
# 约束----
# 每个货物只能放一个位置,且至少是一个位置
subto constraint_0:
forall {i in Items}
sum {b in Bins, s in Sides} x_Item2Bin[i,b,s] == 1;
# 每辆货车总载重上限是10吨 ， 左右两边差异有限制500kg
# 假设左边的重，变换上限和差异，定义左边比右边重
subto constraint_1:
forall {b in Bins}
sum {i in Items} x_Item2Bin[i,b,"left"] * ItemWeight[i] <= BinLoadCapacity_Left;
subto constraint_1_2:
forall {b in Bins}
sum {i in Items} x_Item2Bin[i,b,"right"] * ItemWeight[i] <= BinLoadCapacity_right;
subto constraint_1_3:
forall {b in Bins}
0 <= (sum {i in Items} x_Item2Bin[i,b,"left"]* ItemWeight[i]
- sum {i in Items} x_Item2Bin[i,b,"right"] * ItemWeight[i]) <= BinLeftRightBalanceDiff;
# 每辆货车的每一边的装货长度不超过5米
subto constraint_2:
forall {b in Bins}
forall { s in Sides}
sum { i in Items} x_Item2Bin[i,b,s] * ItemSize[i] <= BinLength;
# 车被用了
subto constraint_3:
forall {b in Bins}
sum {i in Items, s in Sides} x_Item2Bin[i,b,s] >= x_BinUsed[b];
subto constraint_4:
forall {b in Bins}
forall { s in Sides}
forall { i in Items}
x_BinUsed[b] >= x_Item2Bin[i,b,s];
#求解
option solver mindopt; # （可选）指定求解用的求解器，默认是MindOpt
#option mindopt_options 'max_time=600'; #设置求解器参数
option mindopt_options 'print=0'; #设置求解器输出级别，减少过程打印
solve; # 求解
# 结果打印和检查结果
print "----------------结果打印和存文件--------------";
#display;
print "最小需要货车数 = ",sum { b in Bins} x_BinUsed[b];
print "------每辆车----";
print "|车ID|左边货物数|右边货物数|……|左边厚度|右边厚度|……|左边重量|右边重量|左右重量差异|";
print "|--|--|--|--|--|--|--|--|--|--|";
forall { in Bins with x_BinUsed[b] >= 0.5}
print "|{}|{}|{}|……|{}|{}|……|{}|{}|{}|"
%b,
sum{i in Items} x_Item2Bin[i,b,"left"],
sum{i in Items} x_Item2Bin[i,b,"right"],
sum{i in Items} x_Item2Bin[i,b,"left"] * ItemSize[i],
sum{i in Items} x_Item2Bin[i,b,"right"]* ItemSize[i],
sum{i in Items} x_Item2Bin[i,b,"left"] * ItemWeight[i],
sum{i in Items} x_Item2Bin[i,b,"right"]* ItemWeight[i],
sum{i in Items} x_Item2Bin[i,b,"left"] * ItemWeight[i] - sum{i in Items} x_Item2Bin[i,b,"right"]* ItemWeight[i] ;
#打印太长，注释。可以 cmd+/ 来批量注释和解除注释，或者加“and i <=5"来只打印前几行
print "------每个货物----";
print "|货物ID|包装厚度(cm）|重量（kg)|放置车辆|边|";
print "|--|--|--|--|--|";
forall { in Items*Bins*Sides with x_Item2Bin[i,b,s] >= 0.5 and i <=5}
print "|{}|{}|{}|{}|{}|"
%i,ItemSize[i],ItemWeight[i],b,s;
print "|……|……|……|……|……|";
# 存文件
print "车ID,左边货物数,右边货物数,……,左边厚度,右边厚度,……,左边重量,右边重量,左右重量差异" : "车辆装载情况.csv";
close "车辆装载情况.csv";
forall { in Bins with x_BinUsed[b] >= 0.5}
print "{},{},{},……,{},{},……,{},{},{}"
%b,
sum{i in Items} x_Item2Bin[i,b,"left"],
sum{i in Items} x_Item2Bin[i,b,"right"],
sum{i in Items} x_Item2Bin[i,b,"left"] * ItemSize[i],
sum{i in Items} x_Item2Bin[i,b,"right"]* ItemSize[i],
sum{i in Items} x_Item2Bin[i,b,"left"] * ItemWeight[i],
sum{i in Items} x_Item2Bin[i,b,"right"]* ItemWeight[i],
sum{i in Items} x_Item2Bin[i,b,"left"] * ItemWeight[i] - sum{i in Items} x_Item2Bin[i,b,"right"]* ItemWeight[i]
>> "车辆装载情况.csv";
close "车辆装载情况.csv";
print "货物ID,包装厚度(cm),重量(kg),放置车辆,放置边" : "货物装载方案.csv";
close "货物装载方案.csv";
forall { in Items*Bins*Sides with x_Item2Bin[i,b,s] >= 0.5}
print "{},{},{},{},{}"
%i,ItemSize[i],ItemWeight[i],b,s
>> "货物装载方案.csv";
close "货物装载方案.csv";
print "结果文件存储在: [车辆装载情况.csv](车辆装载情况.csv) 和 [货物装载方案.csv](货物装载方案.csv)";
运行上述代码，结果如下：
导入数据-------
总共有100个货物待装载
先假设有20辆货车
数学建模--------------
Running mindoptampl
wantsol=1
print=0
MindOpt Version 1.1.0 (Build date: 20240123)
Copyright (c) 2020-2024 Alibaba Cloud.
Start license validation (current time : 01-MAR-2024 11:22:15).
License validation terminated. Time : 0.011s
OPTIMAL; objective 9.00
Completed.
----------------结果打印和存文件--------------
最小需要货车数 = 9
------每辆车----
|车ID|左边货物数|右边货物数|……|左边厚度|右边厚度|……|左边重量|右边重量|左右重量差异|
|--|--|--|--|--|--|--|--|--|--|
|1|5|5|……|500|462|……|4363|4237|126|
|3|7|10|……|488|411|……|3393|3046|347|
|6|6|5|……|499|474|……|4526|4053|473|
|8|5|5|……|482|499|……|3574|3112|462|
|9|6|4|……|491|488|……|4454|4412|42|
|10|6|5|……|461|488|……|4118|3812|306|
|11|6|6|……|487|440|……|3931|3816|115|
|13|4|5|……|492|469|……|4460|4135|325|
|14|5|5|……|414|492|……|3948|3544|404|
------每个货物----
|货物ID|包装厚度(cm）|重量（kg)|放置车辆|边|
|--|--|--|--|--|
|1|101|505|14|left|
|2|61|305|3|right|
|3|126|630|8|right|
|4|32|256|3|left|
|5|39|429|13|right|
|……|……|……|……|……|
结果文件存储在: [车辆装载情况.csv](车辆装载情况.csv) 和 [货物装载方案.csv](货物装载方案.csv)
求解结果
上面的打印的结果部分为：
总共有100个货物待装载, 先假设有20辆货车
最小需要货车数 = 9
------每辆车----
| 车ID | 左边货物数 | 右边货物数 | …… | 左边厚度 | 右边厚度 | …… | 左边重量 | 右边重量 | 左右重量差异 |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 5 | 5 | …… | 500 | 462 | …… | 4363 | 4237 | 126 |
| 3 | 7 | 10 | …… | 488 | 411 | …… | 3393 | 3046 | 347 |
| 6 | 6 | 5 | …… | 499 | 474 | …… | 4526 | 4053 | 473 |
| 8 | 5 | 5 | …… | 482 | 499 | …… | 3574 | 3112 | 462 |
| 9 | 6 | 4 | …… | 491 | 488 | …… | 4454 | 4412 | 42 |
| 10 | 6 | 5 | …… | 461 | 488 | …… | 4118 | 3812 | 306 |
| 11 | 6 | 6 | …… | 487 | 440 | …… | 3931 | 3816 | 115 |
| 13 | 4 | 5 | …… | 492 | 469 | …… | 4460 | 4135 | 325 |
| 14 | 5 | 5 | …… | 414 | 492 | …… | 3948 | 3544 | 404 |
------每个货物----
| 货物ID | 包装厚度(cm） | 重量（kg) | 放置车辆 | 边 |
|---|---|---|---|---|
| 1 | 101 | 505 | 14 | left |
| 2 | 61 | 305 | 3 | right |
| 3 | 126 | 630 | 8 | right |
| 4 | 32 | 256 | 3 | left |
| 5 | 39 | 429 | 13 | right |
| …… | …… | …… | …… | …… |
另外，上面的结果还存在了文件 车辆装载情况.csv 和 货物装载方案.csv，可在下文的案例地址中查看。
高阶方案尝试
读者可尝试复制本项目，增加难度进行调试。比如：
- 修改data文件夹的item.csv数据，如
data/items-100-2.csv，观察不同数据的结果区别。 - 增加分块的量去测试对速度的影响，如
data/items-200.csv的200个货一起算装载方式。 - 将1010个货物的所有的装载方案计算完。
- 将货车分成不同型号和载重，更改模型进行测试。
- 改进上面的数学建模模型，使得计算更简单快速。
可在云上建模平台进行调试（免安装）,案例地址
