企业官网建设流程全解析

简介：本资源是一份面向算法学习者与C++初学者的线性规划实践项目，聚焦单纯形法这一经典优化算法的完整代码实现，适用于运筹学、数值计算或算法课程设计等场景。压缩包共3个文件（2个cpp源文件+1个h头文件），总大小仅2KB，结构精简：matrix.h封装矩阵运算基础，matrix.cpp实现单纯形表核心逻辑，main.cpp提供标准线性规划问题的输入接口与求解流程，便于理解算法迭代机制与C++面向过程实现细节。已有741人学习下载，适合希望从零掌握单纯形法原理、动手调试算法步骤、并积累小型数值计算项目经验的学习者。代码注释清晰、模块职责明确，可直接编译运行，亦可作为教学示例拓展为带GUI界面或支持不等式约束扩展的进阶版本。

1. C++实现的单纯形算法计算程序：为什么线性规划求解器不能只靠现成库，而要亲手写一遍核心迭代逻辑？

某高校运筹学课程设计里，学生交上来的“线性规划求解器”八成调用GLPK或CBC的 C 接口，输入.lp文件、输出最优值——但当老师追问“第3次迭代时基变量怎么换入换出？检验数向量 λ 是怎么更新的？”时，90% 的人卡壳。这不是考试刁难，而是暴露了一个实操断层：能跑通 ≠ 理解收敛路径，封装越厚，黑匣子越深。C++实现的单纯形算法计算程序，本质不是为了替代商业求解器，而是把“约束矩阵 A、资源向量 b、目标系数 c 如何在每一轮 pivot 中被原地改写”这件事，变成可打断、可打印、可单步调试的内存状态。它解决的是三类真实需求：教学场景中让学生亲眼看见基可行解如何沿顶点移动；嵌入式或资源受限环境（如某实验室边缘设备）下无法部署大型依赖；以及算法工程师验证自定义预处理（如列生成、对偶热启动）与标准单纯形内核耦合时的数值稳定性。本文不讲凸集理论推导，只聚焦一个干净、无第三方依赖、支持等式/不等式混合约束的 C++17 实现——从初始化基矩阵 B 的 LU 分解开始，到最终输出迭代日志、最优解向量和对偶变量，每一步都对应教科书公式，每一行代码都能在 gdb 里 inspect 变量。你不需要数学系背景，但得愿意打开终端编译、修改main.cpp里的约束系数、然后看控制台逐行吐出 pivot 行索引、换入变量下标、当前目标值变化量。

2. 从零构建标准型：约束预处理与初始可行基的构造策略

单纯形法要求问题严格处于标准型：max cᵀx, s.t. Ax = b, x ≥ 0，其中 A ∈ ℝ^(m×n), rank(A) = m, b ≥ 0。但实际输入往往是混合约束：≤、≥、= 夹杂，变量无非负限制。若跳过这步直接套公式，程序会在第一次迭代就因基矩阵奇异或 b 向量含负值而崩溃。常见做法是分两阶段处理：第一阶段引入人工变量求初始可行解，第二阶段用原目标函数优化。但本实现采用更紧凑的“大M法”统一框架，并在预处理阶段完成三件事：约束标准化、变量替换、基矩阵初始化。下面代码块展示核心预处理逻辑。

2.1 约束标准化：将任意线性约束转为等式 + 非负变量

