TransportationProblem 0.0.3
求解产销平衡的运输问题。
License: MIT

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

## Getting Started

Install this package:

```shell
$ pip install TransportationProblem
```

Utilize it to solve a transportation problem:

```python
>>> 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) |

这样的运输问题，是一类具有特殊结构的[线性规划](https://baike.baidu.com/item/线性规划/1661606)问题，所以我们完全可以采用 [单纯形法](https://baike.baidu.com/item/单纯形法/8580570) 求解。

### 单纯形解法

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

$$
\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) 中给出了这一解法的实现：

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

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

### 表上作业法

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

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

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

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

## 表上作业法的实现

### 框架设计

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

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

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

- 闭回路的节点 `ClosedLoopNode`
- 找闭回路的方法 `ClosedLoopMethod`

![屏幕快照 2020-04-17 12.06.27](https://tva1.sinaimg.cn/large/007S8ZIlgy1gdwmwzweyej30u00xv129.jpg)

再看运输问题的表达：

我设计了一个 `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)` 打印出一张漂亮的运量表。

![image-20200417115548306](https://tva1.sinaimg.cn/large/007S8ZIlgy1gdwmkgq1odj30jy09mq3u.jpg)

### 具体算法

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

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

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

![image-20200417155635563](https://tva1.sinaimg.cn/large/007S8ZIlgy1gdwtiyiouaj30xt0u012c.jpg)

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

#### 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` 在导入路径中！）：

```python
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` 方法的实现是带有默认参数的：

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

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

```python
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)
```

**另一个例子**：

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

![IMG_0B78AC4053D1-1](https://tva1.sinaimg.cn/large/007S8ZIlgy1gdx6l6qtzkj30v90jsk0g.jpg)

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

```python
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