堆集装箱翻箱问题的
整数规划模型
（BRP-Ⅰ、BRP-Ⅱ及代码）
系列一
前 言
因为现代供应链系统受到时间的限制，所以如何快速按需访问库存也是当前社会研究的热点之一。在实际生活中，如码头、仓库等，我们都能经常看见箱子或其他物品的堆叠存储和库存访问。
-
今天小编就将为大家介绍在类似场景下，如何移动箱子，以实现在按顺序获取库存的前提下使得箱子的移动次数最少的模型。
对这类问题的研究其实已经开展得十分广泛了，所以小编计划会推出一个集装箱翻箱问题的整数规划模型系列，分别介绍不同文献中的多种模型。
本期目录
01 | 问题描述 |
|---|---|
02 | 建立模型 |
03 | 代码 |
04 | 结果展示 |
问题描述
术语：
1、block：堆叠在某一区域的同质化的物品，可以简单地看成一个个集装箱。
2、slot：在这一区域中，block可以放置的位置。
3、relocation：指的是在这一区域中的一个block的一次移动，即从一个slot移动到另一个slot
4、retrieval：一次访问库存，相当于将区域中的一个block 拿到区域外的目的地中，如移动到卡车上等。
下面是更详细的场景描述：我们考虑一个二维的区域中，有多个相同物品（例如，在集装箱码头的集装箱）堆叠而成的堆组成。我们将这些物品称为blocks，而堆叠区域内block的可能位置称为slots。堆叠面积由其宽度（堆叠数）和高度（层数）定义。relocation实现block从一个slot到另一个slot。retrieval表示block的访问，即从当前的slot到区域外的最终目的地，例如卡车。下图可以帮助我们更直观地理解：
问题满足以下特性：
1、一个block只能从顶部访问，即在一个堆中只有位于最顶层的block才能被获取，
2、堆叠区域中的每个block必须放置在另一个block的顶部或地面（第一层），
3、堆数以限定，即重新堆放形成新的 1堆是不可行的，4、同样，每个堆栈的层数以H限定。
在模型的构建过程中，我们遵循以下假设：
1、在存储blocks时，检索序列是未知的，因此在存储时没有按照顺序堆叠，
2、当我们进行retrieval时，不会有新进入的blocks。假设一就说明了，访问优先级低的block可能会堆叠在优先级高的block上。由于block只能从顶部访问，因此我们的任务是按照给定的顺序，从堆叠区域中获取所有的blocks，而在这个过程中需要经过多次的relocations来达到目标。
到这里
我们对BPR
块重新定位问题 blocks relocation problem
已经有了一个初步的理解啦
关于blocks的retrieval和relocation的文献基本上可以分为两个主要领域，这取决于何时进行检索操作：1、允许retrieval和relocation操作同时发生（本文遵循此条件）2、第二组问题将这两个retrieval和relocation分开，首先执行blocks的relocation，然后按照规定的顺序进行retrieval
两类问题的主要区别在于堆叠二区域中的blocks的数量：而对于第一类问题，由于retrieval操作的执行，区域中的blocks数随着时间的推移而减少，在第二类问题中，区域中blocks的总数保持不变。
BRP是一些更普遍的问题的特例。如，BRP与blocks-world planning(BWP)有关。与BRP一样，BWP也是将block存储在stacks中。与目标是从堆栈中访问所有blocks的BRP相比，BWP允许任何可行的堆栈配置作为目标状态。因此，BWP是一种更通用的方法。
考虑到篇幅限制
本文仅介绍BRP的模型
对其他相关问题感兴趣的小伙伴可以自行查找资料，小编也会在参考文献中放一篇有关BWP的文章供大家参考。
建立模型
本文介绍了文献[1]提出的两种BRP模型。模型一在一定的时间范围内工作，而这个时间范围需要提前进行估计。为了增加模型本身的可用性，文献引入了一些现实的假设，从而产生了一个新模型即模型二，从而显著减少了搜索空间。
下面先介绍模型一
我们给出一个具有个堆和层的二维堆叠区域。考虑到实际情况中的空间限制，区域的最大高度()和最大宽度()是给定的参数。坐标表示区域内的每个slot，block的坐标为即表示此block在这个slot中，其中，，和，，分别表示堆叠区域内的第几堆和第几层。我们假设初始包含N个blocks，标记为1，...，，我们需要，即堆叠区域至少可以容纳个blocks，并且有足够的空间进行relocation（即至少有H-1个空的slots是为了允许通过移除位于目标block顶部的最多H-1个blocks来获取目标block）。此外，让表示获取块的优先级，，，。我们通过引入时间段，，来离散时间范围，其中每个时间段中都包含一次移动（包括relocation和retrieval）。建立模型前需要计算T的值。下面介绍一则引理，因为不是本文的重点，所以这里就不放出这个引理的证明，相关证明在文献[1]中有介绍。引理 1 ：考虑一个具有 N 个blocks、W 个堆和层数 H ≥ N - 1 的 BRP。考虑一个最优解决方案，relocation的次数的上界定义为
其中，
使用这个结果，很容易获得移动次数的上限，作为检索操作的总和加上重定位的上限：UB N。
模型需要以下两组变量：
(i)表示block位置状态的变量和，目的是定义可行的block位置状态配置(ii)表示移动的变量和，，目的是定义可行的block移动，即从一个可行的位置状态移动到下一个可行的位置状态，包括relocation和retrieval。
BRP - Ⅰ模型如下：
约束（1）保证在每个时间段内，每个block必须在堆叠区域内或在目的地区域（即已经被访问）；约束（2）保证在每个时间段内，每个slot最多放置一个block；约束（3）保证每个堆中没有间隙，即，如果在位置的block被移动，则block上方的slot必须是空的。此外，任何block都不能移动到半空的位置；约束（4）保证在每个时间段内，最多允许一次移动；约束（5）保证必须按照规定的顺序获取blocks；约束（6）为“流量平衡约束”：保证在给定时间段内的一个可行的位置状态，通过一组可行的移动，达到在时间段内的新的可行的位置状态；约束（7）保证位置状态变量与移动变量之间的关系其目标是最小化从堆区域获取所有blocks所需的时间。这一目标是通过最大化最后一个block在外部区域所花费的时间来实现的。模型BRP-Ⅰ的缺点是需要一个移动总数的上界。即必须事先计算出。然而，对于实际情况，这个边界可能仍然远离最优边界，因此使用模型BRP-Ⅰ有些不实用。因此，我们接下来提出了一个具有固定时间范围的模型BRP-Ⅱ。
我们首先引入一个附加假设:
当检索目标block时，我们只允许目标block上方的blocks进行relocation操作。例如，在下图中，当要获取block 1时，我们需要relocate block 9、18和7。此外，为了保持可行性，必须遵循LIFO策略。
虽然从实际的角度来看，这一假设是合理的，并且可以减少搜索空间，但它可能屏蔽了最优解，如下面的示例2所示。当遵循假设时，至少需要6次relocation来获取所有的blocks。但是，通过首先将block 2从堆2移动到堆1，可以获得只有4次relocations的解决方案，但是这一解决方案被假设排除在外。而BRP - Ⅰ则考虑了完整的可行域。
两个模型之间最显著的区别在于对时间段的定义和使用不同：在BRP-Ⅰ中，每个时间段中存在一次移动（relocation或retrieval）。在BRP-Ⅱ中，每个时间段t中包含block t所需的relocation（如果需要的话）和retrieval。所以在 BRP-Ⅱ 中，完成操作所需的总周期数 T 是预先知道的并且固定为 N。
在BRP-Ⅱ中，不再是一个决策变量，而是模型的一个参数，定义如下:
BRP-Ⅱ模型如下：
前5个约束与模型BRP-Ⅰ相同
这里不再赘述
约束（8）保证如果在时间段中，block 在block下面，block 被移动，那么在时间段中，block 不能在block下面找到。因为BRP-Ⅱ允许在每个时间段处理多个relocations，所以需要保证 LIFO 策略。如果在时间段 t 中，任何block从移动到，即，那么在时间段 t 中，不允许将位于位置上方的block移动到上方的位置， 即
约束（9）确保只能移动目标block上方的block。在每个时间，下一个要进行retrieval操作的是block 。因此，在每个时间段，如果对于成立，则说明目标block位于堆。在这种情况下，约束（9）不等式的左侧变为零，确保只有堆中的block可以移动；约束（10）保证在同一堆中不允许重新排列；该模型的目标函数是最小化relocation的次数。BRP - Ⅱ允许解决比BRP-Ⅰ大得多的实例。然而，尽管BRP-Ⅱ在计算时间和问题大小方面更有优势，但是BRP-Ⅰ提供了一个更详细的relocation过程，有助于提出实际的方案；且BRP - Ⅰ不需要任何额外的假设。
代码
为了突出重点，这里仅展示调用CPLEX建立模型并进行求解的代码，完整的开源项目的链接会在留言区给出。代码由深圳大学金波老师提供。这个项目是由用Python3.8和CPLEX 20.1编写的。
1. BRP-Ⅰ
from itertools import product
from docplex.mp.model import Model
from bay import Bay
from common import irange
class BRP_I: def __init__(self, bay): bay.validate_full_distinct() self.model = model = Model()
self.W = W = bay.n_stacks self.H = H = bay.n_tiers self.N = N = bay.n_blocks self.T = T = N bay.brp_min_max() #操作次数 retrieval relocation
#定义变量 self.b = b = model.binary_var_dict(product(irange(1, W), irange(1, H), irange(1, N), irange(1, T)), name='b') self.v = v = model.binary_var_dict(product(irange(1, N), irange(1, T)), name='v') self.x = x = model.binary_var_dict(product(irange(1, W), irange(1, H), irange(1, W), irange(1, H), irange(1, N), irange(1, T)), name='x') self.y = y = model.binary_var_dict(product(irange(1, W), irange(1, H), irange(1, N), irange(1, T)), name='y')
# objective model.maximize(model.sum(v[N, t] for t in irange(1, T))) # (1) model.add_constraints(model.sum(b[i, j, n, t] for i, j in product(irange(1, W), irange(1, H))) v[n, t] == 1 for n, t in product(irange(1, N), irange(1, T))) # (2) model.add_constraints(model.sum(b[i, j, n, t] for n in irange(1, N)) <= 1 for i, j, t in product(irange(1, W), irange(1, H), irange(1, T))) # (3) model.add_constraints(model.sum(b[i, j, n, t] for n in irange(1, N)) >= model.sum(b[i, j 1, n, t] for n in irange(1, N)) for i, j, t in product(irange(1, W), irange(1, H - 1), irange(1, T))) # (4) model.add_constraints(model.sum(x[i, j, k, l, n, t] for i, j, k, l, n in product(irange(1, W), irange(1, H), irange(1, W), irange(1, H), irange(1, N))) model.sum(y[i, j, n, t] for i, j, n in product(irange(1, W), irange(1, H), irange(1, N))) <= 1 for t in irange(1, T)) # (5) model.add_constraints(model.sum(v[n, t] for t in irange(1, T)) >= model.sum(v[n 1, t] for t in irange(1, T)) 1 for n in irange(1, N - 1)) # (6) model.add_constraints(b[i, j, n, t] == b[i, j, n, t - 1] model.sum(x[k, l, i, j, n, t - 1] for k, l in product(irange(1, W), irange(1, H))) - model.sum(x[i, j, k, l, n, t - 1] for k, l in product(irange(1, W), irange(1, H))) - y[i, j, n, t - 1] for i, j, n, t in product(irange(1, W), irange(1, H), irange(1, N), irange(2, T))) # (7) model.add_constraints(v[n, t] == model.sum(y[i, j, n, tt] for i, j, tt in product(irange(1, W), irange(1, H), irange(1, t - 1))) for n, t in product(irange(1, N), irange(1, T))) # pre-processing model.add_constraints(b[i, j, n, 1] == int(bay.pri[i - 1][j - 1] == n) for i, j, n in product(irange(1, W), irange(1, H), irange(1, N)))
#每操作一次形成一个bay，便于输出（图形） def get_bays(self): bays = {} for t in irange(1, self.T): conf = [[None] * self.H for _ in range(self.W)] for i, j, n in product(irange(1, self.W), irange(1, self.H), irange(1, self.N)): if round(self.b[i, j, n, t].solution_value) == 1: conf[i - 1][j - 1] = n bays[t] = Bay(self.W, self.H, conf) return bays
def get_n_relos(self): return sum(round(self.x[i, j, k, l, n, t].solution_value) for i, j, k, l, n, t in product(irange(1, self.W), irange(1, self.H), irange(1, self.W), irange(1, self.H), irange(1, self.N), irange(1, self.T)))
#测试数据def test(): conf = [[1, 3, 4], [5], [2]] bay = Bay(3, 3, conf) print('bay') print(bay)
brp_i = BRP_I(bay) if brp_i.model.solve(): print() print('n_relos = {}'.format(brp_i.get_n_relos())) bays = brp_i.get_bays() for t in irange(1, brp_i.T): print('t = {}'.format(t)) print(bays[t])
if __name__ == '__main__': test()2. BRP-Ⅱ
from itertools import product
from docplex.mp.model import Model
from bay import Bayfrom common import irange
class BRP_II: def __init__(self, bay, bug1_fixed=True, bug2_fixed=True): bay.validate_full_distinct() self.model = model = Model()
self.W = W = bay.n_stacks self.H = H = bay.n_tiers self.N = N = bay.n_blocks
self.b = b = model.binary_var_dict(product(irange(1, W), irange(1, H), irange(1, N), irange(1, N)), name='b') self.x = x = model.binary_var_dict(product(irange(1, W), irange(1, H), irange(1, W), irange(1, H), irange(1, N), irange(1, N)), name='x') self.y = y = model.binary_var_dict(product(irange(1, W), irange(1, H), irange(1, N), irange(1, N)), name='y')
self.v = v = {(n, t): int(n < t) for n, t in product(irange(1, N), irange(1, N))}
# objective model.minimize(model.sum(x[i, j, k, l, n, t] for i, j, k, l, n, t in product(irange(1, W), irange(1, H), irange(1, W), irange(1, H), irange(1, N), irange(1, N)))) # (1) model.add_constraints(model.sum(b[i, j, n, t] for i, j in product(irange(1, W), irange(1, H))) v[n, t] == 1 for n, t in product(irange(1, N), irange(1, N))) # (2) model.add_constraints(model.sum(b[i, j, n, t] for n in irange(1, N)) <= 1 for i, j, t in product(irange(1, W), irange(1, H), irange(1, N))) # (3) model.add_constraints(model.sum(b[i, j, n, t] for n in irange(1, N)) >= model.sum(b[i, j 1, n, t] for n in irange(1, N)) for i, j, t in product(irange(1, W), irange(1, H - 1), irange(1, N))) # (6) model.add_constraints(b[i, j, n, t] == b[i, j, n, t - 1] model.sum(x[k, l, i, j, n, t - 1] for k, l in product(irange(1, W), irange(1, H))) - model.sum(x[i, j, k, l, n, t - 1] for k, l in product(irange(1, W), irange(1, H))) - y[i, j, n, t - 1] for i, j, n, t in product(irange(1, W), irange(1, H), irange(1, N), irange(2, N))) # (7) model.add_constraints(v[n, t] == model.sum(y[i, j, n, tt] for i, j, tt in product(irange(1, W), irange(1, H), irange(1, t - 1))) for n, t in product(irange(1, N), irange(1, N))) # (8) Constraints (8) in BRP-II do not allow relocating several blocks to the same destination stack. if bug1_fixed: model.add_constraints((H - 1) * (1 - model.sum(x[i, j, k, l, n, t] for n in irange(1, N))) >= model.sum(x[i, jj, k, ll, n, t] for n, jj, ll in product(irange(1, N), irange(j 1, H), irange(l 1, H))) for i, j, k, l, t in product(irange(1, W), irange(1, H), irange(1, W), irange(1, H), irange(1, N - 1))) else: model.add_constraints(1 - model.sum(x[i, j, k, l, n, t] for n in irange(1, N)) >= model.sum(x[i, jj, k, ll, n, t] for jj, ll, n in product(irange(j 1, H), irange(l 1, H), irange(1, N))) for i, j, k, l, t in product(irange(1, W), irange(1, H), irange(1, W), irange(1, H), irange(1, N - 1))) # (9) Constraints (9) in BRP-II allow relocating blocks below the target block. if bug2_fixed: model.add_constraints((H - 1) * (1 - b[i, j, t, t]) >= model.sum(x[ii, jj, k, l, n, t] for ii in irange(1, W) for jj in irange(1, W if ii != i else j - 1) for k, l, n in product(irange(1, W), irange(1, H), irange(1, N))) for i, j, t in product(irange(1, W), irange(1, H), irange(1, N))) else: model.add_constraints((H - 1) * (1 - model.sum(b[i, j, t, t] for j in irange(1, H))) >= model.sum(x[ii, j, k, l, n, t] for ii in irange(1, W) if ii != i for j, k, l, n in product(irange(1, H), irange(1, W), irange(1, H), irange(1, N))) for i, t in product(irange(1, W), irange(1, N))) # (10) model.add_constraints(x[i, j, i, l, n, t] == 0 for i, j, l, n, t in product(irange(1, W), irange(1, H), irange(1, H), irange(1, N), irange(1, N)))
# pre-processing model.add_constraints(b[i, j, n, 1] == int(bay.pri[i - 1][j - 1] == n) for i, j, n in product(irange(1, W), irange(1, H), irange(1, N)))
def get_bays(self): bays = {} for t in irange(1, self.N): conf = [[None] * self.H for _ in range(self.W)] for i, j, n in product(irange(1, self.W), irange(1, self.H), irange(1, self.N)): if round(self.b[i, j, n, t].solution_value) == 1: conf[i - 1][j - 1] = n bays[t] = Bay(self.W, self.H, conf) return bays
def get_n_relos(self): return self.model.objective_value
#测试数据1def test1(): conf = [[1, 3, 4], [5], [2]] bay = Bay(3, 3, conf) print('bay') print(bay)
brp_ii = BRP_II(bay) if brp_ii.model.solve(): print() print('n_relos = {}'.format(brp_ii.get_n_relos())) bays = brp_ii.get_bays() for t in irange(1, brp_ii.N): print('t = {}'.format(t)) print(bays[t])
brp_ii = BRP_II(bay, bug1_fixed=False) if brp_ii.model.solve(): print() print('# bug 1 unfixed #') print('n_relos = {}'.format(brp_ii.get_n_relos())) bays = brp_ii.get_bays() for t in irange(1, brp_ii.N): print('t = {}'.format(t)) print(bays[t])
#测试数据2def test2(): conf = [[4], [3, 1], [2, 5, 6]] bay = Bay(3, 3, conf) print('bay') print(bay)
brp_ii = BRP_II(bay) if brp_ii.model.solve(): print() print('n_relos = {}'.format(brp_ii.get_n_relos())) bays = brp_ii.get_bays() for t in irange(1, brp_ii.N): print('t = {}'.format(t)) print(bays[t])
brp_ii = BRP_II(bay, bug2_fixed=False) if brp_ii.model.solve(): print() print('# bug 2 unfixed #') print('n_relos = {}'.format(brp_ii.get_n_relos())) bays = brp_ii.get_bays() for t in irange(1, brp_ii.N): print('t = {}'.format(t)) print(bays[t])
if __name__ == '__main__': test1() print() test2()结果展示
更多的算例
也可以在文末留言区获得
1、算例
2. BRP-Ⅰ结果
n_relos = 2
3. BRP-Ⅱ结果
bug 1 fixed
n_relos = 2
bug 1 unfixed
n_relos = 3
参考文献
[1]Caserta, M., Schwarze, S., & Voß, S. (2012). A mathematical formulation and complexity considerations for the blocks relocation problem. European Journal of Operational Research, 219(1), 96-104.
[2]Gupta, N., Nau, D., 1992. On the complexity of blocks-world planning. Artifical Intelligence 56 (2–3), 223–254.
-The End-
文案&代码注释：胡心瑶（华中科技大学管理学院本科三年级）指导老师：秦虎老师 （华中科技大学管理学院） 金波老师（深圳大学管理学院）排版：程欣悦(荆楚理工学院四年级)如对文中内容有疑问，欢迎交流。PS：部分资料来自网络。
如有需求，可以联系：
秦虎老师（华中科技大学管理学院，微信号43340630）
胡心瑶（华中科技大学管理学院本科三年级：1603945923@qq.com）
程欣悦（荆楚理工学院四年级：1293900389@qq.com)
欢迎大家加入数据魔术师粉丝群，我们的活动将会通过粉丝群优先发布， 学习资料将通过粉丝群分享。
欲入群，请转发此文，然后扫描下方二维码联系数据魔术师小助手