// constraint_preprocessor.h struct StandardizedProblem { std::vector<double> c; // 标准化后目标系数（长度 n_new） Matrix A; // 约束矩阵（m × n_new） std::vector<double> b; // 资源向量（长度 m），保证 b[i] >= 0 std::vector<int> original_var_index; // 映射：标准化变量索引 → 原变量索引（用于结果还原） int n_original_vars; // 原始变量数（不含松弛/人工变量） }; StandardizedProblem preprocess_constraints( const std::vector<std::vector<double>>& A_orig, const std::vector<double>& b_orig, const std::vector<double>& c_orig, const std::vector<ConstraintType>& types, // {LE, GE, EQ} const std::vector<bool>& is_unbounded) // 原变量是否无下界 ) { size_t m = A_orig.size(); size_t n_orig = c_orig.size(); std::vector<double> c_ext = c_orig; // 扩展目标系数 Matrix A_ext; std::vector<double> b_ext = b_orig; std::vector<int> var_map; // 步骤1：处理每个约束 for (size_t i = 0; i < m; ++i) { const auto& row = A_orig[i]; double rhs = b_orig[i]; ConstraintType ct = types[i]; if (ct == ConstraintType::LE) { // a_i^T x <= b_i → a_i^T x + s_i = b_i, s_i >= 0 // 添加松弛变量：A 新增一列 [0..0,1,0..0]^T（第i行=1） for (size_t j = 0; j < A_ext.rows(); ++j) { A_ext[j].push_back(0.0); } A_ext[i].back() = 1.0; c_ext.push_back(0.0); // 松弛变量目标系数为0 var_map.push_back(-1); // 松弛变量无原始对应 } else if (ct == ConstraintType::GE) { // a_i^T x >= b_i → a_i^T x - s_i = b_i, s_i >= 0 // 添加剩余变量：系数为 -1 for (size_t j = 0; j < A_ext.rows(); ++j) { A_ext[j].push_back(0.0); } A_ext[i].back() = -1.0; c_ext.push_back(0.0); var_map.push_back(-1); } else { // EQ // 无需添加变量，但需确保 rhs >= 0 if (rhs < 0) { // 两边乘-1，翻转A行符号 for (double& a : A_ext[i]) a = -a; rhs = -rhs; } } b_ext[i] = rhs; } // 步骤2：处理无下界变量 x_j ∈ ℝ // 替换为 x_j = x_j' - x_j'', 其中 x_j', x_j'' >= 0 for (size_t j = 0; j < n_orig; ++j) { if (is_unbounded[j]) { // 在c_ext中，原c[j] 拆为 c[j] * x_j' + (-c[j]) * x_j'' double cj = c_ext[j]; c_ext[j] = cj; // x_j' 系数 c_ext.push_back(-cj); // x_j'' 系数 // A_ext 每行第j列复制一份，第j+n_orig列填 -A[i][j] for (size_t i = 0; i < A_ext.rows(); ++i) { double aij = A_ext[i][j]; A_ext[i][j] = aij; A_ext[i].push_back(-aij); } var_map[j] = static_cast<int>(j); var_map.push_back(static_cast<int>(j)); } else { var_map.push_back(static_cast<int>(j)); } } // 步骤3：确保 b >= 0，对负rhs约束添加人工变量（大M法） for (size_t i = 0; i < m; ++i) { if (b_ext[i] < 0) { // 该约束无自然基变量，必须加人工变量 a_i for (size_t j = 0; j < A_ext.rows(); ++j) { A_ext[j].push_back(0.0); } A_ext[i].back() = 1.0; b_ext[i] = -b_ext[i]; // 将 rhs 置正，同时翻转A行符号（已在上一步隐含） // 人工变量目标系数设为 -M（极大化问题中惩罚） c_ext.push_back(-1e8); // M = 1e8，足够大但避免浮点溢出 var_map.push_back(-2); // -2 标记人工变量 } } return {c_ext, A_ext, b_ext, var_map, static_cast<int>(n_orig)}; }

这段代码的关键在于不做假设、不省步骤。例如，GE约束添加的是“剩余变量”（surplus variable），系数为-1，而非某些教程里错误地写成+1；无下界变量拆分时，不仅扩展c向量，还同步复制并取反A的对应列——这是后续基矩阵 LU 分解不崩溃的前提。var_map数组全程记录每个标准化变量的来源，为最终解还原提供依据。注意c_ext.push_back(-1e8)中的M值：太小（如1e3）会导致人工变量未被彻底驱逐出基，太大（如1e12）则引发double精度丢失，使检验数计算失真。这是个需要根据问题规模动态调整的参数，后文避坑章节会详述。

2.2 初始基矩阵 B 的构造与 LU 分解准备

单纯形迭代的核心是快速解B⁻¹b和B⁻¹A_j。若每次迭代都做完整矩阵求逆，时间复杂度升至 O(m³)，不可接受。工业级求解器用 LU 分解维护基矩阵，本实现采用轻量级的 Doolittle 分解（L 下三角、U 上三角、diag(L)=1）。但分解前必须构造一个非奇异且数值稳定的初始基。常见误区是直接取前 m 列作为初始基——当A前 m 列线性相关时，LU 分解立即失败。本实现采用“最大主元搜索 + 人工变量兜底”策略：

