TransportationProblem 0.0.3

求解产销平衡的运输问题。

TransportationProblem

求解产销平衡的运输问题。

Getting Started

Install this package:

$ pip install TransportationProblem

Utilize it to solve a transportation problem:

>>> import transportation_problem as tp
>>> s = [('A1', 14), ('A2', 27), ('A3', 19)]	# 产地
>>> d = [('B1', 22), ('B2', 13), ('B3', 12), ('B4', 13)]	# 销地
>>> c = [[6, 7, 5, 3], [8, 4, 2, 7], [5, 9, 10, 6]]	# 运价
>>> p = tp.TransportationProblem(s, d, c)
>>> r = p.solve()
>>> print(r)
Transportation problem optimized successfully. Result cost (total): 232.0
      运量	      B1	      B2	      B3	      B4	
      A1	     1.0	     0.0	     0.0	    13.0	
      A2	     2.0	    13.0	    12.0	     0.0	
      A3	    19.0	     0.0	     0.0	     0.0	

产销平衡的运输问题

设某种物品有m个产地A1, A2，· · ·, Am，各产地的产量分别是a1, a2, · · ·, am；有n个销地B1，B2，· · ·，Bn，各个销地的销量分别为b1，b2，· · ·，bn。假定从产地Ai(i=1,2,· · ·,m)向销地Bj(j=1,2,· · ·,n)运输单位物品的运价为c_ij，问怎么调运这些物品才能使总运费最小？

若在这个运输问题中，各产地产量之和等于各销地销量之和，则称该运输问题产销平衡。

要解决一个非产销平衡的运输问题，往往需要人工对问题做出调整，使其成为产销平衡的运输问题，然后才能求解。所以我们只考虑产销平衡的运输问题的求解方法。

则称该运输问题产销平衡。

对于一个产销平衡的运输问题，已知的信息可以用下表表示：

产地\销地	B1	· · ·	Bn	产量
A1	c_11	· · ·	c_1n	a1
:	:	(运价)	:	:
Am	c_m1	· · ·	c_mn	am
销量	b1	· · ·	bn	sum(a_i) ==sum(b_j)

这样的运输问题，是一类具有特殊结构的线性规划问题，所以我们完全可以采用 单纯形法 求解。

单纯形解法

一个产销平衡的运输问题可以转化为如下线性规划问题：

