运筹学之单纯形法（超详细讲解 + MATLAB可运行代码实现）

运筹学之单纯形法（超详细讲解 + MATLAB可运行代码实现）

引用

CSDN

1.

https://m.blog.csdn.net/2301_79571587/article/details/145936826

单纯形法是解决线性规划问题的重要算法，它通过迭代的方式在可行域的顶点间寻找最优解。本文将从单纯形法的定义出发，详细讲解其基本原理，并通过MATLAB代码实现，帮助读者全面理解这一算法。

1. 单纯形法的定义

单纯形法是一种基于几何顶点遍历的线性规划求解算法，其核心逻辑可概括为：

1. 顶点择优原理：

• 线性规划的可行域构成凸多面体，最优解必出现在顶点（基可行解）上。
• 通过有限步迭代，从一个顶点沿目标函数改进方向移动到相邻顶点。

2. 方向选择机制：

• 利用检验数（Reduced Cost）判断当前解的改进潜力。
• 若存在正检验数（最大化问题）或负检验数（最小化问题），则目标函数可继续优化。

3. 终止条件：

• 所有检验数均不满足优化方向时达到最优解。
• 若发现无界方向（存在无限改进空间），则问题无界。

2. 单纯形法的基本原理

2.1 核心步骤

1. 将线性规划的标准型化成线性规划的规范型，来获得一个初始可行解。
2. 对初始基可行解最优性判别，若最优，则停止；否则转下一步。
3. 从初始基可行解向相邻的基可行解转换，使得目标值有所改善，重复第二步和第三步直至找到最优解。

2.2 步骤拆解

为更清楚地了解并操作单纯形法，必须要解决三个问题：

2.2.1 如何确定初始的基可行解？

首先我们先了解四个基本概念：

• 规范型：约束条件矩阵A中有一个单位矩阵I。
• 基：约束系数矩阵的A的秩为m，若A中的B是m*m的可逆矩阵，则称B为线性规划问题的一个基。
• 基解：对于基B，是对应于基B的决策变量（基变量），系数矩阵A中的其它列向量对应的决策变量称为非基变量。如果令非基变量为0，则称为基解。
• 基可行解：非负条件的基解称为基可行解。

认识了四个基本概念后，我们先将线性规划标准型为规范型，其中，系数矩阵中含有一个单位矩阵 I ,不妨设单位矩阵为:

简单示例如下：

找到初始可行基，令非基变量取值为0，便得到一组基可行解。上述示例的基可行解如下：

2.2.2 如何进行解的最优性判别？

首先我们先了解一个概念：

• 检验数：表示将非基变量增加一个单位时，目标函数值的变化率。

其计算公式为：

• ：非基变量在目标函数中的系数
• ：基变量在目标函数中的系数向量
• ：基矩阵的逆矩阵
• ：非基变量在约束中的系数列向量

实例演算：

步骤1：计算初始检验数

实例演算：

步骤2：判别最优性

判别规则：

• 最大化问题（Max）：若所有非基变量的≤ 0，当前解为最优解。
• 最小化问题（Min）：若所有非基变量的≥ 0，当前解为最优解。

实例演算：

所有非基变量的检验数均（最大化问题），说明当前解 Z = 0 非最优，需迭代改进。

2.2.3 如何寻找改进的基可行解？

首先我们了解两个概念：

• 出基变量：从当前基变量中被移除的变量，其离基是为保证新解的可行性（非负约束）
• 入基变量：从非基变量中被选中加入基的变量，其引入将直接改进目标函数值

步骤一：选择入基变量--选择最大正检验数对应的变量

实例演算：

步骤二：选择出基变量，确保新解满足非负约束，通过 最小比值测试（Minimum Ratio Test）选择出基变量，其公式为：

：当前基变量取值
：入基变量在约束中的系数

实例演算：

步骤三：基变换与解更新

操作：

1. 构造新基矩阵：将入基变量的列替换出基变量对应的列
2. 高斯消元：通过行变换将新基变量对应的列转换为单位向量
3. 更新解：令非基变量为0，直接读取基变量取值

实例演算：

第一次操作完以上三步后，重复第二步和第三步直至找到最优解。

3. 用MATLAB实现单纯形法

3.1 单纯形法算法代码

function [xm, fm, exitflag, iterations] = enhanced_simplex(c, A, b, max_iter)
% ENHANCED_SIMPLEX 改进版单纯形法 
% 功能：解决标准形线性规划问题 min c'x s.t. Ax = b, x >= 0 
% 新增特性：两阶段法、LU分解、Bland规则 
 
%% ========== 输入参数校验与初始化 ==========
if nargin == 0  % 测试案例 
    c = [-3; -2; 0; 0];
    A = [2, 1, 1, 0; 1, 2, 0, 1];
    b = [4; 3];
    max_iter = 500;