// basis_initializer.cpp std::vector<int> find_initial_basis(const Matrix& A, const std::vector<double>& b) { size_t m = A.size(); size_t n = A[0].size(); std::vector<int> basis; basis.reserve(m); // Step 1: 优先选取单位向量列（松弛/人工变量天然构成单位阵） std::vector<bool> used(n, false); for (size_t j = 0; j < n; ++j) { bool is_unit_col = true; int unit_row = -1; for (size_t i = 0; i < m; ++i) { if (std::abs(A[i][j]) > 1e-10) { if (unit_row == -1) { unit_row = i; } else { is_unit_col = false; break; } } } if (is_unit_col && unit_row != -1 && !used[unit_row]) { basis.push_back(j); used[j] = true; used[unit_row] = true; // 标记该行已被覆盖 } } // Step 2: 若基未满，用高斯消元思想选列（最大绝对值主元） while (basis.size() < m) { size_t best_col = 0; double max_pivot = 0.0; for (size_t j = 0; j < n; ++j) { if (used[j]) continue; // 计算该列在未覆盖行中的最大绝对值 double col_max = 0.0; for (size_t i = 0; i < m; ++i) { if (!used[i] && std::abs(A[i][j]) > col_max) { col_max = std::abs(A[i][j]); } } if (col_max > max_pivot) { max_pivot = col_max; best_col = j; } } if (max_pivot < 1e-12) { throw std::runtime_error("Cannot find full-rank initial basis: matrix rank deficient"); } basis.push_back(best_col); used[best_col] = true; } return basis; } // 返回 LU 分解后的 L, U 和行置换向量 P（用于部分主元） std::tuple<Matrix, Matrix, std::vector<size_t>> lu_decompose_with_pivot( const Matrix& A_basis, double eps = 1e-12) { size_t n = A_basis.size(); Matrix L(n, std::vector<double>(n, 0.0)); Matrix U = A_basis; // 初始化 U 为基矩阵副本 std::vector<size_t> P(n); for (size_t i = 0; i < n; ++i) P[i] = i; for (size_t k = 0; k < n; ++k) { // 部分主元：在第k列的k~n-1行找最大绝对值 size_t pivot_row = k; double max_val = std::abs(U[k][k]); for (size_t i = k + 1; i < n; ++i) { if (std::abs(U[i][k]) > max_val) { max_val = std::abs(U[i][k]); pivot_row = i; } } if (max_val < eps) { throw std::runtime_error("LU decomposition failed: singular matrix at step " + std::to_string(k)); } // 行交换 if (pivot_row != k) { std::swap(U[k], U[pivot_row]); std::swap(P[k], P[pivot_row]); } // 更新 L 的第k列（k+1 ~ n-1 行） for (size_t i = k + 1; i < n; ++i) { L[i][k] = U[i][k] / U[k][k]; // 更新 U 的第i行：U[i] -= L[i][k] * U[k] for (size_t j = k; j < n; ++j) { U[i][j] -= L[i][k] * U[k][j]; } } } // 设置 L 对角线为1 for (size_t i = 0; i < n; ++i) L[i][i] = 1.0; return {L, U, P}; }

这里有两个关键设计：一是单位向量列优先。松弛变量、人工变量所在列天然构成单位阵的列，它们是数值最稳定、条件数最优的基候选，应无条件优先入选。二是LU 分解内置部分主元（partial pivoting）。单纯形迭代中，基矩阵可能因连续 pivot 而病态，不加主元的 LU 分解会在第 5~10 次迭代后因U[k][k] ≈ 0而崩溃。lu_decompose_with_pivot返回的置换向量P将在后续B⁻¹b计算中用于重排b向量，这是保证数值鲁棒性的硬性要求。不要试图用std::valarray或 Eigen 的简化接口替代——手写 LU 能让你在U[k][k]极小时立刻throw，而不是让后续B⁻¹A_j计算产生 NaN 再静默失败。

3. 单纯形主循环：pivot 选择、基变换与目标值更新的 C++ 实现

标准单纯形法的主循环只有四步：（1）计算检验数 λ = c_BᵀB⁻¹，（2）选换入变量（最大正检验数），（3）选换出变量（最小比值规则），（4）更新基、更新解、更新目标值。本节代码完全对应这四个数学步骤，无任何抽象封装，所有中间向量（如lambda,db,dxB）均显式声明并打印，便于调试。

3.1 检验数 λ 与换入变量选择：为什么用std::max_element而不用手写循环？

检验数向量 λ 定义为 λᵀ = c_BᵀB⁻¹，其中 c_B 是基变量对应的目标系数子向量。计算 λ 的本质是解线性方程组 Uᵀy = c_Bᵀ，再解 Lᵀλ = y（因为 B = LU ⇒ B⁻¹ = U⁻¹L⁻¹ ⇒ c_BᵀB⁻¹ = c_BᵀU⁻¹L⁻¹，故先解 Lᵀy = c_Bᵀ，再解 Uᵀλ = y）。本实现将此过程拆解为两个前代/后代过程：

