用Octave计算运输问题（表上作业法）
广告
{{v.name}}
表上作业法是运输问题的经典手工解法，仅需在运输表上完成三步操作，无需复杂计算，适用于中小规模问题（\(m,n≤5\)），核心逻辑为：找初始可行调运方案→检验方案是否最优→非最优则调整方案，循环至得到最优解。
核心前提
平衡运输问题（不平衡需先转化为平衡），且初始可行方案的基变量数满足 \(m+n-1\)（无退化，退化时需补0虚基变量）。
步骤1：确定初始基本可行调运方案
即找到满足供需约束的调运量 \(x_{ij}\)，共3种常用方法，最小元素法最简单，伏格尔法（VAM） 初始方案更接近最优解（推荐）。
方法1：最小元素法（贪心，从最小单位运费开始调运）
1. 在运输表中找到最小的单位运费\(c_{ij}\)，优先给该单元格分配调运量，调运量为 \(\min(产地A_i的剩余产量, 销地B_j的剩余销量)\)；
2. 若产地\(A_i\)产量耗尽，划去该行；若销地\(B_j\)销量满足，划去该列；若同时耗尽/满足，划去一行一列（仅留其一，避免基变量数不足）；
3. 在剩余的运输表中重复步骤1-2，直到所有产地产量耗尽、销地销量满足。
方法2：伏格尔法（VAM，考虑最小运费与次小运费的差值，减少总运费偏差）
1. 计算运输表中每行的罚数（该行最小单位运费与次小单位运费的差值）、每列的罚数（该列最小与次小的差值），罚数表示“放弃最小运费调运的损失”；
2. 找到最大的罚数（行/列），在该罚数对应的行/列中选择最小的单位运费\(c_{ij}\)，分配调运量 \(\min(剩余产量, 剩余销量)\)；
3. 划去耗尽产量的行/满足销量的列，重复步骤1-2，直到所有供需平衡。
伏格尔法的初始方案总运费通常远小于最小元素法，后续调整次数更少，优先使用。
步骤2：检验初始方案是否为最优解
通过计算非基变量（未分配调运量的单元格）的检验数\(\sigma_{ij}\)判断，若所有检验数\(\sigma_{ij}≥0\)（最小化问题），则当前方案为最优解；若存在\(\sigma_{ij}< 0\)，则方案非最优，需调整。
常用检验方法：闭回路法（直观，手工首选）、位势法（高效，适用于规模稍大的问题）。
方法1：闭回路法（为每个非基变量构建唯一的闭回路，计算检验数）
1. 对于每个非基变量单元格\((i,j)\)，以该单元格为起点，构建闭回路：仅沿水平/垂直方向移动，仅在基变量单元格（已分配调运量）和该非基变量单元格转弯，形成唯一的封闭折线；
2. 闭回路的顶点依次标为\(+,-,+,-...\)（起点非基变量为\(+\)）；
3. 计算检验数 \(\sigma_{ij} = \sum(闭回路上+号顶点的单位运费) - \sum(闭回路上-号顶点的单位运费)\)；
4. 遍历所有非基变量，计算其检验数。
方法2：位势法（引入行位势\(u_i\)、列位势\(v_j\)，避免构建闭回路，计算更快）
1. 对于所有基变量单元格\((i,j)\)，满足 \(u_i + v_j = c_{ij}\)（位势方程）；
2. 令某一行位势\(u_1=0\)（任意行均可），依次求解所有行位势\(u_i\)和列位势\(v_j\)；
3. 对于所有非基变量单元格\((i,j)\)，检验数 \(\sigma_{ij} = c_{ij} - (u_i + v_j)\)；
位势法与闭回路法的检验数结果完全一致，无需构建闭回路，效率更高。
步骤3：调整非最优方案（闭回路调整法）
1. 找到最小的负检验数\(\sigma_{kl}\)（对应非基变量单元格\((k,l)\)），该单元格为调入格（即将该非基变量转为基变量，分配调运量）；
2. 为调入格\((k,l)\)构建闭回路，标上\(+,-\)号（起点\(+\)）；
3. 在闭回路的所有\(-\)号顶点中，找到最小的调运量\(x_{ij}\)，该单元格为调出格（即将该基变量转为非基变量，调运量置0）；
4. 调整闭回路上的调运量：\(+\)号顶点加该最小调运量，\(-\)号顶点减该最小调运量，其余单元格调运量不变；
5. 对调整后的新方案，重复步骤2（检验），直到所有检验数≥0，得到最优解。
关键补充：退化问题处理
若初始方案的基变量数\( < m+n-1\)，称为退化，此时无法构建闭回路/求解位势，需补0虚基变量：在运输表中选择非基变量单元格（通常选单位运费最小的），分配调运量\(x_{ij}=0\)，使其成为虚基变量，保证基变量数=\(m+n-1\)，再进行后续检验。
例子：平衡运输问题
已知条件
3个产地\(A_1,A_2,A_3\)，产量分别为\(30,50,20\)；2个销地\(B_1,B_2\)，销量分别为\(60,40\)；总供应量=总需求量=100（平衡）。
单位运费表（\(c_{ij}\)）：
| 产地\销地 | \(B_1\)（60） | \(B_2\)（40） | 产量\(a_i\) |
|----------|------------|------------|-----------|
| \(A_1\) | 2 | 3 | 30 |
| \(A_2\) | 1 | 4 | 50 |
| \(A_3\) | 5 | 2 | 20 |
目标：求总运费最小的调运方案。
步骤1：伏格尔法找初始可行方案
第一步：计算初始罚数
- 行罚数：\(A_1\):\(\min(2,3)=2\)，次小3→罚数1；\(A_2\):\(\min(1,4)=1\)，次小4→罚数3；\(A_3\):\(\min(5,2)=2\)，次小5→罚数3；
- 列罚数：\(B_1\):\(\min(2,1,5)=1\)，次小2→罚数1；\(B_2\):\(\min(3,4,2)=2\)，次小3→罚数1；
最大罚数为\(A_2/A_3\)（3），选\(A_2\)（最小运费更小），\(A_2\)的最小运费为\(c_{21}=1\)（\(B_1\)），调运量\(\min(50,60)=50\)。
- \(A_2\)产量耗尽，划去\(A_2\)行；\(B_1\)剩余销量=60-50=10。
第二步：剩余表计算罚数
剩余产地\(A_1,A_3\)，剩余销地\(B_1\)(10)、\(B_2\)(40)
- 行罚数：\(A_1\):1，\(A_3\):3；列罚数：\(B_1\):\(\min(2,5)=2\)→罚数3，\(B_2\):\(\min(3,2)=2\)→罚数1；
最大罚数为\(A_3/B_1\)（3），选\(B_1\)列，最小运费\(c_{11}=2\)（\(A_1\)），调运量\(\min(30,10)=10\)。
- \(B_1\)销量满足，划去\(B_1\)列；\(A_1\)剩余产量=30-10=20。
第三步：剩余表仅余\(B_2\)列
剩余调运量：\(A_1\)剩余20，\(A_3\)20，\(B_2\)销量40，直接分配：
- \(A_1→B_2\)：20，\(A_3→B_2\)：20，\(B_2\)销量满足，所有供需平衡。
初始调运方案（基变量数=3+2-1=4，无退化）
\(x_{11}=10, x_{12}=20, x_{21}=50, x_{32}=20\)；总运费 \(Z=10×2+20×3+50×1+20×2=20+60+50+40=170\)
步骤2：位势法检验最优性
1. 设\(u_1=0\)，基变量满足\(u_i+v_j=c_{ij}\)：
- \(x_{11}\)：\(u_1+v_1=2\) → \(v_1=2\)；
- \(x_{12}\)：\(u_1+v_2=3\) → \(v_2=3\)；
- \(x_{21}\)：\(u_2+v_1=1\) → \(u_2=-1\)；
- \(x_{32}\)：\(u_3+v_2=2\) → \(u_3=-1\)；
2. 计算非基变量检验数（\(\sigma_{ij}=c_{ij}-(u_i+v_j)\)）：
- 非基变量\(x_{22}\)：\(\sigma_{22}=4-(-1+3)=2≥0\)；
- 非基变量\(x_{31}\)：\(\sigma_{31}=5-(-1+2)=4≥0\)；
所有检验数≥0，初始方案即为最优解，最小总运费=170。
代码如下：
function [X, totalCost, iter] = bszy(supply, demand, C, method, varargin)
% BSZY 表上作业法（运输问题）
% [X, totalCost, iter] = bszy(supply, demand, C, method, varargin)
% supply: m×1 数组，产地产量（m 个产地）
% demand: n×1 数组，销地销量（n 个销地）
% C: m×n 矩阵，单位运费矩阵（c(i,j)=A_i→B_j的单位运费）
% method: 'vogel' 或 'nw' (默认 'vogel')
if nargin < 4 || isempty(method)
method = 'vogel';
end
% 可选参数：verbose (传 true 打印调试信息)
verbose = false;
if nargin >= 5 && ~isempty(varargin)
if islogical(varargin{1}) && varargin{1}
verbose = true;
end
end
% 快速模式：'auto' 同时尝试 'vogel' 和 'nw'，返回代价更小的结果
if strcmpi(method, 'auto')
[Xv, Cv, itv] = bszy(C, supply, demand, 'vogel', verbose);
[Xn, Cn, itn] = bszy(C, supply, demand, 'nw', verbose);
if Cv <= Cn
X = Xv; totalCost = Cv; iter = itv; return;
else
X = Xn; totalCost = Cn; iter = itn; return;
end
end
[m, n] = size(C);
% 保存原始成本矩阵，用于最后计算总费用并去除虚拟行列
origC = C;
orig_m = m; orig_n = n;
% 平衡供需：若总供给不等于总需求，加入虚拟行或列（成本为0）
sumS = sum(supply);
sumD = sum(demand);
dummyType = 'none';
if abs(sumS - sumD) > 1e-8
if sumS > sumD
% 加入虚拟需求列
diff = sumS - sumD;
C(:,end+1) = 0;
demand(end+1) = diff;
dummyType = 'col';
dummyIdx = size(C,2);
else
% 加入虚拟供给行
diff = sumD - sumS;
C(end+1,:) = 0;
supply(end+1) = diff;
dummyType = 'row';
dummyIdx = size(C,1);
end
[m, n] = size(C);
end
supply = supply(:)';
demand = demand(:)';
if length(supply) ~= m || length(demand) ~= n
error('尺寸不匹配: supply 长度应为 %d, demand 长度应为 %d', m, n);
end
% 初始可行解
switch lower(method)
case 'vogel'
X = vogel_init(C, supply, demand);
otherwise
X = northwest_init(supply, demand);
end
if verbose
fprintf('Initial X matrix:\n');
disp(X);
end
% 验证初始可行解是否满足供需
rowS = sum(X,2)'; colS = sum(X,1);
if verbose
fprintf('Initial X sum rows: [%s]\n', num2str(rowS));
fprintf('Initial X sum cols: [%s]\n', num2str(colS));
fprintf('Expected supply: [%s]\n', num2str(supply));
fprintf('Expected demand: [%s]\n', num2str(demand));
end
if any(abs(rowS - supply) > 1e-6) || any(abs(colS - demand) > 1e-6)
warning('bszy:initial','初始可行解不平衡，尝试使用 Northwest corner 生成初始解');
X = northwest_init(supply, demand);
if verbose
fprintf('Recomputed X (northwest) sum rows: [%s]\n', num2str(sum(X,2)'));
fprintf('Recomputed X (northwest) sum cols: [%s]\n', num2str(sum(X,1)));
end
end
maxIter = 100;
tol = 1e-8;
iter = 0;
while iter < maxIter
iter = iter + 1;
% 基变量掩码（正分配视为基）
basic = X > tol;
% 处理退化：确保基数为 m+n-1
basic = add_degeneracy(basic, m, n, C);
% 求解 u 和 v (势)
[u, v] = solve_potentials(C, basic);
% 计算改进成本（检验最优性）
U = repmat(u(:), 1, n);
V = repmat(v(:)', m, 1);
delta = C - (U + V);
% 将基变量位置设为 +Inf，确保进入变量来自非基单元
temp = delta;
temp(~isfinite(temp)) = Inf;
temp(basic) = Inf;
if verbose
fprintf('delta matrix (Inf means basic or invalid):\n');
disp(temp);
fprintf('basic mask:\n');
disp(basic);
end
% 如果所有非基 delta >= 0 则最优
if all(delta(~basic) >= -tol)
break;
end
% 选择最负的改进成本作为进入变量（只看非基单元）
[minVal, idx] = min(temp(:));
[rEnter, cEnter] = ind2sub([m, n], idx);
if verbose
fprintf('iter=%d minDelta=%.6f enter=(%d,%d)\n', iter, minVal, rEnter, cEnter);
end
if basic(rEnter, cEnter)
error('选择了基变量作为进入变量 — 逻辑错误，详细 delta 已打印');
end
% 在基矩阵中加入进入变量并寻找回路
basic2 = basic;
basic2(rEnter, cEnter) = true;
cycle = find_cycle(basic2, [rEnter, cEnter]);
if isempty(cycle)
error('未找到闭合回路');
end
% cycle 是一系列坐标，最后一项等于第一项。交替 + -
coords = cycle(1:end-1, :);
signs = 1 - 2*((1:size(coords,1))' - 1); % 1, -1, 1, -1,... but we'll use simpler
% 更明确的符号：第一个为 +，第二个为 -，交替
signs = zeros(size(coords,1),1);
for k = 1:length(signs)
if mod(k-1,2) == 0
signs(k) = 1;
else
signs(k) = -1;
end
end
% 找到所有带负号的位置，并计算可调整的最小量
minusIdx = find(signs == -1);
theta = Inf;
for k = minusIdx'
r = coords(k,1); c = coords(k,2);
theta = min(theta, X(r,c));
end
if isinf(theta)
error('调整步长无穷大，异常');
end
if verbose
fprintf('theta=%.6f negPositions=%s\n', theta, mat2str(coords(minusIdx,:)));
end
% 沿回路调整分配
for k = 1:size(coords,1)
r = coords(k,1); c = coords(k,2);
X(r,c) = X(r,c) + signs(k)*theta;
end
% 清除接近零的值
X(abs(X) < tol) = 0;
if verbose
fprintf('after adjust, X nonzeros=%d total=%g\n', nnz(X), sum(X(:).*origC(:))/1);
end
end
% 如果加入了虚拟行/列，去除对应列/行后用原始成本计算总费用
if strcmp(dummyType, 'none')
totalCost = sum(sum(X .* C));
else
if strcmp(dummyType, 'col')
% 原始列数为 orig_n
X_trim = X(:,1:orig_n);
totalCost = sum(sum(X_trim .* origC));
else
X_trim = X(1:orig_m,:);
totalCost = sum(sum(X_trim .* origC));
end
end
end
function X = northwest_init(supply, demand)
% Northwest corner method
m = length(supply); n = length(demand);
X = zeros(m, n);
i = 1; j = 1;
s = supply; d = demand;
while i <= m && j <= n
q = min(s(i), d(j));
X(i,j) = q;
s(i) = s(i) - q;
d(j) = d(j) - q;
if abs(s(i)) < 1e-12
i = i + 1;
elseif abs(d(j)) < 1e-12
j = j + 1;
end
end
end
function X = vogel_init(C, supply, demand)
% Vogel's approximation method
s = supply; d = demand;
m = length(s); n = length(d);
X = zeros(m, n);
tolSmall = 1e-12;
% 每次迭代根据剩余供需重算可用行列，避免死循环
while any(s > tolSmall) && any(d > tolSmall)
rows = s > tolSmall;
cols = d > tolSmall;
% 计算行差与列差
rowPenalty = inf(1,m);
for i = 1:m
if ~rows(i), continue; end
costs = C(i,cols);
if isempty(costs)
continue;
end
if numel(costs) >= 2
cs = sort(costs);
rowPenalty(i) = cs(2) - cs(1);
else
rowPenalty(i) = costs(1);
end
end
colPenalty = inf(1,n);
for j = 1:n
if ~cols(j), continue; end
costs = C(rows,j);
if isempty(costs)
continue;
end
if numel(costs) >= 2
cs = sort(costs);
colPenalty(j) = cs(2) - cs(1);
else
colPenalty(j) = costs(1);
end
end
% 选择最大惩罚
[rmax, ridx] = max(rowPenalty);
[cmax, cidx] = max(colPenalty);
if isinf(rmax) && isinf(cmax)
% 没有可选行或列，退出以防死循环
break;
end
if rmax >= cmax
i = ridx;
% 在可用列中选择最小成本
availableCols = find(cols);
[~, k] = min(C(i,availableCols));
j = availableCols(k);
else
j = cidx;
availableRows = find(rows);
[~, k] = min(C(availableRows,j));
i = availableRows(k);
end
q = min(s(i), d(j));
if q <= tolSmall
% 无法分配更多，避免死循环
break;
end
X(i,j) = X(i,j) + q;
s(i) = s(i) - q;
d(j) = d(j) - q;
end
% 如果还有剩余供需（可能是由于惩罚计算或边界情况），在剩余行列上使用 Northwest 补齐
rows_rem = find(s > tolSmall);
cols_rem = find(d > tolSmall);
if ~isempty(rows_rem) && ~isempty(cols_rem)
s_sub = s(rows_rem);
d_sub = d(cols_rem);
% 使用 Northwest corner 在子问题上分配
Y = northwest_init(s_sub, d_sub);
for ii = 1:length(rows_rem)
for jj = 1:length(cols_rem)
X(rows_rem(ii), cols_rem(jj)) = X(rows_rem(ii), cols_rem(jj)) + Y(ii,jj);
end
end
end
end
function basic = add_degeneracy(basic, m, n, C)
% 如果基数不足 m+n-1，则加入0分配以避免退化
need = m + n - 1 - nnz(basic);
if need <= 0
return;
end
% 选择成本最小的非基单元逐个加入
mask = ~basic;
costList = C;
costList(~mask) = Inf;
costList = costList(:);
[~, order] = sort(costList);
added = 0;
for k = 1:length(order)
if costList(order(k)) == Inf
break;
end
r = ceil(order(k)/n);
c = mod(order(k)-1, n) + 1;
if ~basic(r,c)
basic(r,c) = true;
added = added + 1;
if added >= need
break;
end
end
end
end
function [u, v] = solve_potentials(C, basic)
% 使用基于传播的求解方法：从 u1=0 出发，通过基变量 u_i+v_j=C(i,j) 传播
[m, n] = size(C);
u = nan(m,1);
v = nan(n,1);
u(1) = 0;
% 构建基变量索引列表
[Rs, Cs] = find(basic);
changed = true;
iter = 0;
while changed && iter < (m+n)
changed = false;
iter = iter + 1;
for k = 1:length(Rs)
i = Rs(k); j = Cs(k);
cij = C(i,j);
if ~isnan(u(i)) && isnan(v(j))
v(j) = cij - u(i);
changed = true;
elseif ~isnan(v(j)) && isnan(u(i))
u(i) = cij - v(j);
changed = true;
end
end
end
% 若仍有未被确定的势，设为0（不会影响改进成本的相对比较）
u(isnan(u)) = 0;
v(isnan(v)) = 0;
end
function cycle = find_cycle(basicMask, start)
% 在 basicMask 中寻找包含 start 的一个闭合回路（基于二部图的 BFS）
[m, n] = size(basicMask);
r0 = start(1); c0 = start(2);
% 构建邻接表：节点 1..m 为行，m+1..m+n 为列
Nnodes = m + n;
adj = cell(Nnodes,1);
for i = 1:m
for j = 1:n
if basicMask(i,j)
u = i; v = m + j;
adj{u}(end+1) = v;
adj{v}(end+1) = u;
end
end
end
u0 = r0; v0 = m + c0;
% 临时移除进入边 (u0,v0)，在其余基边中寻找从 u0 到 v0 的路径
adj{u0}(adj{u0} == v0) = [];
adj{v0}(adj{v0} == u0) = [];
% BFS 查找最短路径
prev = zeros(Nnodes,1);
visited = false(Nnodes,1);
queue = zeros(Nnodes,1);
head = 1; tail = 1;
queue(tail) = u0; visited(u0) = true;
found = false;
while head <= tail
cur = queue(head); head = head + 1;
for nb = adj{cur}
if ~visited(nb)
visited(nb) = true;
prev(nb) = cur;
tail = tail + 1; queue(tail) = nb;
if nb == v0
found = true; break;
end
end
end
if found, break; end
end
if ~found
cycle = [];
return;
end
% 重构节点路径，从 u0 到 v0
path_nodes = v0;
cur = v0;
while cur ~= 0 && cur ~= u0
cur = prev(cur);
path_nodes = [cur; path_nodes];
end
% 将节点路径转换为格子坐标（相邻节点对应一个基本单元）
coords = [];
for t = 1:(length(path_nodes)-1)
a = path_nodes(t); b = path_nodes(t+1);
if a <= m && b > m
coords = [coords; a, b - m];
elseif a > m && b <= m
coords = [coords; b, a - m];
end
end
% 在路径末尾加上进入单元，形成闭合回路的最后一条边
coords = [coords; r0, c0];
% 将 coords 旋转，使第一个元素为进入单元 (r0,c0)
idx = find(coords(:,1) == r0 & coords(:,2) == c0, 1);
if isempty(idx)
cycle = [];
return;
end
coords = [coords(idx:end, :); coords(1:idx, :)];
% 确保首尾闭合
cycle = [coords; coords(1,:)];
end
结果如下：
>> [X, totalCost, iter] = bszy([30;50;20], [60;40], [2 3;1 4;5 2])
X =
10 20
50 0
0 20
totalCost = 170
iter = 1