$$ \min z = \sum_{i=1}^m\sum_{j=1}^n c_{ij}x_{ij}\ s.t. \left{\begin{array}{l} \sum_{j=1}^n x_{ij} = a_i\ \sum_{i=1}^m x_{ij} = b_j\ x_{ij} \ge 0 \end{array}\right. $$

对此线性规划问题使用单纯形算法求解即可得到最优运输方案（如果可行）。

由于这个方法只是对单纯形算法的简单应用，在该项目中不做实现。但我在我的单纯形算法项目 SimplexLinprog 中给出了这一解法的实现：

• https://github.com/cdfmlr/SimplexLinprog/blob/master/transportation_problem.py

由于运输问题的特殊性，我们可以将单纯性法简化，用表格来处理。这个解法称为「表上作业法」。

表上作业法

该项目实现了使用表上作业法对产销平衡的运输问题的最优化求解。表上作业法的步骤如下：

	描述	方法
第一步	求初始基行可行解(初始调运方案)	最小元素法、 西北角法、 伏格尔法
第二步	求检验数并判断是否得到最优解. 当非基变量的检验数σ_ij≥0时得到最优解， 若存在检验数σ_ij <0，说明还没有达到最优，转第三步.	闭回路法、 位势法
第三步	调整运量，即基变换. 选一个变量出基，对原运量进行调整得到新的基可行解，转第二步.	闭回路调整法

各步骤中具体的实现方法有：

• 求初始基行可行解
  • 最小元素法：从运价最小的地方开始供应(调运)，然后次小， 直到最后供完为止。
  • 西北角法(或左上角法)：从西北角（左上角）格开始，在格内的右下角标上允许取得的最大数；然后按行（列）标下一格的数；若某行（列）的产量（销量）已满足，则把该行（列）的其他格划去；如此进行下去，直至得到一个基本可行解的方法。
  • 伏格尔法：
    1. 从运价表中分别计算出各行和各列的最小运费和次最小运费的差额，并填入该表的最右列和最下行。
    2. 再从差额最大的行或列中找出最小运价确定供需关系和供需数量（有相同最小运价时，供需数量尽可能大！）。当产地或销地中有一方数量供应完毕或得到满足时，划去运价表中对应的行或列。重复1）和2）， 直到找出初始运输方案为止。
• 求检验数，判断是否得到最优解
  • 闭回路法：这种方法需要对每一个空格寻找一条闭回路，并根据闭回路求出每个空格的检验数。当运输问题中m和n较大时，计算检验数的工作量很大。
  • 位势法：先对初始调运方案求出位势，然后求各空格的检验数。当所有的检验数均为非负时，就得到最优方案。
• 调整运量
  • 闭回路调整法：从检验数为负的空格出发，作闭回路，作进一步调整。

表上作业法的实现

框架设计

首先，将「表上作业法」的步骤抽象如下：

步骤	抽象类	具体实现子类
求初始基行可行解	TransportationIniter	最小元素法MinimumElementIniter; 西北角法NorthwestCornerIniter; 伏格尔法VogelIniter
求检验数， 判断是否达到最优解	TransportationChecker	位势法PotentialChecker; 闭回路法(这个效率不高，没太大意义，所以暂未实现)
调整运量	TransportationOptimizer	闭回路调整法ClosedLoopAdjustmentOptimizer

为了实现闭回路调整法ClosedLoopAdjustmentOptimizer，我们还需要实现「闭回路」的表示：

• 闭回路的节点 ClosedLoopNode
• 找闭回路的方法 ClosedLoopMethod

再看运输问题的表达：

我设计了一个 TransportationProblem 类来表示一个运输问题。它具有属性：

• supply: 产地: [('name', 产量), ...]
• demand: 销地: [('name', 销量), ...]
• costs: 运价: [[1, 2, 3, ...], [4, 5, 6, ...], ...]
• result: 运输问题求解结果，在调用 transportationProblem.solve(...) ，求解成功后才有非 None 值。

TransportationProblem 只有一个 solve() 方法，调用这个方法来对运输问题求解，这个方法需要给定表上作业法的各个步骤要使用的具体方法。也就是说，需要传入 TransportationIniter、TransportationChecker、TransportationOptimizer 的子类各一个作为参数。

若求解成功，则 solve() 会返回一个 TransportationResult 对象，这个对象中包含了:

• 求解结果（运量表）：transportation
• 总的运输成本：total_cost
• 一个原问题的引用： problem

同时，TransportationResult 重写了 __str__()，可以用 print(transportationResult) 打印出一张漂亮的运量表。

具体算法

这些涉及到的大部分算法还是比较简单的。

其中比较复杂的就是「闭回路」的寻找，这个东西还是有很多种方法可以解决的。我写了一个递归回溯的算法（这有点类似于我们做游戏的时候写敌人AI移动脚本里用的一些方法），从各种方面来看，这其实都不是一个很好的算法，但是它比较有趣。

具体各个算法的实现，可以参考这张思维导图：

（图片不清晰可以看如下大纲）

1. 求初始基行可行解

最小元素法:

• 将运价从小到大排序

• 从最小元素开始填运量

  • 1. 获取最小元素坐标
  • 2. 检查有没有"划掉"，即对应supply、demand是否为0
  • 3. 没划掉，则分配这个位置的运量

• 直到所有supply、demand值都为0

伏格尔法:

• 获取行、列差额

  • 1. 找出行/列中最低、次低cost
  • 2. 返回 diff = 次低 - 最低
  • 3. 若某行/列已经被"划掉"，则返回中对应的值为 -1

• 找出最大差额所在的行/列

• 安排这个行/列里最低运价处的运量

西北角法:

• 从西北角，即 (r, c) = (0, 0) 开始处理

• 安排(r, c)处运量

• 安排后如果行满足了(划掉行)则r+1 如果列满足了(划掉列)则c+1

• 直到(r, c)达到可取的最大值（右下角）

2. 求检验数, 判断是否达到最优解

位势法:

• 初始化位势u、v
• 对基变量：$\sigma_{ij} = c_{ij} - (u_i + v_j) = 0$，求出 u 和 v
• 计算非基变量检验数: $\sigma_{ij} = c_{ij} - (u_i + v_j)$
• 若所有的检验数均为非负时，就得到最优方案。

3. 调整运量

闭回路调整法:

• 找出所有<0的检验数， 选其中最小的

• 找这个最小的检验数处开始的一条闭回路

  • 一步一步递归

    • 1. 看看有没走完闭回路
      • 如果上一步的行/列不等于起点的行/列， 而现在这一步的行/列不等于起点的行/列， 就说明“走完了”（可以连起一条闭回路）
      • 走完了就记录当前节点，返回并通知递归的上一层记录它的自己（沿着递归走来的路径把整条闭回路记录下来）
    • 2. 没走完，则求出上一步的「方向」
    • 3. 查找如果不转向继续往前走下去可能的下一个「落脚点」
      • 「落脚点」即从当前点(作业表的某一单元格)向“前方”画射线，射线上若有基变量则可作为一个落脚点
    • 4. 查找如果转90度，可能的「落脚点」
    • 5. 从 3、4 中找到的所有可能落脚点中一个个走了试试看

      • 递归，下一步

      • 等递归返回回来了

        • 若返回中携带表示成功走完的flag， 则记录自己这一步所处的位置， 返回上一层（也携带表示成功走完的flag）
        • 若返回中携带的是表示这条路走不通的失败flag，则继续尝试走下一个可能落脚点
    • 6. 如果尝试了所有可能落脚点，都不可行，则说明上就一步走错了

      • 回溯

        • 返回表示这条路走不通的失败flag

• 在闭回路对运量进行调整

  • 算出闭回路中下标(从0开始，0是检验数<0的非基变量)为奇数处的运量的最小值min_trans
  • 调整，沿着闭回路将运价加/减min_trans. 下标从0开始，所以是偶加奇减 (如：0+, 1-, 2+, 3-)

接口使用

接口的使用还是比较简单的，直接看一个例子：

已知如下产销平衡运输问题，

产地\销地	B1	B2	B3	B4	产量
A1	6	7	5	3	14
A2	8	4	2	7	27
A3	5	9	10	6	19
销量	22	13	12	13

编写驱动脚本（注意保证 transportation_problem 在导入路径中！）：

import transportation_problem as tp

# 定义产地、产量
s = [('A1', 14), ('A2', 27), ('A3', 19)]
# 定义售量、销量
d = [('B1', 22), ('B2', 13), ('B3', 12), ('B4', 13)]
# 运价表
c = [[6, 7, 5, 3], [8, 4, 2, 7], [5, 9, 10, 6]]

# 初始化运输问题对象
p = tp.TransportationProblem(s, d, c)
# 求解，使用西北角法初始化，位势法检验，闭回路法优化调整
r = p.solve(tp.NorthwestCornerIniter, tp.PotentialChecker, tp.ClosedLoopAdjustmentOptimizer)

# 输出解
print(r)

运行，得到输出的最优运量表：

Transportation problem optimized successfully. Result cost (total): 232.0
      运量	      B1	      B2	      B3	      B4	
      A1	     1.0	     0.0	     0.0	    13.0	
      A2	     2.0	    13.0	    12.0	     0.0	
      A3	    19.0	     0.0	     0.0	     0.0	

当然，在日常使用中，我们（一般情况下）不想自己指定方法，所以 solve 方法的实现是带有默认参数的：

def solve(self, initer_class=MinimumElementIniter, checker_class=PotentialChecker, optimizer_class=ClosedLoopAdjustmentOptimizer):
    pass

使用默认参数，代码可以更加简洁：

import transportation_problem as tp

s = [('A1', 14), ('A2', 27), ('A3', 19)]
d = [('B1', 22), ('B2', 13), ('B3', 12), ('B4', 13)]
c = [[6, 7, 5, 3], [8, 4, 2, 7], [5, 9, 10, 6]]

p = tp.TransportationProblem(s, d, c)
r = p.solve()	# 默认最小元素法初始化，位势法检验，闭回路法优化调整
print(r)

另一个例子：

再看一个产销不平衡问题，来自清华大学《􏰄􏰅􏰆􏰄􏰅􏰆运筹学 第四版》的习题。产销不平衡首先要转化为产销平衡问题才能开始求解：

编写代码，求解 表4-28 的最优调运方案：

import transportation_problem as tp

spy = [('A', 400), ('B', 450), ('C', 70)]
dmd = [('甲', 290), ('甲\'', 30), ('乙', 250), ('丙', 270), ('丙\'', 80)]
cst = [[15, 15, 18, 22, 22], [21, 21, 25, 16, 16], [999, 0, 999, 999, 0]]

pbm = tp.TransportationProblem(spy, dmd, cst)
res = pbm.solve(tp.VogelIniter)		# 伏格尔法初始化

print(res)

对于 M，我们编程时可以给定一个对于当前问题来说充分大的值即可。具体到这个问题中，999 已经足够大。

当然，你也可以使用一个无穷 (numpy.inf) 来表示 M，这样也能求出最优方案，不过 inf 值会导致程序无法自动求出最小运价（输出中会显示Result cost (total): nan），所以不推荐这样：

import numpy as np
...
cst = [[15, 15, 18, 22, 22], [21, 21, 25, 16, 16], [np.inf, 0, np.inf, np.inf, 0]]
...

运行，得到结果：

Transportation problem optimized successfully. Result cost (total): 14650.0
      运量	       甲	      甲'	       乙	       丙	      丙'	
       A	   150.0	     0.0	   250.0	     0.0	     0.0	
       B	   140.0	     0.0	     0.0	   270.0	    40.0	
       C	     0.0	    30.0	     0.0	     0.0	    40.0	

输出的运量表即为最优调运方案，最小运价为 14650.0 万元。

开放源代码

MIT License

Copyright (c) 2020 CDFMLR