// simplex_solver.cpp std::vector<double> compute_reduced_costs( const Matrix& L, const Matrix& U, const std::vector<size_t>& P, const std::vector<double>& c_B, const std::vector<double>& c_nonbasic, const std::vector<int>& basis, const Matrix& A, size_t n_nonbasic) { size_t m = L.size(); // Step 1: 解 L^T * y = c_B^T （前代） std::vector<double> y(m, 0.0); for (size_t i = 0; i < m; ++i) { double sum = 0.0; for (size_t j = 0; j < i; ++j) { sum += L[i][j] * y[j]; } y[i] = c_B[i] - sum; } // Step 2: 解 U^T * lambda = y （后代） std::vector<double> lambda(m, 0.0); for (int i = static_cast<int>(m) - 1; i >= 0; --i) { double sum = 0.0; for (size_t j = i + 1; j < m; ++j) { sum += U[j][i] * lambda[j]; // 注意：U 是上三角，U^T 的第i行是U的第i列 } if (std::abs(U[i][i]) < 1e-12) { throw std::runtime_error("Near-singular U diagonal in reduced cost computation"); } lambda[i] = (y[i] - sum) / U[i][i]; } // Step 3: 计算所有非基变量的检验数 c_j - lambda^T * A_j std::vector<double> reduced_costs(n_nonbasic, 0.0); for (size_t j = 0; j < n_nonbasic; ++j) { double dot = 0.0; for (size_t i = 0; i < m; ++i) { // A_j 是第 j 个非基变量的列向量，需按置换 P 索引 size_t orig_row = P[i]; dot += lambda[i] * A[orig_row][j]; } reduced_costs[j] = c_nonbasic[j] - dot; } return reduced_costs; }

注意reduced_costs[j] = c_nonbasic[j] - dot这一行：它严格对应教材公式z_j - c_j = c_BᵀB⁻¹A_j - c_j（极大化问题中，检验数为c_j - c_BᵀB⁻¹A_j，故此处c_nonbasic[j] - dot即为标准检验数）。std::max_element用于找最大正检验数：

