求解产销平衡的运输问题。
Project description
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 中给出了这一解法的实现:
由于运输问题的特殊性,我们可以将单纯性法简化,用表格来处理。这个解法称为「表上作业法」。
表上作业法
该项目实现了使用表上作业法对产销平衡的运输问题的最优化求解。表上作业法的步骤如下:
描述 | 方法 | |
---|---|---|
第一步 | 求初始基行可行解(初始调运方案) | 最小元素法、 西北角法、 伏格尔法 |
第二步 | 求检验数并判断是否得到最优解. 当非基变量的检验数σ_ij≥0 时得到最优解, 若存在检验数σ_ij <0 ,说明还没有达到最优,转第三步. |
闭回路法、 位势法 |
第三步 | 调整运量,即基变换. 选一个变量出基,对原运量进行调整得到新的基可行解,转第二步. | 闭回路调整法 |
各步骤中具体的实现方法有:
- 求初始基行可行解
- 最小元素法:从运价最小的地方开始供应(调运),然后次小, 直到最后供完为止。
- 西北角法(或左上角法):从西北角(左上角)格开始,在格内的右下角标上允许取得的最大数;然后按行(列)标下一格的数;若某行(列)的产量(销量)已满足,则把该行(列)的其他格划去;如此进行下去,直至得到一个基本可行解的方法。
- 伏格尔法:
- 从运价表中分别计算出各行和各列的最小运费和次最小运费的差额,并填入该表的最右列和最下行。
- 再从差额最大的行或列中找出最小运价确定供需关系和供需数量(有相同最小运价时,供需数量尽可能大!)。当产地或销地中有一方数量供应完毕或得到满足时,划去运价表中对应的行或列。重复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. 求初始基行可行解
最小元素法:
-
将运价从小到大排序
-
从最小元素开始填运量
-
- 获取最小元素坐标
-
- 检查有没有"划掉",即对应supply、demand是否为0
-
- 没划掉,则分配这个位置的运量
-
-
直到所有supply、demand值都为0
伏格尔法:
-
获取行、列差额
-
- 找出行/列中最低、次低cost
-
- 返回 diff = 次低 - 最低
-
- 若某行/列已经被"划掉",则返回中对应的值为 -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的检验数, 选其中最小的
-
找这个最小的检验数处开始的一条闭回路
-
一步一步递归
-
- 看看有没走完闭回路
- 如果上一步的行/列不等于起点的行/列, 而现在这一步的行/列不等于起点的行/列, 就说明“走完了”(可以连起一条闭回路)
- 走完了就记录当前节点,返回并通知递归的上一层记录它的自己(沿着递归走来的路径把整条闭回路记录下来)
-
- 没走完,则求出上一步的「方向」
-
- 查找如果不转向继续往前走下去可能的下一个「落脚点」
- 「落脚点」即从当前点(作业表的某一单元格)向“前方”画射线,射线上若有基变量则可作为一个落脚点
-
- 查找如果转90度,可能的「落脚点」
-
- 从 3、4 中找到的所有可能落脚点中一个个走了试试看
-
递归,下一步
-
等递归返回回来了
- 若返回中携带表示成功走完的flag, 则记录自己这一步所处的位置, 返回上一层(也携带表示成功走完的flag)
- 若返回中携带的是表示这条路走不通的失败flag,则继续尝试走下一个可能落脚点
-
- 如果尝试了所有可能落脚点,都不可行,则说明上就一步走错了
-
回溯
- 返回表示这条路走不通的失败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
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distribution
Hashes for TransportationProblem-0.0.3.tar.gz
Algorithm | Hash digest | |
---|---|---|
SHA256 | 8eebb2c1f4f90a5d2b06d49b03deb3de29f553d3fda4505dd78e56ef4b9796d3 |
|
MD5 | 236027d0ec0c2cab6f7b000165cb3112 |
|
BLAKE2b-256 | e76b097441a1173642ec0ea0e48fb838d022cb994b35eb79463e904fedb3cb92 |