基于整数规划的方形件排样和组批优化问题研究
常见的板式产品如玻璃，PCB板，铝合金门窗等产品因其结构定制化程度高的特点，相关生产制造的企业往往采用“多品种小批量”的个性化生产模式。通过对客户订单的组批，预先规划好各个产品的排样方式，由于所有的产品均由原板材切割生成，因此这种生产模式直接引起生产效率同生产成本之间的矛盾，若能协调好生产订单组批与排样优化难题可以为企业计划排产的顺利进行提供有效的决策支持。针对此类问题，本文从组批与排样过程中的实际约束入手，建立了整数规划模型，合理配置原材料资源和订单组批，从时间和空间两个维度提高原材料利用率和生产效率。
针对问题一，首先在满足各个子集生产订单需求和三阶段齐头切等约束的条件下建立排样整数规划模型，在尽可能少的母板上切割出要求的产品。基于所建模型，设计了一种三阶段启发式搜索树算法(3-Stage Heuristic Search Tree,3-SHST )对模型进行求解。基于宽度非递增原则，采用贪婪式排列算法对原始数据进行排列组合，在FFDH算法的基础上利用后序遍历查找产品项可能插入的位置节点，并生成切割树，遍历所有树节点找到最优插入位置。最后在所给数据集 dataA1~dataA4 中进行算法验证，最终求得四个数据集上所需原板材的个数分别为89块、89块、89块和87块，板材利用率分别为93.87%、93.12%、94.08%和94.08%，排样算法运行时间分别为0.1064s、0.0917s、0.1289s和0.1112s，所有数据集的排样示意图和相关排样信息表格分别以.jpg和.csv文件格式输出。
针对问题二，首先在问题一所建整数规划模型的基础上，增加对每个批次产品项总数和面积总和的约束及相关决策变量，建立组批整数规划模型，在材料利用率最大化的优化目标下平衡生产效率和材料利用率的关系。设计了一种基于生产条件约束的订单层次聚类算法，以生产条件约束为前提，利用杰卡德距离定义各个订单簇之间材料种类的区别，最后利用凝聚层次聚类算法得到最优的组批结果。针对各个批次的排样问题，采用问题一中的排样优化算法得到最优的组批排样方案。最后在所给数据集 dataB1~dataB5 中进行算法验证，最终求得五个数据集上组批次个数分别为 55个、41个、51 个、44 个和 70个，使用板材个数分别为3726块、2498块、2497块、2565块和4004块，板材利用率分别为79.93%、77.14%、77.44%、79.15%和77.17%，组批算法运行时间分别为33.441s、13.2627s、13.7917s、11.3077s和47.0652s。所有数据集的排样示意图和相关组批信息表格分别以.jpg和.csv文件格式输出。
关键词：二维库存切割；排样优化；组批优化；整数规划
一、问题背景与重述
1.1问题背景
随着数字化、网络化、智能化的工业浪潮席卷全球，“智能制造”已被列为“中国制造2025”的主攻方向。在机械制造、纺织加工、印刷作业等领域经常会遇到这样的问题：从一组大的矩形原料上按需求切割下一定尺寸一定需求量的小矩形，按照最朴素的节约资源降低生产成本的想法，希望在使用原材料最少的前提下完成订单需求，这就是典型的二维切削库存（2DCS，two-dimensional cutting stock）问题[1]，其变体问题为二维装箱（2DBP，two-dimensional bin packing）问题[2]。题干中介绍的三阶段的齐头切模式常被用于作为解决上述问题的限制条件之一，这是因为这种切割方式在材料利用率和切割过程的复杂性之间取得了很好的平衡[3]。整数线性规划（ILP，Integer linear programming）常被用于二阶段和三阶段的齐头切2DCS或2DBP问题[4]，文献[2]将Lodi等人的混合整数规划（MILP，The Mixed-Integer Linear Programming）模型[5]扩展到三阶段模式并用于列生成以解决2DBP问题。同时若在原材料利用率最大的要求上加上生产批次与交货时间约束、切割次数最少约束等附加条件则会使优化目标增多，约束条件也更加复杂，求解也变得更为困难。本赛题正是在此背景下要解决小批次多种类的订单的二维切削库存问题的排样优化和订单组批问题。
批量订单数据
| 订单组批优化 生成订单分批方案 | 组批 |
组批方案
| 订单排样优化 生成订单排样方案 | 排样 |
排样方案
输出排样结果
图 1组批与排样流程
3
1.2问题重述
基于上述研究背景，题目提供了两个独立的数据集A和B，分别包含产品项数量、需求、尺寸、产品id、订单号等信息，本文需要解决以下两个问题：
问题1：排样优化问题。对数据集A中的每个子集，在满足各子集生产订单需求和三阶段齐头切约束条件下，建立相关的模型并设计求解算法，使用最少的相同尺寸的原板材
来完成订单生产需求，使板材利用率达到最大化，并按格式输出排样图和排样信息表格。
此问题中未对订单号作相关要求，且同一子集中产品的材料也相同，因此这两项数据可以不做考虑。这是一个典型的三阶段齐头切2DCS或2DBP问题。
问题2：订单组批问题。对数据集B中的每个子集，在问题一约束的基础上，增加以下进一步约束：
1) 每份订单当且仅当出现在一个批次中；
2) 每个批次中的相同材质的产品项（item）才能使用同一块板材原片进行排样； 3) 为保证加工环节快速流转，每个批次产品项（item）总数不能超过1000； 4) 因工厂产能限制，每个批次产品项（item）的面积总和不能超过250m 2；
要根据上述所有条件，建立规划模型并设计相应的求解算法，对数据集B中的全部订单进行组批，然后对每个批次进行独立排样，使得板材原片的用量最少。本问题虽然最终的优化目标与问题一相同，但由于加上了同批次 item 数量，同批次 item 面积和以及同材质item排样三个条件，实质上是在问题一的基础上上升了一个维度，提出了一个子问题，即如何组批才能平衡生产效率和材料利用率的关系，使得材料利用率最大化。同时数据集B的数据量比A要大的多，这也就对模型和算法的求解效率提出了需求。
1.3技术路线图
问题一：方形件排样优化问题
生成订单分批方案
图 2本文技术路线图
二、基本假设和符号说明
2.1基本假设
基于题干和实际求解中可能考虑到的因素，在求解之前我们需做出以下假设以保证模
型和算法的逻辑完备：
（1）本次赛题所有订单的交货期均相同，不做区分；
（2）只考虑齐头切的切割方式（直线切割、切割方向垂直于板材一条边，并保证每次直线切割板材可分离成两块）；
（3）切割阶段数不超过3，同一个阶段切割方向相同；
（4）排样方式为精确排样；
（5）假定板材原片仅有一种规格且数量充足；
（6）排样方案不用考虑锯缝宽度（即切割的缝隙宽度）影响； （7）数据集A和B中的各子集数据是独立的，应分别处理。
2.2符号说明
针对贯穿全文的一些主要变量及其符号含义列表如下。
表 1符号及变量定义说明
| 符号 | 符号含义 | 单位 | |
| Area i | 数据集Ai中item面积总和 | m² | |
| L | 所使用板材原片的长 | mm | |
| W | 所使用板材原片的宽 | mm | |
| il | item的长 | mm | |
| iw | item的宽 | mm | |
| d | 所有stripe的尺寸类别数 | 个 | |
| t | 能生成stripe的stack的个数 | 个 | |
| m | |||
| 所有stripe的最大顺序编号 | |||
| kq a | 0/1变量，顺序编号为k的stack是否可以被编号为a的stripe包含 | ||
| 个 | |||
| aiy | 0/1变量，顺序编号为i的item是否可以被编号为a的stack包含 | ||
| ax | 正整数变量，在顺序编号为a的stack中包含item的个数 | ||
| az | 正整数变量，在顺序编号为a的stripe中包含stack的个数 | 个 | |
| J A B | 杰卡德相似系数 | ||
| J ( , ) A B | 杰卡德距离 |
三、问题一：方形件排样优化问题求解
3.1数据分析与思考
首先对给出的 dataA1-dataA4 四组数据进行粗略的观察和分析，可分别计算出每组数据中item的面积之和，以dataA1数据为例，item面积之和为：
| Area 1 | = | 752 i = 0 | length i | * | width i | = | 2.4868561455 10 ( | m | 2 | ) | (3.1) |
用面积除以原板材面积可得极限情况下最小使用的板材数：
| Area 1 | | (1220 | | 2440 10 ) | = | 83.54 | | 84 | (3.2) |
以 item 面积之和除以极限情况下使用板材面积之和可得 dataA1 所提供数据的原材料极限利用率：
| Area 1 / (84*2440*1220 10 ) | = | 99.45% | (3.3) |
当然，由于题干中给出的不能拼接等条件以及三阶段齐头切方式所带来的必然的材料
浪费，实际的原材料利用率不肯可能达到以上计算值，同样的，如果模型及算法的计算结
果显示利用率超过此值那说明结果也是错误的。对其他三组数据进行同样的处理得出的结果如表2所示：
表 2极限情况下材料利用率分析
| 数据集 | dataA1 | dataA2 | dataA3 | dataA4 |
| 极限情况下所需 | 84 | 83 | 84 | 82 |
| 原材料块数 | ||||
| 极限利用率 | 99.45% | 98.66% | 99.68% | 97.44% |
继续对数据进行观察，发现每个数据集中的方形件的需求量都是1且很少有长宽完全一样的两个方形件，但是若以长度或者宽度来进行排序，可以发现有不少的方形件具有同样的长度或者宽度，这也符合题干中所说的同一个 stack 中 item 的长或宽需相等的条件。为了更加直观的看出方形件的外观尺寸分布情况，利用数据集中的数据在 origin 软件中画出如图3所示散点频数统计图，X、Y坐标分别是方形件的长和宽，XY平面内是方形件的长宽散点图，在 XZ 和 YZ 平面内我们作出了某一长度或宽度区间段内方形件的频数统计直方图，可以看出在某些区间内方形件的长或宽特征体现出高度集中的特征。
6
（1）在相同栈（stack）里的产品项（item）的宽度（或长度）应该相同； （2）最终切割生成的产品项是完整的，非拼接而成。
优化的目标是要满足尽可能少的板材原片数量，由于本问题为典型的 NP-hard 问题，一般没有polynomal复杂度的算法，算法复杂度随着变量的增加呈指数级增长（指数爆炸），因此在上述题干约束的基础上引申做出以下约束：
（1）切割方式固定为第 1 阶段采用横向切割生成 stripe（条带），第 2 阶段采用纵向切割生成stack（栈），第三阶段采用横向切割生成item（产品项）；
（2）每个stack中的第一个(最上面)元素是stack中最宽的元素；
（3）每个stripe中的第一个(最左边)元素是stripe中最宽的元素；
（4）所有被选择的item的宽度按照不递减排列；
为解释模型需要，对问题一涉及的符号作进一步说明：
（2）假设最终可能会产生d种不同尺寸的stripe，对于每一个stripe，都可以通过他们
（3）最终产生的stack个数为m，每个stack的产生都源于最初插入了一个顺序编号为i的 item，记此 stack 的顺序编号为a；与此类似，假定最终可能会产生m个 stripe，每一个stripe的产生都源于最初插入了一个顺序编号为k的stack，记此stripe的顺序编号为a。这里符号不进行区分只用来表明元素的产生与其组成之间存在上述关系。
基于上述所有约束条件，我们选择通过整数规划对该问题进行建模，模型如下： （一）决策变量：
代码
1. class TreeNode():
2. def __init__(self, length, width, depth = 0, sum = 0, item = None) -> None:
3. self.length = length # 节点长度信息
4. self.width = width # 节点宽度信息
5. self.depth = depth # 节点所处树的位置，root为-1
6. self.sum = sum # 节点长度或宽度和 依据所处位置
7. self.item = item # 只有叶子节点有该属性，表示处于叶子节点处的实际item
8. self.childs = list()
9.
10. def insert_node(self, node):
11. for child in self.childs:
12. if child.insert_node(node):
13. return True
14. if self.depth == 3:
15. return False
# 当节点深度为2的时候，子节点的切法是横着切，只在当前节点下方加一个新节点
18. elif self.depth == 2:
19. cur_pos = self.childs[len(self.childs) - 1].sum + node.width
20. if cur_pos <= self.width and node.length == self.length: # 查看是否有空间 放置(暂时不考虑旋转)
21. node.sum = cur_pos
22. node.depth = 3
23. self.add_child(node)
24. return True
25. # 当节点深度为1的时候， 子节点的切法是竖着切，在当前节点下方要加两个新节点(一个 stack和一个item)
26. elif self.depth == 1:
27. cur_pos = self.childs[len(self.childs) - 1].sum + node.length
28. if node.width <= self.width and cur_pos <= self.length:
29. stack = TreeNode(node.length, self.width, depth = 2 ,sum = cur_pos)
30. self.add_child(stack)
31. node.sum = node.width
32. node.depth = 3
33. stack.add_child(node)
34. return True
35. # 当处理节点深度为0，即该块是某个原件，子节点的切法是横着切，在当前节点下方要加入三个 新节点(一个strip，一个stack和一个item)
36. elif self.depth == 0:
37. cur_pos = self.childs[len(self.childs) - 1].sum + node.width
38. if cur_pos <= self.width:
39. strip = TreeNode(self.length, node.width, depth = 1, sum = cur_pos)
40. self.add_child(strip)
41. stack = TreeNode(node.length, strip.width, depth = 2, sum = node.lengt
h)
42. strip.add_child(stack)
43. node.sum = node.width
44. node.depth = 3
45. stack.add_child(node)
46. return True
47. # 当处理节点深度为-1时，即当前无原件可以放置当前块，则要拿一块新原件放置该块，在根节点下方要加入四个新节点(一个sheet， 一个strip，一个stack和一个item)
48. else:
49. sheet = TreeNode(self.length, self.width, depth = 0)
50. self.add_child(sheet)
51. strip = TreeNode(sheet.length, node.width, depth = 1, sum = node.width)
52. sheet.add_child(strip)
53. stack = TreeNode(node.length, strip.width, depth = 2, sum = node.length)
54. strip.add_child(stack)
5. node.sum = node.width
56. node.depth = 3
57. stack.add_child(node)
58. return True
59.
60. return False
61.
62. def add_child(self, node):
63. self.childs.append(node)
64. # sorted(self.childs)
65.
66. def __str__(self):
67. return "length: {}, width: {}, depth: {}".format(self.length, self.width, self
.depth
可视化排样方案程序：
1. from PIL import Image, ImageDraw
2. from cut_tree_node import Item, TreeNode
3.
4. def visualize_sheet(node, save_file, fill = (135, 206, 235), outline = "black", width =
2):
5. # 可视化图，返回方形件被占面积以及方形件中所有被切割元件信息
6. # fill，outline，width为矩形的填充，边框等信息
7. # 深度优先搜索
8. def dfs_node(node):
9. info = {}
10. node_index = 0
11. if node.depth == 2:
12. node_index = node.sum - node.length
13. elif node.depth == 1 or node.depth == 3:
14. node_index = node.sum - node.width
15. info[node_index] = list()
16. for child in node.childs:
17. info[node_index].append(dfs_node(child))
18.
19. if node.depth == 3:
20. return {"y": node_index, "item": node.item}
21. return info
22.
23. img = Image.new('RGB', (node.length, node.width), (255, 255, 255))
24. draw = ImageDraw.Draw(img)
25. info = dfs_node(node)[0] # 通过深度优先搜索遍历树来得到被切割块所处x, y坐标
26.
27. # 返回值
28. sum_area = 0 # 当前方形件被占面积