auto max_it = std::max_element(reduced_costs.begin(), reduced_costs.end()); if (*max_it <= 1e-8) { // 所有检验数 ≤ 0，已达最优 return SimplexStatus::OPTIMAL; } int entering_idx = std::distance(reduced_costs.begin(), max_it);

用 STL 算法而非手写循环，不是为了炫技，而是避免for (int i=0; i<n; i++)中i类型与size_t不匹配导致的无符号整数回绕（unsigned integer wraparound）——这是 C++ 新手在n=0时i--陷入死循环的经典翻车点。

3.2 最小比值规则与换出变量选择：如何安全处理退化情形？

换出变量由最小比值规则确定：θ = min{ x_B[i] / (B⁻¹A_enter)[i] | (B⁻¹A_enter)[i] > 0 }。难点在于：（1）B⁻¹A_enter计算涉及 LU 求解，（2）退化时多个x_B[i] = 0，比值为 0，需防除零，（3）浮点误差可能导致本应为 0 的分母计算为1e-16，触发错误换出。本实现采用三重防护：

// Solve B * dB = A_entering_column → get dB = B^{-1} A_enter std::vector<double> solve_B_times_x( const Matrix& L, const Matrix& U, const std::vector<size_t>& P, const std::vector<double>& A_enter_col) { size_t m = L.size(); std::vector<double> b_permuted(m); for (size_t i = 0; i < m; ++i) { b_permuted[i] = A_enter_col[P[i]]; // Apply permutation P to RHS } // Solve L * z = b_permuted (forward substitution) std::vector<double> z(m, 0.0); for (size_t i = 0; i < m; ++i) { double sum = 0.0; for (size_t j = 0; j < i; ++j) { sum += L[i][j] * z[j]; } z[i] = b_permuted[i] - sum; } // Solve U * dB = z (backward substitution) std::vector<double> dB(m, 0.0); for (int i = static_cast<int>(m) - 1; i >= 0; --i) { double sum = 0.0; for (size_t j = i + 1; j < m; ++j) { sum += U[i][j] * dB[j]; } if (std::abs(U[i][i]) < 1e-12) { throw std::runtime_error("Near-zero U diagonal in solve_B_times_x"); } dB[i] = (z[i] - sum) / U[i][i]; } return dB; } // Find leaving variable index int find_leaving_variable( const std::vector<double>& xB, const std::vector<double>& dB, double eps = 1e-10) { double min_ratio = std::numeric_limits<double>::max(); int leaving_idx = -1; for (size_t i = 0; i < xB.size(); ++i) { if (dB[i] > eps) { // Only consider positive denominators double ratio = xB[i] / dB[i]; // Handle degeneration: if xB[i] is near zero, ratio is near zero // But we still need to pick one — use lexicographic rule implicitly via loop order if (ratio < min_ratio - eps) { min_ratio = ratio; leaving_idx = static_cast<int>(i); } else if (std::abs(ratio - min_ratio) < eps) { // Tie-breaking: choose smallest i (Bland's rule to prevent cycling) if (leaving_idx == -1 || static_cast<size_t>(leaving_idx) > i) { leaving_idx = static_cast<int>(i); } } } } if (leaving_idx == -1) { throw std::runtime_error("Unbounded problem: no positive dB[i] found"); } return leaving_idx; }

关键点有三：第一，dB[i] > eps中的eps必须与前面 LU 分解的eps一致（此处为1e-10），否则同一数值在不同模块被判定为“正”和“非正”，逻辑断裂。第二，Bland's rule（选最小下标打破比值相等）被硬编码在else if分支中，这是防止单纯形法在退化问题中无限循环的必要措施。第三，xB[i] / dB[i]计算前不检查xB[i]是否为负——因为xB是基可行解，理论要求xB ≥ 0，若出现负值，说明前面 LU 分解或solve_B_times_x有数值误差，应在preprocess_constraints阶段通过b >= 0断言捕获。

3.3 基变换与解更新：原地更新还是新建矩阵？

pivot 操作的本质是：将换入变量x_e替换换出变量x_l，即基索引数组basis中第l个元素由old_j改为e，同时更新基解xB和非基解xN。高效做法是原地更新xB和A的列，而非重建整个矩阵：

// Update xB: xB = xB - theta * dB for (size_t i = 0; i < xB.size(); ++i) { xB[i] -= theta * dB[i]; } xB[leaving_idx] = theta; // New basic variable takes value theta // Update nonbasic solution: xN[entering_idx] was 0, now becomes theta? No — it's now basic. // So set it to 0 in xN (though it will be removed from xN next iteration) xN[entering_idx] = 0.0; // Update basis vector basis[leaving_idx] = entering_col_global_idx; // entering_col_global_idx 是原始A中的列号 // Update A matrix columns: for each nonbasic column j, do A_j := A_j - (dB[i]/dB[l]) * A_l // But we don't store full A, only need to update future dB calculations — so skip // Instead, we'll recompute dB from scratch next iteration using new basis

注意最后注释：“我们不在本次迭代中更新A矩阵的列，而是在下次迭代中用新basis重新提取列”。这是有意为之的权衡：存储A的完整副本会占用 O(mn) 内存，而单纯形通常迭代次数远小于n，因此每次solve_B_times_x时从原始A提取所需列，比维护一个动态A更省内存且不易出错。这也是为什么preprocess_constraints返回的A_ext必须是完整矩阵——它是唯一真相源。

4. 避坑指南：单纯形 C++ 实现中 4 个血泪经验换来的数值陷阱

写完上述代码，你以为就能跑通教科书例题了？不。我在某跨平台系统中调试这个算法时，花了 3 天定位一个NaN，又花了 2 天修复退化导致的无限循环。以下是真实踩过的坑，按现象→原因→解决排列，每一条都带可复现的测试用例。

4.1 现象：迭代 7 次后xB出现-1.23e-18，后续检验数全为NaN

原因：LU 分解中U[k][k]因浮点累积误差变为1e-17，1/U[k][k]溢出为inf，inf * finite = NaN。根本原因是未在 LU 分解前对A_basis做列归一化（column scaling）。当某列元素量级为1e6，另一列为1e-3时，U矩阵条件数暴增。
解决：在lu_decompose_with_pivot前插入列缩放：

// 在 find_initial_basis 后，构造 A_basis 前 std::vector<double> col_scales(A_basis[0].size(), 1.0); for (size_t j = 0; j < A_basis[0].size(); ++j) { double max_abs = 0.0; for (size_t i = 0; i < A_basis.size(); ++i) { max_abs = std::max(max_abs, std::abs(A_basis[i][j])); } if (max_abs > 1e-12) col_scales[j] = 1.0 / max_abs; } // 然后构造缩放后的 A_basis_scaled[i][j] = A_basis[i][j] * col_scales[j] // LU 分解后，解向量需反向缩放

缩放不改变解的数学性质，但将U对角线拉回O(1)量级，1/U[k][k]安全。

4.2 现象：求解max x1 + x2, s.t. x1 + x2 <= 1, x1, x2 >= 0时，迭代 100+ 次仍未停止

原因：退化（degeneracy）导致连续多次theta = 0，基在顶点处“原地踏步”。虽然xB不变，但basis数组在换入换出，检验数微小波动（1e-15级别）使算法误判未达最优。
解决：强制启用 Bland's rule（已在 3.2 节代码中实现），且在compute_reduced_costs中对检验数做阈值截断：

// 在 compute_reduced_costs 末尾添加 for (double& rc : reduced_costs) { if (std::abs(rc) < 1e-12) rc = 0.0; // 彻底清零微小检验数 }

阈值1e-12比1e-8更激进，因为退化问题中1e-8仍可能触发虚假换入。

4.3 现象：输入含大系数（如c = [1e7, 1e-7]）时，最优值与scipy.optimize.linprog相差 10%

原因：“大M法”中人工变量惩罚系数M = 1e8与目标系数1e7量级相近，导致人工变量未被完全驱逐，残留在最终基中，污染解。
解决：动态计算M，取M = 10 * (sum of absolute values of all original c coefficients)：

double M = 0.0; for (double cj : c_orig) M += std::abs(cj); M *= 10.0; // 确保 M 不溢出 if (M > 1e10) M = 1e10;

并在preprocess_constraints中用此M替换c_ext.push_back(-M)。

4.4 现象：gdb调试时dB向量全为0，min_ratio无法计算

原因：A_enter_col提取错误。entering_idx是非基变量在nonbasic_indices数组中的下标，但A_enter_col应从原始A_ext矩阵中提取第nonbasic_indices[entering_idx]列，而非第entering_idx列。
解决：在solve_B_times_x调用前，明确列索引映射：

int global_enter_col = nonbasic_indices[entering_idx]; // 关键！ std::vector<double> A_enter_col; for (size_t i = 0; i < A_ext.size(); ++i) { A_enter_col.push_back(A_ext[i][global_enter_col]); }

这个错误极隐蔽，因为当nonbasic_indices是[2,3,4]时，entering_idx=0对应global_enter_col=2，若误用0，则提取了第 0 列（可能是基变量列），B⁻¹A_0本应为单位向量，但A_0不在基中，结果必为垃圾值。

提示：所有浮点比较必须用std::abs(a-b) < eps，禁用a == b。eps值在项目中全局定义为const double EPS = 1e-12;，LU 分解、比值判断、检验数清零全部使用同一EPS，避免多标准混乱。

5. 结果解析与对偶信息提取：不只是最优值，还要对偶变量和敏感性分析

单纯形表的最终形态不仅给出最优解x*，还隐含对偶问题的最优解y* = c_BᵀB⁻¹（即检验数计算中的lambda向量），以及资源约束的影子价格（shadow price）。本节展示如何从终止状态中安全提取这些信息，并做基础敏感性分析——比如“如果b[0]增加 1，目标值增加多少”。

5.1 从最终基中提取对偶变量 y* 和影子价格

对偶变量y*就是最终迭代中计算出的lambda向量（见 3.1 节）。但要注意：lambda是c_BᵀB⁻¹，其分量y*[i]表示第i个约束的影子价格，即b[i]每增加 1 单位，目标值的增量（在允许范围内）。提取代码如下：

// After optimal status is reached std::vector<double> y_star = compute_lambda_from_final_basis( L_final, U_final, P_final, c_B_final); // y_star[i] is shadow price for constraint i std::cout << "Shadow prices:\n"; for (size_t i = 0; i < y_star.size(); ++i) { std::cout << "Constraint " << i << ": " << y_star[i] << "\n"; }

compute_lambda_from_final_basis复用 3.1 节的compute_reduced_costs前两步（解Lᵀy = c_B和Uᵀlambda = y），但只返回lambda，不计算检验数。这里的关键是：影子价格仅在b[i]的当前邻域内有效。若b[i]变化过大，基可能改变，影子价格失效。

5.2 资源敏感性分析：b[i] 的允许变化范围

对偶理论指出，当b[i]变化时，只要新b' = b + Δb_i * e_i满足B⁻¹b' ≥ 0，则当前基仍最优，只需更新xB' = B⁻¹b'。

本文还有配套的精品资源，点击获取