elseif nargin < 4 
    max_iter = 500;
end 
 
[m, n] = size(A);
exitflag = -1; xm = []; fm = 0; iterations = 0;
 
% 维度一致性检查 
if length(b) ~= m || length(c) ~= n 
    error('输入维度不匹配: A是%dx%d, b是%dx1, c是%dx1', m,n,length(b),length(c));
end 
%% ========== 第一阶段：可行性检查 ==========
% 构造人工变量 
phase1_A = [A, eye(m)];
phase1_c = [zeros(n,1); ones(m,1)];
artificial_basis = (n+1 : n+m)';  % 必须转为列向量 
% 调用核心算法 
[x_phase1, ~, exitflag_phase1] = core_simplex(phase1_c, phase1_A, b, artificial_basis, max_iter);
% 判断可行性 
if exitflag_phase1 ~= 1 || any(abs(x_phase1(n+1:end)) > 1e-6)
    disp('无可行解');
    return;
end 
 
%% ========== 第二阶段：优化原目标 ==========
basis = find(x_phase1(1:n) > 1e-6);  % 提取有效基 
if length(basis) ~= m 
    error('基变量数量不足，可能遇到退化问题');
end 
[xm, fm, exitflag, iterations] = core_simplex(c, A, b, basis, max_iter);
end  % 主函数结束 
 
%% ========== 核心算法子函数 ==========
function [x_opt, f_opt, exitflag, iter] = core_simplex(c, A, b, basis, max_iter)
% CORE_SIMPLEX 核心单纯形算法实现 
%% 初始化 
[m, n] = size(A);
tol = 1e-10;
x_opt = zeros(n,1);
exitflag = -1;
iter = 0;
% 初始基校验 
A_b = A(:,basis);
if rank(A_b) < m 
    error('初始基矩阵秩不足');
end 
[L, U, P] = lu(A_b);  % LU分解 
%% 主循环 
for iter = 1:max_iter 
    % 计算检验数 
    lambda = L' \ (U' \ c(basis));
    reduced_cost = c' - lambda' * A;
    % 最优性检查 
    if all(reduced_cost >= -tol)
        x_b = U \ (L \ (P*b));
        x_opt(basis) = x_b;
        f_opt = c' * x_opt;
        exitflag = 1;
        return;
    end 
    % Bland规则选择入基变量 
    s = find(reduced_cost < -tol, 1, 'first');
    % 计算方向向量 
    As = A(:,s);
    direction = U \ (L \ (P*As));
    % 无界判断 
    if all(direction <= tol)
        exitflag = 0;
        return;
    end 
    % 计算离基变量 
    x_b = U \ (L \ (P*b));
    theta = x_b ./ direction;
    theta(direction <= tol) = Inf;
    [~, q] = min(theta);
    % 更新基变量 
    basis(q) = s;
    % 更新LU分解 
    A_b = A(:,basis);
    [L, U, P] = lu(A_b);
end 
exitflag = -1;  % 超过迭代限制 
end  % 子函数结束   

3.2 测试案例

% ========== 测试案例执行 ==========
c = [-3; -2; 0; 0];
A = [2, 1, 1, 0; 
     1, 2, 0, 1];
b = [4; 3];
[x, fval, flag] = enhanced_simplex(c, A, b);
 
% ========== 输出优化 ==========
disp('════════════════ 单纯形法求解结果 ════════════════');
fprintf('算法状态: %s\n', getExitStatus(flag));  % 状态解释 
fprintf('迭代次数: %d\n', length(x));          % 显示实际迭代次数 
disp('--------------------------------------------');
 
% 解向量分列显示（含变量类型标注）
var_types = {'x₁', 'x₂', 's₁', 's₂'};  % s表示松弛变量 
disp('[最优解]');
for i = 1:length(x)
    fprintf('%s = %8.4f\n', var_types{i}, x(i));
end 
 
fprintf('\n目标函数值: %.4f → 原始问题等效为最大利润: %.2f\n',...
        fval, -fval);  % 转换为最大化问题解释 
disp('══════════════════════════════════════════════');
 
% ========== 辅助函数：状态码转文字 ==========
function status = getExitStatus(flag)
    switch flag 
        case 1,  status = '找到最优解 ✅';
        case 0,  status = '问题无界解 ⚠️ (目标值趋向无穷)';
        case -1, status = '无可行解 ❌';
        otherwise, status = '未知状态';
    end 
end   

输出的结果：

**════════════════ 单纯形法求解结果 ════════════════
算法状态: 找到最优解
迭代次数: 4

[最优解]
x₁ = 1.6667
x₂ = 0.6667
s₁ = 0.0000
s₂ = 0.0000**

目标函数值: -6.3333 → 原始问题等效为最大利润: 6.33
══════════════════════════════════════════════