简介：二阶段法与单纯形法的MATLAB实现代码包，专为学习线性规划算法的学生、科研人员及工程师设计，重点解决初始解不可行情况下如何求出可行解并进一步寻优的问题。代码包共3个文件，主程序以.m脚本实现完整的两阶段法流程，配套.asv文件为编辑过程中生成的自动备份，另有.txt说明文档帮助理解算法步骤与函数调用方式。整个包仅8KB，结构精炼，适合直接阅读和调试。内容涉及如何定义线性规划目标函数与约束、如何通过第一阶段构造人工变量得到初始可行解、如何转入第二阶段调用单纯形法迭代求最优值，并涵盖无界解、无解等异常情况的处理思路。目前已有374人学习浏览，对于希望借助MATLAB快速掌握两阶段法和单纯形法原理的读者而言，这是一份可直接运行、便于对照算法流程的入门与进阶参考资料。
1. 二阶段法matlab：单纯形表里没有单位列时，才是真正入口
手写单纯形法时，教科书默认你手里有一张现成的单纯形表——单位列排好，初始可行基直接可见。换成 matlab 自己搭模型就会发现，"≤" 约束补出的松弛变量确实带来单位列，"≥" 约束补出的剩余变量系数是 −1，等式约束什么也补不出来，表的左侧根本没有单位矩阵。二阶段法就是为缺基问题设计的标准解法：第一阶段用人工变量和辅助目标强行制造可行基，第二阶段卸掉人工变量、恢复原目标继续迭代。下面按原理、实现、参数、排错四个层次展开，适合运筹学课程、算法项目以及想绕开 linprog 黑盒的工程师。
2. 二阶段单纯形法原理：人工变量、辅助目标与可行基判定
2.1 为什么没有初始可行基：松弛变量只解决一半问题
单纯形法的每次迭代都要求当前基是可行基，也就是基变量取值全部非负。模型若全是 "≤" 约束，加松弛变量后天然得到一个单位矩阵，初始可行基零成本获得。实际建模里 "≥" 约束和等式约束出现频率很高：运输问题里的下限约束、配方问题的比例约束、逻辑约束转换出来的等式，都会让系数矩阵缺列。
以这个最小化问题为例：
min z = x1 + x2 s.t. 2x1 + x2 ≥ 4 x1 + 3x2 ≥ 6 x1, x2 ≥ 0
补剩余变量 x3、x4 后得到标准形：
2x1 + x2 − x3 = 4 x1 + 3x2 − x4 = 6
x3、x4 的系数是 −1，不是 +1 的单位列，两行都没有可用的初始基。单纯形法要求初始基可行，这一步卡住，后续的转轴公式就全部无从谈起。matlab 优化工具箱里的 linprog 默认对偶单纯形能自己处理这种情况，但手写实现时，两阶段法是最直接、数值上也最干净的选择。
2.2 第一阶段：把目标函数换成人工变量之和
给缺基的每一行补一个人工变量，系数取 +1，构造辅助问题：
min w = x5 + x6 2x1 + x2 − x3 + x5 = 4 x1 + 3x2 − x4 + x6 = 6 x1..x6 ≥ 0
基变量取 x5、x6，初始可行基立即可得。第一阶段的目标只有一个判定作用：原问题有可行解，当且仅当最优的 w = 0。w 是人工变量之和，天然非负；w = 0 说明所有人工变量恰好取 0，剩下的 x1..x4 满足原约束，人工列只是第一阶段迭代用的脚手架。
回到例子。初始解 x5 = 4、x6 = 6，w = 10。把 w 写成非基变量的表达式：
w = 10 − 3x1 − 4x2 + x3 + x4
min 问题取检验数最负的变量进基，x2 的检验数 −4 最负，进基。最小比值 min(4/1, 6/3) = 2，x6 出基。转轴后 w = 2 − (5/3)x1 + x3 − (1/3)x4 + (4/3)x6，x1 的检验数 −5/3 最负，进基；最小比值 min(2/(5/3), 6/(1/3)) = 6/5，x5 出基。此时 x1 = 6/5、x2 = 8/5，w = 0，第一阶段结束。
| 迭代 | 进基变量 | 出基变量 | 最小比值说明 | w 值 |
|---|---|---|---|---|
| 1 | x2（检验数 −4） | x6 | min(4/1, 6/3) = 2 | 2 |
| 2 | x1（检验数 −5/3） | x5 | min(2/(5/3), 6/(1/3)) = 6/5 | 0 |
若最优 w > 0，说明至少一个人工变量必须取正值才能满足约束，原问题无可行解，直接返回。这也是二阶段法和大 M 法的主要差别：大 M 法把人工变量乘一个很大的 M 塞进原目标，M 选取不当会让数值问题病态化；二阶段法把可行性判定拆成独立的第一阶段，不引入额外大数，数值行为更可预期。
2.3 第二阶段：丢弃人工列，恢复原目标继续迭代
上例第一阶段结束时基是 {x1, x2}，没有人工变量残留。删掉人工列，目标换回 z = x1 + x2，用同一个基继续迭代。解 B'λ = cB 得到 λ = (0.4, 0.2)，非基变量检验数 r3 = 0.4、r4 = 0.2 全部非负，已经满足 min 问题的最优性判定，不需要再转轴。
这里有个容易被忽略的细节：第一阶段的最优基里可能残留值为 0 的人工变量，原因通常是退化或冗余约束。此时不能直接删人工列，否则基矩阵少一列。标准做法是找该行一个非零且不在基内的原变量列做一次退化转轴，把人工变量换出；如果整行系数全是 0，说明该约束是其他行的线性组合，整行删掉。这个逻辑在下一章的代码里单独占一个步骤，也是两阶段法手写实现里报错率最高的地方。
注意：第二阶段换回原目标后，检验数必须按新目标从头重算。沿用第一阶段残留的检验数是最常见的错误，因为对偶变量 λ 依赖基变量的目标系数 cB，目标一变，全部判别数失效。
3. matlab两阶段法完整实现：simplex_core 与人工变量清基
3.1 输入约定：标准形、b ≥ 0 与不等式转换
代码按标准形 min c'x、Ax = b、x ≥ 0 书写。调用前需要自己完成两步：不等式补松弛变量或剩余变量，再整理成等式；右端项正负由函数内部处理。函数对 b < 0 的行自动翻号，调用方少一个出错点。翻号只改变该行的符号，可行域和解集不变。
3.2 核心迭代：simplex_core
第一阶段和第二阶段共用同一个单纯形核心：给定初始可行基，反复做"算检验数、判定最优、算进基方向、最小比值出基"。检验数用对偶变量计算，不显式求逆：
λ 满足 B'λ = cB，检验数 r = c − A'λ。
function [x, z, basis, iters] = simplex_core(A, b, c, basis, tol) % 单纯形法核心迭代：给定初始可行基，求解 min c'*x % 输入 A(m,n)、b(m,1)、c(n,1)、basis(m,1) 基变量列号、tol 容差 % 输出 x(n,1) 最优解、z 最优值、basis 最终基、iters 实际转轴次数 [m, n] = size(A); iters = 0; while true B = A(:, basis); xB = B \ b; % 基变量取值，同时是可行性判定依据 lambda = B' \ c(basis); % 解 B'*lambda = cB，得到对偶变量 r = c - A' * lambda; % 所有列的检验数 reduced cost r(basis) = 0; % 数值上强制基变量检验数为 0 [rmin, q] = min(r); if rmin >= -tol % min 问题：检验数全非负即最优 x = zeros(n, 1); x(basis) = xB; z = c' * x; return; end d = B \ A(:, q); % 进基方向 if max(d) <= tol error('two_phase:unbounded', '检验数为负但进基方向无正分量，问题无界'); end ratio = xB ./ d; % 最小比值法则 ratio(d <= tol) = Inf; % 非正方向不参与比值 [~, p] = min(ratio); % p 行出基，q 列进基 basis(p) = q; iters = iters + 1; if iters > 10000 error('two_phase:maxiter', '迭代超过 10000 次，疑似退化循环'); end end end几个关键点说明：
B \ b和B' \ c(basis)走的是 MATLAB 的 LU 分解求解路径，不要写成inv(B) * b。显式求逆在维数超过几十、条件数偏大的时候，误差会直接污染后续转轴。r(basis) = 0是必要的数值修正。浮点运算下基变量检验数会残留 1e-16 量级的噪声，不置零可能把最优基误判为仍有负检验数，导致多余的转轴甚至循环。- 无界判定依据是：存在负检验数，但该方向 d 的分量全部小于等于 0，增大进基变量不会碰到任何约束上界。这里用
max(d) <= tol而不是all(d <= 0)，容差同样施加在这一步。
3.3 主函数：两阶段衔接与人工变量清基
主函数把"补人工变量、第一阶段、清基、第二阶段"串起来。第一阶段的目标向量只在人工变量列位置置 1，其余位置为 0。
function [x, z, info] = two_phase_simplex(A, b, c, tol) % 二阶段法 matlab 实现：min c'*x，s.t. A*x = b，x >= 0 % 输入 A 等式约束系数 m*n，b 右端项 m*1（负分量自动翻号），c 目标系数 n*1 % 输出 x 最优解（无可行解时为空），z 最优值（不可行时为 Inf） % info.status = 'optimal' / 'infeasible' % info.phase1_iters / info.phase2_iters 记录两阶段转轴次数 if nargin < 4, tol = 1e-9; end [m, n] = size(A); info = struct('phase1_iters', 0, 'phase2_iters', 0, 'status', 'optimal'); % ---- 1. b 负行翻号 ---- neg = find(b < 0); A(neg, :) = -A(neg, :); b(neg) = -b(neg); % ---- 2. 找天然单位列 ---- basis = zeros(m, 1); for i = 1:m for j = 1:n if abs(A(i,j) - 1) < tol && ... abs(sum(abs(A(:,j))) - 1) < tol && ... ~ismember(j, basis) basis(i) = j; break; end end end % ---- 3. 给缺基行补人工变量，构造第一阶段矩阵 ---- artRows = find(basis == 0); k = numel(artRows); A1 = [A, zeros(m, k)]; for t = 1:k A1(artRows(t), n + t) = 1; basis(artRows(t)) = n + t; end c1 = [zeros(n, 1); ones(k, 1)]; % 第一阶段目标：人工变量之和 % ---- 4. 第一阶段迭代 ---- [~, w, basis1, it1] = simplex_core(A1, b, c1, basis, tol); info.phase1_iters = it1; if w > tol % 辅助目标没压到 0 => 原问题不可行 info.status = 'infeasible'; x = []; z = Inf; return; end % ---- 5. 清掉基里残留的人工变量（值为 0 的退化情形） ---- artCols = n + (1:k); delRows = []; for t = find(ismember(basis1, artCols))' % 找该行非零、且不在基内的原变量列 cand = find(abs(A1(t, 1:n)) > tol & ~ismember(1:n, basis1)); if isempty(cand) delRows = [delRows, t]; % 该行线性相关，冗余约束，整行删掉 else basis1(t) = cand(1); % 退化转轴：人工变量出基 end end if ~isempty(delRows) A1(delRows, :) = []; b(delRows) = []; basis1(delRows) = []; end % ---- 6. 第二阶段：丢人工列，换回原目标 ---- A2 = A1(:, 1:n); [x, z, ~, it2] = simplex_core(A2, b, c, basis1, tol); info.phase2_iters = it2; end3.4 三个容易写错的细节
第一，第 2 步用sum(abs(A(:,j))) ≈ 1判单位列，等价于"该列只有一个非零元素且模为 1"，配合A(i,j) ≈ 1就锁定了那一个非零元素所在的行。比逐行索引做差集更简洁，也避免 m = 1 时索引向量为空导致&&运算报错。
第二，第 5 步必须先收集delRows再统一删除，不能边遍历边删行。边删边遍历会让后续行号错位，这是手写两阶段法最常见的隐性 bug，症状是"第一次运行正常，第二次结果完全不对"。
第三，第 5 步对 tol 的依赖方式和第 2 步不一样：第 2 步要求列近似单位列，第 5 步只要求系数非零。如果把 tol 调大，第 5 步可能把微小非零系数当成零，误判冗余约束而删行——这个错误比第 2 步的误判隐蔽得多。所以默认 tol 维持 1e-9 不要轻易动。
4. 二阶段法求解算例与参数调整：tol、迭代上限与 Bland 法则
4.1 跑通最小算例，核对输出
把两个函数存成 m 文件，运行测试脚本：
% 例：min z = x1 + x2 % 约束 2x1 + x2 >= 4；x1 + 3x2 >= 6；x1,x2 >= 0 % 标准形：补剩余变量 x3、x4，等价 form 见正文 2.1 节 A = [2 1 -1 0; 1 3 0 -1]; b = [4; 6]; c = [1; 1; 0; 0]; [x, z, info] = two_phase_simplex(A, b, c); fprintf('x = [%g %g %g %g]\n', x); fprintf('z = %g\n', z); fprintf('phase1 pivots = %d, phase2 pivots = %d\n', ... info.phase1_iters, info.phase2_iters);期望输出：x = [1.2 1.6 0 0]，z = 2.8，phase1 pivots = 2，phase2 pivots = 0。第二阶段零转轴是正常现象：第一阶段结束时的基有时恰好让原目标的检验数全部非负。但绝不能因此省略第二阶段——换目标后检验数必须从头重算，这一条是两阶段法正确性的底线。验证时重点看 w 是否为 0、x 是否满足 A*x = b，这两点比目标值更先暴露问题。
4.2 三个必调参数：容差、迭代上限、Bland 法则
手写单纯形法要调的参数不多，但每个都能让结果从"对"变成"错"。
| 参数 | 代码位置 | 作用 | 建议值 |
|---|---|---|---|
| tol | two_phase_simplex 第 4 个入参 | 检验数、最小比值、可行性判定共用 | 1e-9 |
| 迭代上限 | simplex_core 底部iters > 10000 | 防止退化循环死转 | 10000 |
| Bland 开关 | 进基/出基选择处手动替换 | 消除退化循环 | 出现循环时开启 |
tol 太大（如 1e-6）会吃掉真实的负检验数，输出一个接近最优但不是最优的解；tol 太小（如 1e-12）会把浮点噪声当成合法转轴方向，在最优基附近反复横跳。1e-9 对绝大多数双精度算例够用，碰到量级差异极大的矩阵，先做行缩放而不是改容差。
退化循环的现场特征是：某次转轴后目标值不再下降，但检验数里仍有负值。Bland 法则（最小下标规则）是教科书上保证终止的进基/出基规则——在负检验数里选下标最小的进基，在最小比值并列的行里选基变量列号最小的出基。替换 simplex_core 中的两处选择：
% 替换进基选择：负检验数里取最小下标 negIdx = find(r < -tol); q = negIdx(1); % 替换出基选择：比值并列时取基变量列号最小 [minVal, ~] = min(ratio); ties = find(abs(ratio - minVal) < tol); [~, idx] = min(basis(ties)); p = ties(idx);注意 Bland 法则在浮点实现里必须基于 tol 过滤后的负检验数列表，否则 1e-16 的负噪声会破坏"最小下标"的语义。代码里保留 10000 次上限，是因为很多退化实例靠随机抖动能自己脱离，上限只当保险丝用。
4.3 随机算例交叉验证：与 linprog 逐个比对
验证两阶段法实现没有写错，最可靠的办法是随机生成一批已知可行的问题，和 linprog 对比最优值：
rng(42); bad = 0; for t = 1:50 m = 4; n = 8; A = randn(m, n); x0 = rand(n, 1) + 0.1; % 随机正解，保证原问题可行 b = A * x0; c = rand(n, 1) + 0.5; % 目标系数取正，配合 x>=0 保证有下界 [x1, z1, info] = two_phase_simplex(A, b, c); opt = optimoptions('linprog', 'Display', 'off'); [~, z2] = linprog(c, [], [], A, b, zeros(n, 1), [], [], opt); if ~strcmp(info.status, 'optimal') || abs(z1 - z2) > 1e-6 bad = bad + 1; fprintf('第 %d 组不一致：self=%g linprog=%g\n', t, z1, z2); end end fprintf('不一致实例数：%d / 50\n', bad);这段验证的要点：b 用 A*x0 生成，保证可行域非空，第一阶段必然收敛到 w = 0；约束全是等式，正好覆盖二阶段法最需要人工变量的场景。目标系数取正，配合 x ≥ 0 保证最优值有下界，避免随机出无界实例干扰比对。若出现不一致，先打印该组的 A、b、c 手工检查，多数问题出在容差设置而不是算法逻辑。
提示：simplex_core 返回的 iters 只统计实际转轴次数，最后一次纯最优性校验不计数。调试时如果 phase2_iters 为 0，不代表第二阶段没执行，只是没有发生转轴。
5. 二阶段法排错：无界、不可行与人工变量残留的定位
5.1 三类报错的快速定位
无界错误来自 simplex_core 抛出的two_phase:unbounded。先查模型抄写：等式约束的符号、决策变量非负限制有没有遗漏。无界多半是少写了一个约束，而不是代码问题。第二阶段报无界而第一阶段正常，重点检查原目标 c 的符号是否写反。
不可行状态来自第一阶段 w > tol。先把 tol 调到 1e-10 重跑一次，排除容差误判；仍然不可行再用 linprog 跑同一模型复核。矛盾约束如 x1 ≥ 5 和 x1 ≤ 3 同时出现是最常见的来源，这类问题第一阶段结束时的 w 会明显大于 0，不是贴边的那种。
5.2 收敛断言放在 return 之前
在 simplex_core 的 return 前加两行断言，把"解不可行"和"目标值不对"两类错误分开：
resid = norm(A * x - b) / (1 + norm(b)); assert(resid < 1e-8, '可行性断言失败：残差 %g', resid); assert(all(x > -1e-9), '可行性断言失败：存在负分量');残差过大说明基变量取值和约束对不上，问题出在转轴过程；只有负分量说明解的符号出了问题，多半是目标方向或翻号逻辑写反。这两条断言在任何规模的算例上都适用，是比对照最优值更前置的检查。
5.3 规模变大时的两个改动
数据规模上去之后，两个位置最先遇到瓶颈。一是simplex_core里每次转轴都调用三次B \ b、B' \ c(basis)、B \ A(:,q)，应该改为在一次[L,U] = lu(B)后复用分解结果，每次转轴只做回代。二是现在维护的是整张表格式单纯形，稀疏大规模问题应当改成修正单纯形（revised simplex），只存基矩阵，检验数按列逐个计算而不是更新整张表。两阶段法的结构不变，变的只是线性方程组的求解方式和检验数的计算粒度。
本文还有配套的精品资源，点击获取
