线性规划

线性规划（Linear Programming, LP）是运筹学的开山鼻祖，也是数学建模竞赛中处理「在有限资源下做最优决策」的最基础、最实用的工具。凡是「目标与约束都能写成决策变量的线性式子」的优化问题——生产计划、资源分配、运输调运、投资组合、混合配料——几乎都可以用线性规划建模求解。它理论完备（有精确的最优性证明）、求解极快（百万级变量也能秒出结果），还自带「副产品」：影子价格与灵敏度分析，可以直接回答「哪项资源最值钱」「参数能变多少方案不失效」这类评审最爱问的问题。本文从原理、适用场景、指标、可视化到可运行代码（手写单纯形法 + scipy 对照），完整梳理线性规划的竞赛实战用法。

一、算法含义

1.1 通俗理解

某工厂每天用三种资源（原料 A、原料 B、工时）生产两种产品 P1、P2：生产一件 P1 赚 3 万元，生产一件 P2 赚 2 万元；每件 P1 消耗 2 kg 原料 A、1 kg 原料 B、1 小时工时，每件 P2 消耗 1 kg 原料 A、2 kg 原料 B、1 小时工时；每天只有 100 kg 原料 A、150 kg 原料 B、80 小时工时。问：每天各生产多少件，利润最大？

设决策变量 $x_1, x_2$ 为两种产品的日产量，问题可以完整写成：

$\max\; z = 3x_1 + 2x_2$

$\text{s.t.} \quad 2x_1 + x_2 \le 100 \quad (\text{原料A}), \qquad x_1 + 2x_2 \le 150 \quad (\text{原料B})$

$x_1 + x_2 \le 80 \quad (\text{工时}), \qquad x_1 \ge 0,\ x_2 \ge 0 \quad (\text{非负})$

这就是一个线性规划问题：目标函数是决策变量的线性函数（利润），所有约束也是线性等式或不等式，求一组非负决策变量使目标最优。所谓「线性」，是指表达式中每个变量都只出现一次、没有平方/乘积/指数/取对数等运算——$3x_1 + 2x_2$ 是线性的，$3x_1^2 + 2x_1 x_2$ 就不是。别小看这个限制：正因为「线性」，可行域才有漂亮的几何结构（凸多面体），最优解才会落在顶点上，算法才能又快又保证最优。

1.2 一般形式与标准型

线性规划的一般形式（本文采用的记号，与第六节代码一致）：

$\min_{x}\; c^\top x \ \text{s.t. } Ax \le b,\ x \ge 0$

其中 $x \in \mathbb{R}^n$ 是决策变量向量，$c \in \mathbb{R}^n$ 是目标系数向量，$A \in \mathbb{R}^{m \times n}$ 是约束系数矩阵，$b \in \mathbb{R}^m$ 是右端项（资源量）。若题目要求最大化（如利润），利用恒等式

$\max\; c^\top x \iff \min\; -c^\top x$

即可转化为最小化形式（最优解相同，最优值相差一个负号）。任何线性规划都能化成标准型（等式约束 + 非负变量 + 非负右端项）：

$\min\; c^\top x \quad \text{s.t.} \quad Ax = b,\ x \ge 0,\ b \ge 0$

常见转化技巧如下：

原形式	转化方法
$\max\; c^\top x$	等价于 $\min\; -c^\top x$，最后目标值取负号还原
$a_i^\top x \ge b_i$	两边乘 $-1$：$-a_i^\top x \le -b_i$
$a_i^\top x = b_i$	拆成两个不等式：$a_i^\top x \le b_i$ 且 $-a_i^\top x \le -b_i$
$x_j$ 无符号限制（自由变量）	令 $x_j = x_j^+ - x_j^-$，$x_j^+, x_j^- \ge 0$
$a_i^\top x \le b_i$	引入松弛变量 $s_i \ge 0$：$a_i^\top x + s_i = b_i$

松弛变量的物理意义是「资源剩余量」：约束 $2x_1 + x_2 \le 100$ 写成 $2x_1 + x_2 + s_1 = 100$ 后，$s_1 = 100 - (2x_1 + x_2)$ 就是没用完的原料 A。它在单纯形法与互补松弛理论中都扮演关键角色。

1.3 几何意义：凸多面体与顶点最优

可行域是所有约束「同时满足」的点集：

$S = \{ x \mid Ax \le b,\ x \ge 0 \}$

每个不等式 $a_i^\top x \le b_i$ 是一个半空间（二维中是直线的一侧，三维中是平面的某一侧），可行域就是 $m + n$ 个半空间的交集，是一个凸多面体。凸性的含义是：可行域内任意两点的连线整段仍在可行域内（$\lambda x + (1-\lambda) y \in S$），这保证「局部最优 = 全局最优」，也是单纯形法不会陷进「局部坑」的根本原因。

目标函数 $z = c^\top x$ 的**等值线（面）是一族互相平行的直线（超平面）：$c^\top x = k$，法向量就是 $c$。让 $k$ 沿梯度方向 $c$ 平移，等值线扫过可行域时，最后离开的那个点就是最优解。由于可行域是「有棱有角」的多面体，最优解必然在多面体的某个顶点（极点）**上（若可行域有界且最优解存在）。这就是线性规划的核心定理：

顶点最优定理：若线性规划有有限最优值，则至少在一个顶点处达到最优。顶点 = 若干条约束边界（$Ax = b$ 中的独立面 + 坐标面 $x_j = 0$）的交点，即「基本可行解」（basis feasible solution）。

直观理解：把可行域想象成一块凸多边形的板子，目标等值线像一把尺子沿 $c$ 方向平移，尺子最后离开板子的位置一定是一个「角」而不是边的中点（除非等值线与某条边平行——此时整条边都是最优解，称为多解/退化情形）。顶点的个数有限（最多 $C_{m+n}^{m}$ 个候选），所以「在所有顶点中挑最好的」是一个有限搜索问题——这正是单纯形法的思路。

1.4 单纯形法（Simplex Method）思想与迭代骨架

单纯形法是 Dantzig 于 1947 年提出的顶点间「爬山」算法：从一个顶点出发，沿可行域的一条边走到目标值更优的相邻顶点，重复直到无法改进为止。每一轮迭代分四步：

1. 选基（得到当前顶点）：把变量分成基变量 $x_B$（$m$ 个）与非基变量 $x_N$（$n$ 个，取 0）。由等式约束 $Ax = b$ 解出：

$x_B = B^{-1} b \ge 0, \qquad x_N = 0$

其中 $B$ 是由基变量对应列组成的 $m \times m$ 可逆矩阵（基矩阵）。只要 $B^{-1}b \ge 0$，$(x_B, x_N)$ 就是一个顶点（基本可行解）。

2. 算检验数（判断是否最优）：对最小化问题，非基变量的**检验数（简约成本）**为

$r_N^\top = c_N^\top - c_B^\top B^{-1} N$

若所有 $r_j \ge 0$，则当前顶点最优，停止；否则取 $r_j < 0$ 的列 $j$ 作为入基变量（目标值沿该方向下降；Bland 规则选最小下标可防止循环）。

3. 比值检验（决定出基变量）：入基变量从 0 增大，最先被「顶到 0」的基变量出基。最大可增量为

$\theta = \min_{i} \left\{ \frac{(B^{-1}b)_i}{(B^{-1}a_j)_i} :\ (B^{-1}a_j)_i > 0 \right\}$

取到最小值的行 $r$ 对应的基变量出基。若所有 $(B^{-1}a_j)_i \le 0$，比值检验全为 $+\infty$，说明目标值可以无限改善——问题无界。

4. 换基（移动到相邻顶点）：交换入基列与出基列，更新 $B$，回到第 2 步。实际计算中整个过程用单纯形表（高斯-若尔当消元）完成：枢轴行除以主元，其余行消去主元列，目标行同时更新。目标函数值恰好出现在表的右下角。

初始化（第 1 步的初始顶点）一般靠两阶段法或大 M 法：当约束是 $Ax \le b,\ b \ge 0$ 时，松弛变量 $s = b$ 天然给出初始顶点（本文实例即如此）；一般情况先解阶段一辅助问题 $\min \sum_i a_i$（$a_i$ 为人工变量），其最优值为 0 说明原问题可行，并得到一个可行基；阶段二再以它出发求原问题最优。若阶段一最优值 $> 0$，说明不存在满足全部约束的点——原问题不可行。第六节代码实现了完整的两阶段单纯形法。

1.5 对偶理论与影子价格

每个线性规划（原问题）都伴随一个「镜像」问题（对偶问题）。对本文的 max 形式：

$\text{原问题:} \quad \max\; c^\top x \quad \text{s.t.} \quad Ax \le b,\ x \ge 0$

$\text{对偶问题:} \quad \min\; b^\top y \quad \text{s.t.} \quad A^\top y \ge c,\ y \ge 0$

对偶变量 $y \in \mathbb{R}^m$ 与 $m$ 个资源约束一一对应。以 1.1 的实例为例，对偶问题为：

$\min\; w = 100y_1 + 150y_2 + 80y_3 \quad \text{s.t.} \quad 2y_1 + y_2 + y_3 \ge 3,\quad y_1 + 2y_2 + y_3 \ge 2,\quad y \ge 0$

三条核心性质：

• 弱对偶：任意可行解都有 $c^\top x \le b^\top y$，即对偶问题的最优值永远是原问题最优值的「上界」；
• 强对偶：若原问题有最优解，则对偶问题也有最优解且二者相等：$c^\top x^* = b^\top y^*$。对偶间隙 $gap = w^* - z^*$ 等于 0 就是「已达全局最优」的数学证明；
• 互补松弛：最优时 $(b_i - a_i^\top x^*) \cdot y_i^* = 0$——某项资源有剩余（$b_i - a_i^\top x^* > 0$）则其影子价格必为 0；影子价格 $> 0$ 的资源必然恰好用完。

影子价格 $y_i$ 的经济含义是第 $i$ 项资源的边际价值：在最优基不变的前提下，资源 $i$ 每增加 1 个单位，最优目标值近似增加 $y_i$ 个单位（严格地说是最优值曲线在该点的斜率）：

$y_i = \frac{\partial z^*}{\partial b_i} = (c_B^\top B^{-1})_i$

本文实例中 $y = (1, 0, 1)$：每多 1 kg 原料 A 利润多 1 万元，工时每多 1 小时利润多 1 万元，而原料 B 每天剩 10 kg 用不完，所以 $y_2 = 0$——「值钱的资源值得买，用不完的资源买来也白搭」，这就是影子价格对决策的指导意义。

1.6 灵敏度分析

竞赛论文中评审必问：「你的方案对数据波动稳健吗？系数变多少方案会失效？」灵敏度分析回答的正是这个问题。它的结论是在「当前最优基不变」的前提下，各系数的允许变化区间：

• 右端项 $b_i$ 的区间：$b$ 变化时最优解 $x_B = B^{-1}b$ 随之线性变化，只要 $B^{-1}(b + \Delta b\, e_i) \ge 0$（基变量仍非负），最优基就不变。于是 $\Delta b_i$ 的允许范围由不等式组解出，记为 $b_i \in [\underline{b_i}, \overline{b_i}]$。区间内影子价格不变、$z^*(b)$ 是斜率为 $y_i$ 的线性函数；区间外最优基改变（图 3 中的转折点）。
• 目标系数 $c_j$ 的区间：$c$ 变化不改变可行域，只改变检验数。让所有非基变量检验数保持非负（min 形式），即得 $c_j \in [\underline{c_j}, \overline{c_j}]$。区间内最优解 $x^*$ 不变（最优值按 $x_j^*$ 的比例线性变化）；区间外需重新求解。

灵敏度区间只保证「解结构不变」，是局部结论；超出区间的外推（如把影子价格线性外推到资源翻倍）是常见的误用，务必注意（见 7.2）。

1.7 优缺点

优点：

• 模型结构简单直观，「资源—消耗—目标」三张表即可建模，评审一眼看懂；
• 理论完备：强对偶给出精确最优性证明（$gap = 0$），不会像启发式算法那样「不知道解好不好」；
• 求解极快：单纯形法在实践中平均迭代 $O(m)$ 次，HiGHS 等现代求解器可解百万变量的实际问题，竞赛数据规模完全无压力；
• 自带决策信息：影子价格（哪些资源值得扩充）、灵敏度区间（方案稳健性）、简约成本（该不该生产某产品）都是论文的加分项；
• 是整数规划、网络流、多目标规划等高级方法的基础（08 篇的整数规划就在 LP 上加「取整」约束）。

缺点：

• 要求目标与约束均为线性：现实中的规模效应（边际成本递减）、价格弹性、$x_1 \cdot x_2$ 型耦合都是非线性的，需改用非线性规划或分段线性化近似；
• 变量必须连续可分割：产品件数、员工人数、是否投资这类整数决策会得到 $x_1 = 3.7$ 这类不可执行解，需改用整数规划（直接四舍五入既不保证可行也不保证最优）；
• 是静态单期模型：不考虑时间先后与随机因素，多阶段决策需动态规划、随机数据需随机规划/鲁棒优化；
• 最优解可能出现多解/退化（等值线平行于某条约束边），报告最优方案时需要说明不唯一；
• 灵敏度结论是局部的，区间外必须重新求解，不能外推。

二、何时使用（适用场景与条件）

2.1 适用场景

1. 资源分配：若干稀缺资源（资金、人力、原料、设备）在若干竞争性活动之间的最优分配。这是 LP 最经典的形态，本文生产计划实例即属此类。
2. 生产计划：多产品、多工序、多资源限制下的产量组合与排产（可加库存约束推广为多周期）。
3. 运输问题：$m$ 个产地供应 $n$ 个销地，求运费最小的调运方案（LP 的特例，有专门的表上作业法，见第八节）。
4. 投资组合：在风险、行业分散、流动性等线性约束下最大化期望收益（Markowitz 均值-方差模型本身是二次的，但「线性收益 + 线性约束」的版本是 LP）。
5. 混合配料：饲料配比、汽油调和、合金熔炼——在满足营养成分/质量标准（线性不等式）下最小化成本。
6. 下料问题：把标准长度的原材料切割成订单要求的小段，使余料最少（变量是各切割模式的份数；若要求整数份则升级为整数规划）。
7. 人员排班：各班次人数需求（线性不等式）下安排员工数量使总人力成本最小。
8. 网络流/最短路/最大流：都是 LP 的特殊结构，可用 LP 求解器统一处理（数据包络分析 DEA、公平分配问题同理）。

2.2 数学建模竞赛中的典型题目

• 调度与分配类：机场航班机位分配、车队调度、志愿者/救援物资分配——本质都是「资源受限下的最优指派」，LP 是最直接的建模起点；
• 生产与库存类：多周期生产计划（本周期产量 + 上周期库存 − 本期需求 = 本期库存）、工厂选址配送一体化（若加 0/1 选址变量则转整数规划）；
• 经济与能源类：电力系统机组出力分配（经济调度）、碳排放配额分配；
• 数据驱动类：DEA 效率评价、把某些回归问题松弛为 LP（最小一乘回归 $L_1$ 回归可化为 LP，可对比 01~06 篇的最小二乘）。

竞赛实战提醒：很多看起来「不像优化」的题目，一旦写出「目标 + 线性约束」的数学形式就是 LP。建模阶段把 LP 写清楚，即使最终需要整数/非线性推广，LP 也是基线与对照——先解松弛版拿到理论上界，再讨论差距，是评审最喜欢的叙事结构。

2.3 使用前提（建模前检查清单）

检查项	具体要求	检查手段
目标线性	目标是决策变量的线性组合 $c^\top x$（可求和、可加权平均）	写出目标式，确认无 $x_i x_j$、$x_i^2$、$\ln x_i$ 等
约束线性	每个约束都是线性等式/不等式（资源消耗量 = 系数 × 产量）	逐条约束写成 $a^\top x \le b$ 形式
变量连续	变量可任意分割（产量、资金、时间；件数/人数不行）	决策变量的物理含义
确定性数据	系数 $c, A, b$ 已知且确定（或至少能给出名义值/均值）	数据来源是否可靠
单目标	只有一个目标函数（多目标需先加权/ε-约束转化，见 2.5）	题目的「指标」是否唯一
约束不矛盾	可行域非空（多个约束可能互斥）	求解后检查 status（不可行则回到建模检查约束方向）

竞赛提示：以上检查不必在论文中逐条罗列，但至少把目标函数与全部约束以数学公式完整写出，并注明每个系数的单位与来源——「公式清晰」是优化类论文得分的第一要素。

2.4 不适用情形（考虑替代方法）

1. 非线性目标/约束：边际成本递增、规模经济、价格随产量变化（$p(x) \cdot x$ 型收益）、化学反应速率等。改用非线性规划（04 篇思路），或对光滑非线性做分段线性近似后回到 LP（需在论文中说明近似误差）。
2. 整数解要求：产品件数、车辆台数、项目「投或不投」、选址「建或不建」。此时 LP 的松弛解可能不可执行，应使用整数规划/0-1 规划（08 篇，scipy.optimize.milp 可直接求解）。
3. 多目标权衡：既要利润最大又要污染最小又要就业最多，单目标 LP 无法表达。先用加权求和、ε-约束法或目标规划转成单目标，再求解。
4. 随机/不确定数据：需求、价格、资源量是随机变量时，确定性 LP 的最优解可能对波动极脆弱。升级为随机规划、鲁棒优化（LP 的鲁棒对应是锥规划），或至少配合灵敏度分析与情景分析。
5. 决策相互嵌套的博弈/动态过程：多方博弈用博弈论（可化为 LP 的零和博弈除外），多阶段决策用动态规划，连续动态过程用最优控制。

2.5 与整数规划、非线性规划、多目标规划的选择

方法	适用情形	与线性规划的关系
线性规划 LP	线性目标 + 线性约束 + 连续变量	本文主角，一切优化方法的基石
整数规划 ILP/MILP	部分变量必须取整数（台数、0/1 决策）	LP + 取整约束；先解 LP 松弛得上界，再用分支定界（08 篇）
非线性规划 NLP	目标或约束含非线性函数	LP 的推广；无全局最优保证，需初值与凸性判断
多目标规划	多个互相冲突的目标	加权法把目标合并成 $\sum w_k f_k$ 后仍可用 LP 求解器
目标规划 GP	目标有优先级、允许偏差（软约束）	LP 变体：变量换成正负偏差量 $d^+, d^-$，约束加 $\sum d$ 项
网络流/最短路/最大流	网络结构问题	LP 的特例，也可用专门算法（效率更高）
随机规划/鲁棒优化	数据不确定	LP 的对偶理论在鲁棒/随机对偶中仍起作用

选择原则（竞赛实战）：先尝试写成 LP → 检查变量是否连续、目标约束是否线性 → 连续且线性，直接 LP；要取整，升级整数规划；非线性，先试分段线性化（LP 可解）再考虑 NLP；多目标，先加权/ε-约束。论文按「LP 基线 → 根据实际约束升级」的递进逻辑写，每一步升级都要给出「为什么必须升级」的证据。

三、算法指标

线性规划的输出不只是「最优解」，而是一整套决策信息。下面按第六节代码的输出顺序，给出每个指标的中文名、公式、含义与解读；符号约定见第五节。

3.1 最优目标函数值 $z^*$ 与最优解 $x^*$

$z^* = c^\top x^* = \max_x \left\{ c^\top x \mid Ax \le b,\ x \ge 0 \right\}$

含义：最优方案 $x^*$ 及其对应的目标值。解读：$x^*$ 各分量就是「每种产品生产多少」的行动方案，$z^*$ 是该方案下的最大利润/最小成本。报告时务必带单位（本文实例 $x^* = (20, 60)$ 件/日、$z^* = 180$ 万元/日），并给出最优基（哪些约束起「卡脖子」作用）。注意 max 问题经 $c \to -c$ 转换求解后，求解器返回的函数值要取负号还原（scipy linprog 即如此）。

3.2 对偶间隙 gap

$gap = w^* - z^* = b^\top y^* - c^\top x^* \ge 0$

含义：对偶问题最优值 $w^*$ 与原问题最优值 $z^*$ 之差。由弱对偶定理 gap 恒非负；强对偶定理保证原问题有最优解时 $gap = 0$。解读：$gap = 0$（数值上小于 $10^{-9}$ 量级即可）是「当前解确为全局最优」的严格数学证明，比任何启发式算法的「收敛了」都硬气——竞赛论文中一句「对偶间隙为 0，解为全局最优」直接封死质疑。若 gap 显著大于 0 且求解器未报错，通常是数值问题或模型未收敛。

3.3 影子价格（对偶变量）$y$

$y_i = \frac{\partial z^*}{\partial b_i} = (c_B^\top B^{-1})_i$

含义：第 $i$ 项资源的边际价值——资源每增加 1 单位，最优目标值近似增加 $y_i$。解读：$y_i > 0$ 说明该资源是瓶颈（恰好用完、约束是紧的），值得购买扩充；$y_i = 0$ 说明该资源有剩余（约束松），白送也不用。由互补松弛，$y_i$ 与松弛量 $s_i$ 至少一个为 0。注意影子价格只在灵敏度区间（3.4）内有效。求解途径有三种，可互相印证：① 从最优单纯形表目标行中松弛变量列读出；② 求解器返回的边际值（HiGHS 的 ineqlin.marginals，注意 min 形式要取负号）；③ 直接解对偶问题。

3.4 灵敏度区间（允许变化范围）

右端项 $b_i$ 的允许范围（保持当前最优基不变）：

$\Delta b_i \in [\Delta^-,\ \Delta^+] \quad \text{由} \quad B^{-1}(b + \Delta b_i\, e_i) \ge 0 \quad \text{解出，区间为} \ b_i \in [\underline{b_i}, \overline{b_i}]$

目标系数 $c_j$ 的允许范围（最优基不变 ⟺ 所有非基检验数保持非负）：

$c_j \in [\underline{c_j}, \overline{c_j}] \quad \text{由} \quad r_N(c_j) \ge 0 \ (\text{min 形式}) \quad \text{解出}$

含义与解读：区间内「哪些约束卡脖子、哪个顶点最优」的结构不变——$b_i$ 区间内影子价格 $y_i$ 不变（最优值随 $b_i$ 线性变化），$c_j$ 区间内最优解 $x^*$ 完全不变。区间外最优基改变，需重新求解。竞赛论文中的标准话术：「当原料 A 拥有量在 $[90, 160]$ kg 内波动时，当前最优生产方案的结构不变」——这是回答稳健性的核心证据。

3.5 简约成本（reduced cost）$r_j$

max 形式下（与对偶变量挂钩）：

$r_j = c_j - y^\top a_j$

最优时所有 $r_j \le 0$：基变量 $r_j = 0$（正在生产的产品，利润已被充分挖掘）；非基变量 $r_j < 0$，$|r_j|$ 的经济含义是「强行生产 1 单位该产品会造成的机会成本损失」。松弛变量的简约成本恰为 $-y_i$（用不完的资源，其「成本」就是影子价格）。min 形式下符号相反（$r_j = c_j - c_B^\top B^{-1}a_j \ge 0$）。解读示例：5 变量实例中 $r_{x_2} = -1.1667$，意味着「每强行生产 1 件 x2，总利润将损失 1.1667 万元」——这回答「为什么不生产 x2」的质疑。

3.6 状态判定：不可行 / 无界

• optimal（最优）：可行域非空且有界的目标方向，返回最优解；
• infeasible（不可行）：约束互相矛盾、可行域为空（两阶段法中阶段一最优值 $> 0$，或求解器 status=2）。典型原因：约束方向写反、限制过严（如「产量 ≤ 100 且 ≥ 200」）；
• unbounded（无界）：目标值可无限改善（如最大化利润但漏写了资源约束；两阶段法中比值检验全为 $+\infty$，或求解器 status=3）。注意：无界是「建模错误」的信号而不是「利润无限」的好消息——一定是漏约束了；
• 退化（degenerate）：某个基变量取值为 0（多于 $n$ 个约束边界交于一点），单纯形法可能出现循环（用 Bland 规则规避），最优解通常不唯一。退化不是错误，但要在论文中说明「最优解不唯一」。

3.7 求解时间

手写单纯形法（纯 Python、教学实现）与工业级求解器（HiGHS）的耗时可分别计时报告。竞赛数据规模下差异不大（毫秒级），但作为方法论说明很有价值：「本模型规模为 2 变量 3 约束，HiGHS 求解耗时约 1 ms；模型可扩展至更大规模」。注意不同机器计时会有波动，论文中报数量级即可。

3.8 指标汇总表

指标	公式	含义与解读	本文实例取值
最优解 $x^*$ / 最优值 $z^*$	$z^* = c^\top x^*$	最优行动方案与目标值	$x^* = (20, 60)$ 件/日，$z^* = 180$ 万元/日
对偶间隙 gap	$gap = w^* - z^*$	$=0$ 即全局最优的证明	$w^* = 180$，gap $= 0$
影子价格 $y$	$y_i = \partial z^* / \partial b_i = (c_B^\top B^{-1})_i$	资源边际价值；$>0$ 稀缺，$=0$ 有剩余	$y = (1, 0, 1)$ 万元/单位
$b$ 的灵敏度区间	$B^{-1}(b + \Delta b_i e_i) \ge 0$	最优基不变时资源的波动范围	$b_1 \in [90, 160]$，$b_2 \in [140, +\infty)$，$b_3 \in [50, 83.33]$
$c$ 的灵敏度区间	检验数不变号	最优解不变时利润系数的波动范围	$c_1 \in [2, 4]$，$c_2 \in [1.5, 3]$
简约成本 $r_j$	$r_j = c_j - y^\top a_j$（max 形式）	基变量为 0；非基 $\|r_j\|$ = 强行生产的机会成本	$r_x = (0, 0)$，$r_s = (-1, 0, -1)$
状态判定	status ∈ {optimal, infeasible, unbounded}	诊断建模错误	optimal（无退化）
求解时间	计时统计	报告数量级即可（随机器波动）	手写 ≈ 0.1 ms，HiGHS ≈ 1.0 ms

四、可视化图表

线性规划的论文图承担三件事：展示可行域与最优点的几何关系、展示算法如何找到最优点、展示最优值随参数变化的规律。以下 4 张图是竞赛 LP 的标准配置，均由第六节代码生成（保存于 figures/ 目录）。

4.1 四张标准图汇总

图名	用途	关键解读点
① lp_feasible_region.png（二维可行域 + 目标等值线 + 最优点）	图解法核心：一张图讲完「可行域—目标—最优解」	浅蓝填充的多边形是可行域（5 个顶点，凸）；三条彩色直线是资源约束边界（原料A 橙色、原料B 绿色、工时 紫色），箭头方向为可行侧；灰色虚线族是目标等值线 $3x_1 + 2x_2 = k$（$k = 60, 120, 180$），沿梯度方向 $c = (3,2)$ 平移目标值增大；红色五角星标出最优解 $(20, 60)$，位于「原料A 与工时」两条约束边的交点——两条资源恰好用完，几何上解释了为什么 $y_1, y_3 > 0$
② lp_simplex_path.png（单纯形迭代顶点路径）	展示算法的顶点移动过程	同一可行域上画出橙色折线：从初始顶点 $(10, 70)$（Phase I 求得，$z = 170$）沿边移动到最优顶点 $(20, 60)$（$z = 180$）；箭头注释标出每步换基前后的目标值变化；路径只走「顶点—边—顶点」，且每步目标值严格上升，直观演示单纯形法的爬山机制与有限步终止
③ lp_b_curve.png（右端项 $b_1$ 变化时最优值曲线）	展示影子价格与灵敏度区间	横轴 $b_1$（原料A 拥有量）从 70 到 175 kg，纵轴最优利润 $z^*(b_1)$；曲线是分段线性的，共 4 段、3 个转折点（$b_1 = 75, 90, 160$）；每段斜率分别是 $2, 4/3, 1, 0$，各段斜率 = 该区段原料A 的影子价格；灰色竖虚线标出当前点 $b_1 = 100$（斜率 = 1 = $y_1$）；最右段斜率为 0 说明原料A 多到不再稀缺（瓶颈转移到工时）；转折点 90 与 160 正是当前最优基灵敏度区间 $[90, 160]$ 的端点——区间内线性、区间外拐弯
④ lp_shadow_price_bar.png（各资源影子价格柱状图）	回答「哪项资源最值钱」	三根柱子对应原料A、原料B、工时，高度为影子价格 $y = (1, 0, 1)$；正影子价格柱用蓝色、0 值柱用灰色（并标注「剩余 10 kg，增加投入不改变最优利润」），一眼区分瓶颈资源与冗余资源；柱顶标数值，纵轴单位「万元/单位资源」；竞赛答辩时指这张图说「建议优先扩充原料A 与工时」

4.2 好图与异常图的特征

好图的特征

• 图①：可行域填充色柔和、约束线颜色区分且带约束名称标注、等值线方向与梯度箭头一致、最优点用醒目标记并与等值线最高值位置吻合；顶点坐标全部标出，方便评审核验「最优解在顶点」；
• 图②：路径沿可行域边界（边）走、每步标注目标值且单调上升、起终点名称清晰（初始顶点 → 最优顶点）；
• 图③：曲线分段线性、转折点（kinks）用圆点标出并注明横坐标、每段斜率有标注、当前工作点用竖虚线标出；转折点数量与灵敏度区间边界一致；
• 图④：正影子价格与零影子价格视觉区分（颜色或灰度）、每根柱标数值、对零值给出经济解释（哪项资源有剩余）。

异常图的特征（出现即需要排查）

• 图①可行域「缺角」或最优解画在可行域外 → 约束方程写错或绘图代码中的边界函数与模型不一致；
• 图①等值线方向与梯度箭头相反（箭头指向目标值减小的方向）→ 梯度方向写反，或 max/min 混淆；
• 图②路径「穿膛而过」（不沿边）或目标值下降 → 单纯形实现有误（比值检验或入基规则错），需回到迭代公式逐行核对；
• 图③曲线不线性（弯曲）→ 分段区间取点过稀或求解失败（status 非 0）；转折点漏标 → 灵敏度区间未与图核对；
• 图④出现负的柱子 → 影子价格符号取反（min/max 转换时边际值忘了取负号），负影子价格在标准 LP 中不可能出现。

五、符号说明

符号	含义	示例/单位
$n$	决策变量个数	本文小实例 $n = 2$（两种产品）
$m$	资源约束个数	本文小实例 $m = 3$（三种资源）
$x$	决策变量向量，$n \times 1$	$x = (x_1, x_2)$ 件/日
$x_j$	第 $j$ 个决策变量	$x_1$ = 产品 P1 日产量（件）
$c$	目标系数向量，$n \times 1$	$c = (3, 2)$ 万元/件
$A$	约束系数矩阵，$m \times n$	$a_{ij}$ = 第 $i$ 种资源对第 $j$ 种产品的单耗
$a_j$	$A$ 的第 $j$ 列（产品 $j$ 的消耗向量）	$a_1 = (2, 1, 1)^\top$
$b$	右端项（资源拥有量），$m \times 1$	$b = (100, 150, 80)$
$z$	目标函数值	$z = c^\top x$，万元
$z^*$、$x^*$	最优目标值、最优解	$z^* = 180$，$x^* = (20, 60)$
$s$	松弛变量向量（资源剩余量），$m \times 1$	$s_2 = 10$ kg（原料B 剩余）
$S$	可行域（凸多面体）	$S = \{x \mid Ax \le b, x \ge 0\}$
$B$	基矩阵（$m \times m$ 可逆，基变量对应列）	由 $x_1, x_2$ 与 $s_2$ 的列组成
$N$	非基矩阵（其余 $n$ 列）	与 $B$ 互补
$x_B$、$x_N$	基变量、非基变量	$x_B = B^{-1}b$，$x_N = 0$
$r_j$	简约成本/检验数	max 形式 $r_j = c_j - y^\top a_j$
$\theta$	比值检验的最小比值（决定出基变量）	$\theta = \min_i\{(B^{-1}b)_i/(B^{-1}a_j)_i\}$
$y$	对偶变量（影子价格）向量，$m \times 1$	$y = (1, 0, 1)$ 万元/单位资源
$w$	对偶问题目标值	$w = b^\top y$，最优时 $w^* = z^*$
gap	对偶间隙	$gap = w^* - z^* \ge 0$，最优时 $= 0$
$[\underline{b_i}, \overline{b_i}]$	右端项 $b_i$ 的灵敏度区间	$b_1 \in [90, 160]$ kg
$[\underline{c_j}, \overline{c_j}]$	目标系数 $c_j$ 的灵敏度区间	$c_1 \in [2, 4]$ 万元/件
$\mathrm{status}$	求解状态：optimal / infeasible / unbounded	scipy 返回码 0 / 2 / 3
$\min / \max$	目标方向	$\max c^\top x \iff \min -c^\top x$
$e_i$	第 $i$ 个单位向量	灵敏度分析中表示「只动 $b_i$」

六、可运行程序（完整代码）

运行环境：Python 3.12，仅依赖 numpy / scipy / matplotlib（本机若未安装可用 pip install numpy scipy matplotlib 安装）。以下所有代码块按顺序拼接保存为一个 .py 文件即可直接运行（如 lp_demo.py），无需任何外部数据文件；控制台打印第三节全部指标，并在 figures/ 子目录生成第四节 4 张图（文件名与 4.1 节表格一致）。

数据设计说明：主实例为 2.1 节的生产计划问题——max z = 3x1 + 2x2，三种资源约束（原料A、原料B、工时）。数字特意选得「好画」：可行域只有 5 个顶点，最优解 $(20, 60)$ 是原料A 直线与工时直线的整点交点，手算即可核对（$z^* = 180$）。另构造一个 5 变量 4 约束实例演示多变量求解（最优解 $x^* = (10, 0, 0, 10, 0)$，$z^* = 80$）。手写部分实现两阶段单纯形法（Bland 规则防循环），并用 scipy.optimize.linprog（HiGHS 求解器）对照；影子价格用三条途径交叉核对；灵敏度区间由最优表反推；最后绘制全部 4 张图。

# -*- coding: utf-8 -*-
# ============================================================
# 线性规划（Linear Programming）完整示例脚本
# 环境要求：Python 3.12 + numpy / scipy / matplotlib
# 运行方式：python lp_demo.py（在本文档目录运行，自动创建 figures/ 保存图片）
# 数据：脚本内构造的生产计划实例，无需任何外部文件
# 内容：手写两阶段单纯形法 → scipy.optimize.linprog(HiGHS) 对照
#       → 对偶问题与影子价格三途径核对 → 灵敏度分析 → 顶点枚举
#       → 不可行/无界判定演示 → 5 变量实例 → 四张标准图
# ============================================================
import os
import time
import warnings
from itertools import combinations

import numpy as np
import matplotlib

# ========== 0. 基础设置：matplotlib 后端 + 中文字体 ==========
# 先尝试交互式后端（本地运行 plt.show() 会弹出图像窗口）；
# 若当前环境无图形界面（服务器 / CI / 沙箱），自动回退到 Agg 后端，
# 此时 plt.show() 为空操作，脚本照常保存图片、不会报错。
try:
    import matplotlib.pyplot as _plt_probe
    _plt_probe.figure()          # 探测：能否真正创建图形窗口
    _plt_probe.close("all")
except Exception:
    matplotlib.use("Agg", force=True)

import matplotlib.pyplot as plt
from scipy.optimize import linprog

# matplotlib 中文显示设置（macOS 用 PingFang SC，Windows 自动回退 SimHei）
plt.rcParams["font.sans-serif"] = ["PingFang SC", "Arial Unicode MS", "SimHei"]
plt.rcParams["axes.unicode_minus"] = False   # 正常显示负号

FIG_DIR = "figures"                          # 图片输出目录
os.makedirs(FIG_DIR, exist_ok=True)          # 不存在则自动创建

# ========== 1. 问题建模：某厂生产计划（贴近竞赛的小实例） ==========
# 决策变量: x1 = 产品 P1 日产量(件), x2 = 产品 P2 日产量(件)
# 目标(max 利润, 万元/日):   max z = 3*x1 + 2*x2
# 约束:
#   原料A:  2*x1 + 1*x2 <= 100  (kg/日)
#   原料B:  1*x1 + 2*x2 <= 150  (kg/日)
#   工时:   1*x1 + 1*x2 <=  80  (h/日)
#   非负:   x1 >= 0, x2 >= 0
# 手算参考: 联立 2*x1+x2=100 与 x1+x2=80 得交点 (20, 60), 代入目标 z = 180;
#           程序将用单纯形法独立求解并核对。

C = np.array([3.0, 2.0])                  # 目标系数(单位利润, 万元/件)
A = np.array([[2.0, 1.0],                 # 原料A 消耗系数
              [1.0, 2.0],                 # 原料B 消耗系数
              [1.0, 1.0]])                # 工时消耗系数
B = np.array([100.0, 150.0, 80.0])        # 资源拥有量(右端项)
RES_NAMES = ["原料A", "原料B", "工时"]

print("=" * 62)
print("【问题】某厂生产计划(2 产品 x 3 资源)")
print("  max z = 3*x1 + 2*x2")
print("  s.t. 2*x1 + 1*x2 <= 100 (原料A),  1*x1 + 2*x2 <= 150 (原料B),")
print("       1*x1 + 1*x2 <=  80 (工时),  x1, x2 >= 0")
print("  手算参考: x* = (20, 60), z* = 180")

# ========== 2. 手写两阶段单纯形法（完整实现） ==========
def _pivot(T, basis, row, col):
    """高斯-若尔当主元消去: 让第 row 行第 col 列变为 1, 其余行该列消为 0。"""
    T[row] = T[row] / T[row, col]
    for i in range(T.shape[0]):
        if i != row and abs(T[i, col]) > 1e-15:
            T[i] = T[i] - T[i, col] * T[row]
    basis[row] = col


def two_phase_simplex(c, A, b, maximize=True, tol=1e-9):
    """
    两阶段单纯形法: 求解  max/min c^T x  s.t.  A x <= b, x >= 0 (b 可为任意符号)。
    【第一步】化成标准型: 每个 "<=" 约束引入松弛变量 s, 得  A x + s = b, x, s >= 0;
              若某行 b_i < 0, 先整行乘 -1 (此时该行松弛变量系数为 -1, 相当于"剩余变量");
              max 问题等价于 min d^T x (d = -c)。
    【第二步】Phase I 找初始基可行解: 每行再引入人工变量 a, 解辅助问题 min 1^T a;
              若辅助问题最优值 > 0, 说明原问题无可行解(不可行);
              本例约束全是 "<=" 且 b >= 0, 松弛变量本身即可构成可行基。
    【第三步】Phase II 求最优: 删除人工变量列, 恢复原目标 d, 继续换基迭代直到最优。
    换基规则(Bland 规则防循环):
              入基: 选"下标最小"的负简约成本列 (min 问题中 r_j < 0 才可能改进目标);
              出基: 最小比值检验 theta = min{ x_Bi / (B^{-1} a_j)_i }, 只对正系数行取比,
                    平局时选基变量下标最小的行。
    返回 dict: x(最优解), z(最优值), status, iterations, iters_p1, iters_p2,
              basis(最优基变量下标), path(Phase II 顶点轨迹, 画图用), tableau(最终单纯形表)
    """
    c = np.asarray(c, dtype=float)
    A = np.asarray(A, dtype=float).copy()
    b = np.asarray(b, dtype=float).copy()
    d = -c if maximize else c.copy()       # 标准型目标: min d^T x
    m, n = A.shape
    # 预处理: b_i < 0 的行整行乘 -1, 该行松弛变量系数记为 -1 (剩余变量)
    sign = np.ones(m)
    for i in range(m):
        if b[i] < 0:
            A[i] = -A[i]
            b[i] = -b[i]
            sign[i] = -1.0
    n_tot = n + 2 * m                      # 变量列: x(1..n) | s(n+1..n+m) | a(n+m+1..n+2m)
    # 初始表: 约束行 [A | diag(sign) | I | b], 目标行(Phase I): min 1^T a
    T = np.zeros((m + 1, n_tot + 1))
    T[:m, :n] = A
    T[:m, n:n + m] = np.diag(sign)         # 松弛/剩余变量列
    T[:m, n + m:n_tot] = np.eye(m)         # 人工变量列
    T[:m, -1] = b
    T[m, n + m:n_tot] = 1.0                # Phase I 目标: min 1^T a
    basis = list(range(n + m, n_tot))      # 初始基 = 人工变量
    for i in range(m):                     # 消元: 基变量在目标行系数置 0
        T[m] = T[m] - T[i]

    iters_p1 = 0
    # ---------- Phase I: 求可行基 ----------
    while True:
        ent = -1
        for j in range(n_tot):             # Bland: 最小下标的负简约成本列入基
            if T[m, j] < -tol:
                ent = j
                break
        if ent < 0:
            break
        ratio = np.full(m, np.inf)
        for i in range(m):                 # 比值检验定出基行
            if T[i, ent] > tol:
                ratio[i] = T[i, -1] / T[i, ent]
        lev = int(np.argmin(ratio))
        if np.isinf(ratio[lev]):           # 所有系数 <= 0: 目标可无限减小 -> 无界
            return {"status": "unbounded", "iterations": iters_p1, "iters_p1": iters_p1}
        rmin = ratio[lev]
        cands = [i for i in range(m) if abs(ratio[i] - rmin) < 1e-8]
        lev = min(cands, key=lambda i: basis[i])   # Bland 平局规则
        _pivot(T, basis, lev, ent)
        iters_p1 += 1

    if -T[m, -1] > 1e-7:                   # Phase I 最优值 = 人工变量之和 > 0 -> 不可行
        return {"status": "infeasible", "iterations": iters_p1, "iters_p1": iters_p1}
    # 退化情形: 若有人工变量仍留在基中(取值必为 0), 把它换出或删除冗余行
    rows_to_drop = []
    for i in range(m):
        if basis[i] >= n + m:              # 该行基变量是人工变量
            col = -1
            for j in range(n + m):         # 找非人工列的非零元作枢轴
                if abs(T[i, j]) > tol:
                    col = j
                    break
            if col >= 0:
                _pivot(T, basis, i, col)
                iters_p1 += 1
            else:                          # 整行系数为 0: 冗余约束, 删除该行
                rows_to_drop.append(i)
    if rows_to_drop:
        keep = [i for i in range(m) if i not in rows_to_drop]
        T = np.vstack([T[keep], T[m:]])
        basis = [basis[i] for i in keep]
        m = len(basis)
    # ---------- Phase II: 删除人工列, 恢复原目标 ----------
    T = np.hstack([T[:, :n + m], T[:, -1:]])
    T[-1] = 0.0
    T[-1, :n] = d                          # 原目标系数(松弛变量成本为 0)
    for i in range(m):                     # 消元: 基变量在目标行系数置 0
        jb = basis[i]
        if jb < n:
            T[-1] = T[-1] - d[jb] * T[i]
    path = []
    iters_p2 = 0
    while True:
        x_now = np.zeros(n)                # 记录当前顶点坐标(画单纯形路径图用)
        for i in range(m):
            if basis[i] < n:
                x_now[basis[i]] = T[i, -1]
        path.append(x_now.copy())
        ent = -1
        for j in range(n + m):             # Bland: 找最小下标的负简约成本列
            if T[-1, j] < -tol:
                ent = j
                break
        if ent < 0:                        # 所有简约成本 >= 0: 已最优
            break
        ratio = np.full(m, np.inf)
        for i in range(m):
            if T[i, ent] > tol:
                ratio[i] = T[i, -1] / T[i, ent]
        lev = int(np.argmin(ratio))
        if np.isinf(ratio[lev]):           # 比值检验失败 -> 无界
            return {"status": "unbounded", "iterations": iters_p1 + iters_p2,
                    "iters_p1": iters_p1, "iters_p2": iters_p2}
        rmin = ratio[lev]
        cands = [i for i in range(m) if abs(ratio[i] - rmin) < 1e-8]
        lev = min(cands, key=lambda i: basis[i])
        _pivot(T, basis, lev, ent)
        iters_p2 += 1
    x = np.zeros(n)
    for i in range(m):                     # 从最优表读解
        if basis[i] < n:
            x[basis[i]] = T[i, -1]
    z = T[-1, -1] if maximize else -T[-1, -1]  # 目标行 RHS = -min 值 = max 值
    return {"x": x, "z": z, "status": "optimal", "iterations": iters_p1 + iters_p2,
            "iters_p1": iters_p1, "iters_p2": iters_p2,
            "basis": basis.copy(), "path": path, "tableau": T, "n_slack": m}


def basis_names(basis, n):
    """把基变量下标翻译成可读名称: 下标 < n 为 x, 否则为松弛变量 s。"""
    return [("x%d" % (j + 1)) if j < n else ("s%d" % (j - n + 1)) for j in basis]

# ========== 3. 手写单纯形法求解 + scipy.optimize.linprog 对照 ==========
# scipy 的 linprog 只接受 min 形式, 故 max c^T x 传入 c=-C;
# 返回的 res.fun 是 min 问题的最优值 = -z*, 需再取负号还原为利润。
t0 = time.perf_counter()
sol = two_phase_simplex(C, A, B, maximize=True)
t_simplex = time.perf_counter() - t0

t0 = time.perf_counter()
res = linprog(c=-C, A_ub=A, b_ub=B, bounds=[(0, None)] * len(C), method="highs")
t_scipy = time.perf_counter() - t0

print("=" * 62)
print("【指标 1】最优解与最优目标函数值")
print("  手写单纯形法: x* = %s,  z* = %.10g,  status = %s,  Phase I %d 次 / Phase II %d 次换基"
      % (np.round(sol["x"], 6), sol["z"], sol["status"], sol["iters_p1"], sol["iters_p2"]))
print("  Phase II 顶点迭代路径: %s" % [np.round(p, 4).tolist() for p in sol["path"]])
print("  最优基(变量): %s" % basis_names(sol["basis"], len(C)))
print("  scipy linprog: x* = %s,  z* = %.10g,  status = %d (%s)"
      % (np.round(res.x, 6), -res.fun, res.status, res.message.strip()))
print("  两种方法对照: x 最大偏差 = %.3e,  最优值偏差 = %.3e"
      % (np.max(np.abs(res.x - sol["x"])), abs(-res.fun - sol["z"])))
print("【指标 7】求解时间:  手写单纯形 %.3f ms | scipy(HiGHS) %.3f ms"
      % (t_simplex * 1e3, t_scipy * 1e3))

# ========== 4. 对偶问题与影子价格（三条途径交叉核对） ==========
# 对偶问题: min w = 100*y1 + 150*y2 + 80*y3  s.t.  A^T y >= c,  y >= 0
# 理论: 最优单纯形表目标行中"松弛变量列"的系数就是影子价格 y_i (i = 1..m)
Tf = sol["tableau"]
m, n = len(B), len(C)
y_tableau = Tf[-1, n:n + m]                              # 途径1: 从最终表读出
if hasattr(res, "ineqlin") and hasattr(res.ineqlin, "marginals"):
    # 途径2: linprog 解的是 min -c^T x, scipy 返回的 ineqlin.marginals 满足
    #        d(min)/db_i = -marginals_i, 故 max 问题的影子价格 = -marginals
    y_scipy = -res.ineqlin.marginals
else:
    y_scipy = None
res_dual = linprog(c=B, A_ub=-A.T, b_ub=-C,             # 途径3: 直接解对偶问题
                   bounds=[(0, None)] * m, method="highs")
y_dual = res_dual.x
print("=" * 62)
print("【指标 3】影子价格(对偶变量) 三条途径核对")
print("  途径1 最优表读出: y = %s" % np.round(y_tableau, 6))
print("  途径2 scipy边际值取负: y = %s" % (np.round(y_scipy, 6) if y_scipy is not None else "当前版本不可用"))
print("  途径3 解对偶问题: y = %s" % np.round(y_dual, 6))
for k, name in enumerate(RES_NAMES):
    print("  %s: y%d = %.6f  -> 该资源每增加 1 单位, 最优利润约增加 %.2f 万元"
          % (name, k + 1, y_tableau[k], y_tableau[k]))
print("【指标 2】对偶间隙: 对偶最优值 w* = %.10g,  原问题最优值 z* = %.10g,  gap = w* - z* = %.3e"
      % (res_dual.fun, sol["z"], res_dual.fun - sol["z"]))
# 互补松弛验证: 最优时 (b_i - a_i^T x*) * y_i = 0, 即"有剩余的资源影子价格必为 0"
s_use = A @ sol["x"]
s_rem = B - s_use
print("  互补松弛验证: 资源用量 = %s, 剩余 = %s" % (np.round(s_use, 4), np.round(s_rem, 4)))
print("    -> 原料B 剩余 %.1f kg > 0, 故 y2 = 0; 原料A 与工时恰好用完, 故 y1, y3 > 0" % s_rem[1])

# ========== 5. 灵敏度分析: 目标系数 c 与右端项 b 的允许变化区间 ==========
# 由最优表可得 B^{-1} = 最优表约束行中"初始单位阵列"(松弛列)对应的列
Binv = Tf[:m, n:n + m]
# (a) 右端项 b 的区间: 基保持可行  <=>  B^{-1} b >= 0
b_ranges = []
for i in range(m):
    e = np.zeros(m); e[i] = 1.0
    alpha = Binv @ e                       # b_i 的单位变化对基变量取值的影响
    beta = Binv @ B                        # 当前基变量取值
    lo, hi = -np.inf, np.inf
    for k in range(m):
        if abs(alpha[k]) > 1e-9:
            bound = -beta[k] / alpha[k]
            if alpha[k] > 0:
                lo = max(lo, bound)
            else:
                hi = min(hi, bound)
    b_ranges.append((B[i] + lo, B[i] + hi))
# (b) 目标系数 c 的区间: 最优基不变 <=> 所有非基变量简约成本 >= 0 (min 形式)
d = -C                                     # min 形式的目标系数
d_b = np.array([d[j] if j < n else 0.0 for j in sol["basis"]])
c_ranges = []
for j in range(n):
    lo, hi = -np.inf, np.inf
    pos = sol["basis"].index(j) if j in sol["basis"] else None
    for k in range(n + m):
        if k in sol["basis"]:
            continue
        col = Tf[:m, k]                    # 最终表第 k 列 = B^{-1} a_k
        rk = (d[k] if k < n else 0.0) - d_b.dot(col)   # 当前简约成本(>= 0)
        if pos is not None:                # c_j 变化会牵动所有非基列检验数
            coeff = col[pos]
            if abs(coeff) > 1e-9:
                if coeff > 0:
                    lo = max(lo, -rk / coeff)
                else:
                    hi = min(hi, -rk / coeff)
        elif k == j:                       # c_j 非基: 只影响自己那一列
            lo = max(lo, -rk)
    c_ranges.append((C[j] + lo, C[j] + hi))
print("=" * 62)
print("【指标 4】灵敏度区间(最优基不变的前提下)")
for j in range(n):
    print("  目标系数 c%d ∈ [%.4f, %.4f]  (当前值 %.1f)"
          % (j + 1, c_ranges[j][0], c_ranges[j][1], C[j]))
for i in range(m):
    print("  右端项 b%d(%s) ∈ [%.4f, %s]  (当前值 %.1f)"
          % (i + 1, RES_NAMES[i], b_ranges[i][0],
             "%.4f" % b_ranges[i][1] if np.isfinite(b_ranges[i][1]) else "+inf", B[i]))
# 简约成本: r_j = c_j - y^T a_j (max 形式; 最优时所有 r_j <= 0)
print("【指标 5】简约成本 r_j = c_j - y^T a_j  (max 形式, 最优时应 <= 0)")
for j in range(n):
    rj = C[j] - y_tableau.dot(A[:, j])
    print("  x%d: r = %.4f" % (j + 1, rj))
for i in range(m):
    print("  s%d(资源%s): r = %.4f   (注: 松弛变量简约成本 = -影子价格)"
          % (i + 1, RES_NAMES[i], -y_tableau[i]))

# ========== 6. 可行域顶点枚举(验证"最优解在顶点")与状态判定 ==========
def enumerate_vertices(A, b, n):
    """枚举二维可行域全部顶点: 两两联立边界直线(含 x1=0, x2=0)求交点并检验可行性。"""
    lines = np.vstack([A, np.eye(n)])      # 约束边界: A x = b 与坐标轴 x = 0
    rhs = np.hstack([b, np.zeros(n)])
    verts = []
    for (i, j) in combinations(range(len(lines)), 2):
        M = np.vstack([lines[i], lines[j]])
        if abs(np.linalg.det(M)) < 1e-9:   # 两线平行/重合, 无唯一交点
            continue
        p = np.linalg.solve(M, rhs[[i, j]])
        if np.all(p >= -1e-8) and np.all(A @ p <= b + 1e-8):  # 可行性检验
            if all(np.linalg.norm(p - q) > 1e-6 for q in verts):   # 去重
                verts.append(p)
    return verts

verts = enumerate_vertices(A, B, n)
verts = sorted(verts, key=lambda v: (v[0], v[1]))
print("=" * 62)
print("【灵敏度相关】可行域全部顶点及各顶点目标值")
for v in verts:
    print("  顶点 (%6.2f, %6.2f)  ->  z = %.1f" % (v[0], v[1], C.dot(v)))
print("  验证: 最优值 %.1f 在顶点 (%s) 取得, 与单纯形法一致"
      % (sol["z"], np.round(sol["x"], 4)))
print("【指标 6】状态判定: status = %s (feasible & bounded -> 存在最优解)" % sol["status"])
# 退化性检查: 最优基中是否有基变量取 0
degenerate = any(abs(Tf[i, -1]) < 1e-8 for i in range(m))
print("  退化性检查: 最优基变量取值 = %s, %s" % (np.round(Tf[:m, -1], 4), "存在退化" if degenerate else "无退化"))
# 状态判定演示: 不可行与无界两种异常情况
sol_inf = two_phase_simplex(np.array([1.0, 1.0]),
                            np.array([[1.0, 1.0], [-1.0, -1.0]]),
                            np.array([5.0, -10.0]), maximize=True)   # x1+x2<=5 且 x1+x2>=10
res_inf = linprog(c=[-1, -1], A_ub=[[1, 1], [-1, -1]], b_ub=[5, -10],
                  bounds=[(0, None)] * 2, method="highs")
sol_unb = two_phase_simplex(np.array([1.0, 0.0]),
                            np.array([[0.0, 1.0]]),
                            np.array([5.0]), maximize=True)          # 只有 x2<=5, x1 可无限增大
res_unb = linprog(c=[-1, 0], A_ub=[[0, 1]], b_ub=[5],
                  bounds=[(0, None)] * 2, method="highs")
print("  不可行示例: 手写 status = %s | scipy status = %d -> 两者一致"
      % (sol_inf["status"], res_inf.status))
print("  无界示例  : 手写 status = %s | scipy status = %d -> 两者一致"
      % (sol_unb["status"], res_unb.status))

# ========== 7. 稍大实例: 5 变量生产计划(演示多变量求解) ==========
C5 = np.array([3.0, 2.0, 4.0, 5.0, 1.0])   # 5 种产品单位利润
A5 = np.array([[2.0, 1.0, 3.0, 2.0, 1.0],  # 原料A 消耗
               [1.0, 3.0, 1.0, 4.0, 2.0],  # 原料B 消耗
               [3.0, 1.0, 2.0, 1.0, 3.0],  # 工时消耗
               [1.0, 2.0, 2.0, 3.0, 1.0]]) # 设备机时消耗
B5 = np.array([40.0, 50.0, 60.0, 45.0])
print("=" * 62)
print("【稍大实例】5 变量生产计划: max z = 3x1+2x2+4x3+5x4+x5")
sol5 = two_phase_simplex(C5, A5, B5, maximize=True)
res5 = linprog(c=-C5, A_ub=A5, b_ub=B5, bounds=[(0, None)] * 5, method="highs")
print("  手写单纯形: x* = %s, z* = %.6f, Phase I %d 次 / Phase II %d 次换基, 基 = %s"
      % (np.round(sol5["x"], 6), sol5["z"], sol5["iters_p1"], sol5["iters_p2"],
         basis_names(sol5["basis"], 5)))
print("  scipy linprog: x* = %s, z* = %.6f" % (np.round(res5.x, 6), -res5.fun))
print("  对照: x 最大偏差 = %.3e, z 偏差 = %.3e"
      % (np.max(np.abs(res5.x - sol5["x"])), abs(-res5.fun - sol5["z"])))
y5 = sol5["tableau"][-1, 5:5 + 4]          # 5 变量实例的影子价格
print("  影子价格 y = %s" % np.round(y5, 6))
for j in range(5):
    rj = C5[j] - y5.dot(A5[:, j])
    tag = "  <- 非基, |r| 即强制生产的单位利润损失" if sol5["x"][j] < 1e-8 else "  (基变量)"
    print("  简约成本 r_x%d = %+.4f%s" % (j + 1, rj, tag))
print("  资源使用: 用量 = %s, 剩余 = %s"
      % (np.round(A5 @ sol5["x"], 4), np.round(B5 - A5 @ sol5["x"], 4)))

# ========== 8. 四张可视化图(对应第四节) ==========
def line1(x): return 100 - 2 * x           # 2x1 + x2 = 100
def line2(x): return (150 - x) / 2         # x1 + 2x2 = 150
def line3(x): return 80 - x                # x1 + x2 = 80

poly = np.array([[0, 0], [50, 0], [20, 60], [10, 70], [0, 75]])   # 可行域顶点(逆时针)

# ---- 图1: 二维可行域 + 目标等值线 + 最优点 ----
fig, ax = plt.subplots(figsize=(7.4, 6.2))
ax.fill(poly[:, 0], poly[:, 1], color="#cde2fb", alpha=0.85, label="可行域 S (凸多边形)")
ax.plot(np.r_[poly[:, 0], poly[0, 0]], np.r_[poly[:, 1], poly[0, 1]],
        color="#2a78d6", lw=2.0, alpha=0.9)
for fn, col in [(line1, "#eb6834"), (line2, "#1baf7a"), (line3, "#4a3aa7")]:
    xs = np.linspace(0, 62, 400)
    ys = fn(xs)
    msk = ys >= -1
    ax.plot(xs[msk], ys[msk], color=col, lw=1.8, zorder=2)
# 目标等值线: 3x1 + 2x2 = k (k = 60, 120, 180), 沿梯度方向目标值增大
for k in [60, 120, 180]:
    xs = np.linspace(0, 62, 2)
    ys = (k - 3 * xs) / 2
    msk = ys >= -1
    ax.plot(xs[msk], ys[msk], "--", color="#8a8a8a", lw=1.2, zorder=1)
    ax.text(xs[msk][-1] + 0.5, ys[msk][-1], "z=%d" % k, fontsize=9, color="#555555")
ax.arrow(38, 38, 9, 6, head_width=1.6, head_length=2.4, fc="#333333", ec="#333333",
         length_includes_head=True, zorder=4)
ax.text(50, 26, "梯度方向\nc=(3,2)", fontsize=9, ha="center")
ax.text(16, 73, "2x1 + x2 = 100\n(原料A)", fontsize=9, color="#c2541e", rotation=-63)
ax.text(6, 67, "x1 + 2x2 = 150\n(原料B)", fontsize=9, color="#0e8a5f", rotation=-27)
ax.text(30, 53, "x1 + x2 = 80\n(工时)", fontsize=9, color="#3a2e8c", rotation=-45)
ax.plot(*sol["x"], marker="*", ms=17, color="#e34948", zorder=5, label="最优解 (20, 60)")
ax.annotate("最优解 x* = (20, 60)\nz* = 180", xy=sol["x"], xytext=(33, 78),
            fontsize=10, arrowprops=dict(arrowstyle="->", color="#e34948"), color="#e34948")
for v in verts:
    ax.plot(v[0], v[1], "o", ms=6, color="#2a78d6", zorder=3)
    ax.annotate("(%g, %g)" % (v[0], v[1]), xy=v, xytext=(v[0] + 1.5, v[1] - 6),
                fontsize=8, color="#1c5cab")
ax.set_xlim(-1, 63); ax.set_ylim(-2, 92)
ax.set_xlabel("x1 (产品P1产量/件)"); ax.set_ylabel("x2 (产品P2产量/件)")
ax.set_title("图1 可行域、目标等值线与最优点")
ax.legend(loc="lower left", fontsize=9)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, "lp_feasible_region.png"), dpi=200, bbox_inches="tight")
plt.close(fig)

# ---- 图2: 单纯形法顶点迭代路径 ----
fig, ax = plt.subplots(figsize=(7.4, 6.2))
ax.fill(poly[:, 0], poly[:, 1], color="#cde2fb", alpha=0.85)
ax.plot(np.r_[poly[:, 0], poly[0, 0]], np.r_[poly[:, 1], poly[0, 1]],
        color="#2a78d6", lw=2.0)
for fn, col in [(line1, "#eb6834"), (line2, "#1baf7a"), (line3, "#4a3aa7")]:
    xs = np.linspace(0, 62, 400)
    ax.plot(xs, np.maximum(fn(xs), -1), color=col, lw=1.8)
path = np.array(sol["path"])
ax.plot(path[:, 0], path[:, 1], "-o", color="#eb6834", lw=2.0, ms=8, zorder=4,
        label="单纯形迭代路径")
ax.annotate("初始顶点 (10, 70)\n(Phase I 求得可行基)\nz = 170", xy=path[0],
            xytext=(1, 56), fontsize=9, color="#c2541e",
            arrowprops=dict(arrowstyle="->", color="#c2541e"))
if len(path) > 1:
    ax.annotate("第1次换基\nz: 170 -> 180", xy=((path[0, 0] + path[1, 0]) / 2, (path[0, 1] + path[1, 1]) / 2),
                xytext=(26, 46), fontsize=9, color="#c2541e",
                arrowprops=dict(arrowstyle="->", color="#c2541e"))
ax.plot(*sol["x"], marker="*", ms=17, color="#e34948", zorder=5, label="最优解 (20, 60)")
ax.set_xlim(-1, 63); ax.set_ylim(-2, 92)
ax.set_xlabel("x1 (产品P1产量/件)"); ax.set_ylabel("x2 (产品P2产量/件)")
ax.set_title("图2 单纯形法顶点迭代路径: (10,70) -> (20,60)")
ax.legend(loc="lower left", fontsize=9)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, "lp_simplex_path.png"), dpi=200, bbox_inches="tight")
plt.close(fig)

# ---- 图3: 右端项 b1 变化时最优值变化曲线(分段线性, 展示影子价格转折) ----
b1_grid = np.arange(70.0, 176.0, 0.5)
z_curve = []
for b1v in b1_grid:
    b_tmp = B.copy(); b_tmp[0] = b1v
    r = linprog(c=-C, A_ub=A, b_ub=b_tmp, bounds=[(0, None)] * n, method="highs")
    z_curve.append(-r.fun)
z_curve = np.array(z_curve)
slopes = np.diff(z_curve) / np.diff(b1_grid)          # 相邻两点间斜率 = 该区段影子价格
kinks = [float(b1_grid[i]) for i in range(1, len(slopes)) if abs(slopes[i] - slopes[i - 1]) > 1e-3]
print("=" * 62)
print("【图3 支持】b1 在 [70, 175] 内最优值曲线的转折点: b1 = %s" % kinks)
fig, ax = plt.subplots(figsize=(7.6, 5.4))
ax.plot(b1_grid, z_curve, color="#2a78d6", lw=2.2, label="最优值 z*(b1)")
for kb in kinks:
    zi = z_curve[np.argmin(np.abs(b1_grid - kb))]
    ax.plot(kb, zi, "o", ms=8, color="#eb6834", zorder=4)
    ax.annotate("转折点 b1=%g" % kb, xy=(kb, zi), xytext=(kb - 18, zi + 8),
                fontsize=9, color="#c2541e",
                arrowprops=dict(arrowstyle="->", color="#c2541e"))
ax.axvline(100, color="#8a8a8a", ls="--", lw=1.2)
ax.text(101.5, 156, "当前 b1 = 100\n影子价格 y1 = 1 (斜率)", fontsize=9)
# 各段斜率标注(斜率 = 该区段的影子价格)
for (s, e, lbl) in [(70, 75, "斜率=2"), (82, 87, "斜率=4/3"), (115, 125, "斜率=1"),
                    (165, 172, "斜率=0\n(资源不再稀缺)")]:
    i0 = np.argmin(np.abs(b1_grid - s)); i1 = np.argmin(np.abs(b1_grid - e))
    ax.text(b1_grid[(i0 + i1) // 2], z_curve[(i0 + i1) // 2] + 4, lbl, fontsize=9,
            ha="center", color="#1c5cab")
ax.set_xlabel("右端项 b1 (原料A 拥有量/kg)")
ax.set_ylabel("最优利润 z* (万元)")
ax.set_title("图3 右端项 b1 变化时的最优值曲线(分段线性)")
ax.legend(loc="upper left", fontsize=9)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, "lp_b_curve.png"), dpi=200, bbox_inches="tight")
plt.close(fig)

# ---- 图4: 各资源影子价格柱状图 ----
fig, ax = plt.subplots(figsize=(7.6, 4.8))
y_all = np.r_[y_tableau[:3]]
colors = ["#2a78d6" if v > 1e-9 else "#c9c9c9" for v in y_all]
bars = ax.bar(RES_NAMES, y_all, color=colors, width=0.55, edgecolor="#1c5cab", lw=1.0)
for bbar, v in zip(bars, y_all):
    ax.text(bbar.get_x() + bbar.get_width() / 2, v + 0.03, "%.2f" % v, ha="center",
            fontsize=12, color="#0b0b0b")
ax.axhline(0, color="#8a8a8a", lw=1.0)
ax.text(1, 0.5, "原料B 影子价格 = 0:\n剩余 10 kg, 增加投入\n不改变最优利润", fontsize=9,
        ha="center", color="#52514e")
ax.set_ylim(0, 1.25)
ax.set_ylabel("影子价格 y_i (万元/单位资源)")
ax.set_title("图4 各资源影子价格(资源边际价值)柱状图")
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, "lp_shadow_price_bar.png"), dpi=200, bbox_inches="tight")
plt.close(fig)

# 交互式环境下弹出全部窗口; 无图形界面(如 MPLBACKEND=Agg)时静默跳过
with warnings.catch_warnings():
    warnings.simplefilter("ignore", UserWarning)
    plt.show()
print("=" * 62)
print("四张图已保存到 figures/ 目录: lp_feasible_region.png, lp_simplex_path.png, lp_b_curve.png, lp_shadow_price_bar.png")

七、结果解读与注意事项

7.1 以第六节代码输出为例的完整解读

运行第六节代码，控制台输出 【问题】、【指标 1】~【指标 7】、【灵敏度相关】、【稍大实例】、【图3 支持】 各节。以下按论文叙述顺序把每个数字「翻译」成人话。

(1) 最优解与最优值（指标 1）

手写单纯形法与 scipy（HiGHS）都给出 $x^* = (20, 60)$、$z^* = 180$，两种方法偏差为 0——最优生产方案是每天生产 P1 20 件、P2 60 件，最大利润 180 万元。手写算法 Phase I 换基 3 次（找初始可行基）、Phase II 换基 1 次（从初始顶点 $(10, 70)$ 走到最优顶点 $(20, 60)$，利润 $170 \to 180$）。最优基为 $\{x_1, s_2, x_2\}$：基变量是两种产品产量和原料B 的松弛量，说明卡脖子的约束是原料A 与工时（对应两个非基松弛变量 $s_1 = s_3 = 0$），原料B 的松弛量 $s_2 = 10$ 留在基里——原料B 每天用掉 $20 + 120 = 140$ kg，剩余 10 kg。

(2) 影子价格（指标 3）与互补松弛

三条途径（最优表读出、HiGHS 边际值取负、直接解对偶问题）一致给出 $y = (1, 0, 1)$：

• $y_1 = 1$：原料A 每多 1 kg，最优利润增加 1 万元（例如原料A 变成 101 kg 时，最优解变为 $(21, 59)$，利润 181 万元）；
• $y_2 = 0$：原料B 每天剩余 10 kg 用不完，白送 1 kg 也不改变利润——互补松弛的体现（剩余量与影子价格至少一个为 0）；
• $y_3 = 1$：工时每多 1 小时，最优利润增加 1 万元（80 小时全用完）。

决策含义：若要扩充产能，原料A 与工时的投入产出比都是 1 万元/单位，原料B 不必购买——这就是「影子价格指导决策」的标准叙事。

(3) 对偶间隙（指标 2）

对偶问题最优值 $w^* = 180$ 与原问题最优值 $z^* = 180$ 完全相等，gap $= 0$。论文写法：「由强对偶定理，对偶间隙为 0，故所得解为全局最优解。」这是线性规划相对启发式算法的「降维打击」：最优性是被数学定理保证的。

(4) 灵敏度区间（指标 4）

• 目标系数：$c_1 \in [2, 4]$（P1 利润在 24 万元/件之间波动时，最优方案仍是每天 20+60 件）；$c_2 \in [1.5, 3]$（P2 利润在 1.53 万元/件之间波动时同理）。区间外最优解会变：例如 P2 利润涨到 3.5 万元/件，最优解将移动到 $(10, 70)$。
• 右端项：$b_1 \in [90, 160]$（原料A 在 90160 kg 内波动，最优基不变、影子价格恒为 1）；$b_2 \in [140, +\infty)$（原料B 不少于 140 kg 即可，再多都是剩余——$+\infty$ 上界正是「资源不稀缺」的数学表达）；$b_3 \in [50, 83.33]$（工时 5083.33 小时，影子价格恒为 1）。其中 $b_1$ 区间的两个端点 90 与 160 正是图 3 中最优值曲线斜率变化的转折点（见下）；$b_2$ 的上界为 $+\infty$，表示原料B 再充足也不会改变最优基。

(5) 简约成本（指标 5）

主实例中两个决策变量都是基变量，$r_{x_1} = r_{x_2} = 0$（正在生产的产品，利润已被充分挖掘）；三个松弛变量的简约成本分别为 $-1, 0, -1$，恰好等于各自影子价格的相反数。5 变量实例中 $r_{x_2} = -1.1667$（x2 是非基变量，产量为 0）：「强行生产 1 件 x2，总利润将损失 1.1667 万元」——回答「为什么不生产 x2」。

(6) 顶点枚举与状态判定（指标 6）

枚举可行域全部 5 个顶点：$(0,0)$、$(50,0)$、$(20,60)$、$(10,70)$、$(0,75)$，目标值分别为 0、150、180、170、150——最优值 180 在顶点 $(20, 60)$ 取得，用最朴素的方式验证了顶点最优定理。退化性检查：最优基变量取值 $(20, 10, 60)$ 全为正，无退化，最优解唯一。两个异常示例中，手写实现与 scipy 一致地给出 infeasible（status=2，$x_1 + x_2 \le 5$ 与 $x_1 + x_2 \ge 10$ 矛盾）与 unbounded（status=3，只约束 $x_2$ 而 $x_1$ 可无限增大、利润可无限提高——典型漏约束）。

(7) 5 变量实例

5 产品 4 资源问题的最优解 $x^* = (10, 0, 0, 10, 0)$、$z^* = 80$：只生产 x1（10 件）与 x4（10 件），其余不生产（简约成本为负，强产必亏）；原料A、原料B 恰好用完（用量 40/50，剩余 0/0），工时与机时分别剩余 20、5——对应影子价格 $y = (1.1667, 0.6667, 0, 0)$：只有前两种资源值钱。手写单纯形与 scipy 结果一致（偏差 $10^{-15}$ 量级）。这个例子说明：变量再多，LP 的「解法结构」完全一样，而影子价格/简约成本自动告诉你「该生产什么、该买什么资源」。

(8) 四张图的解读

• 图1（可行域 + 等值线 + 最优点）：灰色虚线是等值线 $z = 60, 120, 180$，沿箭头（梯度 $c = (3,2)$）方向利润增大；最外面的等值线 $z = 180$ 恰好「擦过」可行域顶点 $(20, 60)$——几何上，最优解就是「等值线沿梯度平移、最后离开可行域的那个点」。$(20, 60)$ 是原料A 直线与工时直线的交点，两条资源同时用完，与 $y_1 = y_3 = 1 > 0$ 呼应。
• 图2（单纯形路径）：橙色路径显示算法从顶点 $(10, 70)$（Phase I 求得，$z = 170$）沿「原料B 约束边」走到 $(20, 60)$（$z = 180$），一步到位。路径始终沿可行域边界、目标值单调上升——单纯形法就是「顶点间的爬山」。
• 图3（b1 变化曲线）：$z^*(b_1)$ 是分 4 段的分段线性函数，转折点 $b_1 = 75, 90, 160$ 标记着最优基的三次更替，其中 90 与 160 正是当前最优基灵敏度区间 $[90, 160]$ 的两个端点；$[90, 160]$ 段斜率 $= 1 = y_1$（区间内影子价格不变）；$b_1 \in [75, 90]$ 时原料B 的约束接管（与原料A 一起卡脖子，斜率变为 $4/3$），$b_1 < 75$ 时最优解退到「只生产 P2」（斜率 2）；$b_1 > 160$ 时原料A 不再稀缺，瓶颈转移到工时（斜率 0）。竖虚线标出当前工作点 $b_1 = 100$。这张图把「影子价格 = 导数、灵敏度区间 = 线性段」讲得明明白白。
• 图4（影子价格柱状图）：蓝色柱（原料A、工时，高度 1）是瓶颈资源，灰色柱（原料B，高度 0）是冗余资源——一眼看出「钱该往哪花」。

7.2 常见坑（务必避开）

1. max / min 没有换算就丢给求解器。scipy.optimize.linprog 只做最小化：max 问题要传 c = -c_original，且返回的 fun 要取负号才是原问题的最大值。忘了取负号，最优值符号就反了。
2. 约束方向写反。$\ge$ 型约束要整体乘 $-1$ 变成 $\le$（A_ub, b_ub 约定为「小于等于」）；把「$x_1 + x_2 \ge 10$」写成「$\le 10$」会让模型漏掉一半可行域甚至报不可行。程序报 infeasible 时，第一件事就是逐条检查约束方向与数值。
3. 把影子价格当「全局导数」用。影子价格只在灵敏度区间内成立。原料A 翻倍（100 → 200 kg）时，$y_1 = 1$ 已失效（图 3 显示 160 kg 之后斜率变为 0），利润绝不会「线性翻倍」。外推影子价格是灵敏度分析最常见的误用。
4. 把 $y_i = 0$ 误读为「资源没用」。$y_i = 0$ 只说明「当前再增加该资源不会改进目标」——可能是它真的冗余（如本例原料B），也可能是另一项资源卡得更死导致它暂时用不完。判断方法：看松弛量是否为正（互补松弛）。
5. 连续解直接四舍五入。LP 给出 $x_1 = 3.7$ 件时，取 3 或 4 都不保证可行或最优（可能违反约束，也可能错过更好的整数方案）。需要整数解时老老实实用整数规划（08 篇的 scipy.optimize.milp），并在论文中说明「LP 松弛给出上界，整数解与其差距为 xx」。
6. 把无界当喜讯。报 unbounded 意味着「利润可以无限大」——数学上成立、经济上荒谬，一定是漏写了约束（忘了资源限制或非负条件）。同理，报 infeasible 意味着约束矛盾或写反。两种状态都是建模 bug 的信号，先修模型再谈求解。
7. 忽视退化与多解。若最优基变量中有 0（退化），最优解往往不唯一（整条边都是最优）；论文中应写明「最优解不唯一，本文给出其中一个」或利用多解空间追求次级目标（如选择库存最小的最优解）。手写单纯形法遇退化可能循环，务必使用 Bland 规则。
8. 单位混乱导致数值病态。系数相差 $10^6$ 倍（如「元」与「亿元」混用）会使求解器数值不稳。统一单位（本文全部用「万元」「kg」「小时」）再建模。
9. 灵敏度区间解读越界。区间只保证「最优基（解的结构）不变」：$c_j$ 区间内最优解数值不变，但最优值会按 $x_j^* \Delta c_j$ 变化；$b_i$ 区间内最优值按 $y_i \Delta b_i$ 变化但解本身会变。两种「不变」的对象不同，论文中别写混。
10. 只用求解器、不交叉验证。竞赛论文建议「手写/自编求解 + 求解器对照」（如本文两阶段单纯形 vs HiGHS）或「原问题与对偶问题互相验证（gap = 0）」，这是展示算法功底的加分项，也能防住 c 取负号之类的手滑错误。

7.3 竞赛论文写作建议（话术模板）

优化类论文的标准结构是「建模 → 求解 → 验证 → 分析」，线性规划的段落建议按下面的话术组织，直接可用的模板句：

建模句：「设决策变量 $x_1, x_2$ 分别为产品 P1、P2 的日产量。以日利润最大为目标、以原料与工时限制为约束，建立线性规划模型：$\max z = 3x_1 + 2x_2$，$\text{s.t. } 2x_1 + x_2 \le 100, x_1 + 2x_2 \le 150, x_1 + x_2 \le 80, x_1, x_2 \ge 0$。模型中目标函数与全部约束均为决策变量的线性函数，故问题为线性规划，可求全局最优解。」

求解句：「采用两阶段单纯形法求解，并与 HiGHS 求解器结果交叉验证，两者一致：最优方案为 $x^* = (20, 60)$，即每天生产 P1 20 件、P2 60 件，最大利润 $z^* = 180$ 万元。进一步由强对偶定理，对偶问题最优值与原问题相等、对偶间隙为 0，证明所得解为全局最优。」

影子价格句：「求解得各资源影子价格为 $(1, 0, 1)$：原料 A 与工时的边际价值均为 1 万元/单位，且二者在最优方案中恰好用完，是产能瓶颈；原料 B 每天剩余 10 kg，影子价格为 0。因此建议优先扩充原料 A 采购量与工时，而非原料 B。」

灵敏度句：「灵敏度分析表明：当原料 A 拥有量在 $[90, 160]$ kg、工时在 $[50, 83.33]$ 小时、产品 P1 单位利润在 $[2, 4]$ 万元内波动时，最优解的结构保持不变；原料 B 不低于 140 kg 即可。当前参数距各区间边界均有较大余量，方案对数据波动具有稳健性。」

通用模板：「针对 ……（题目背景），以 …… 为目标、以 …… 为约束建立线性规划模型；使用单纯形法（HiGHS）求得全局最优解 ……；影子价格分析显示 …… 为瓶颈资源，灵敏度分析表明方案在参数波动 ±…% 内保持稳定。」

写作细节建议：① 模型公式务必完整写出（目标 + 每条约束 + 非负条件），这是优化论文的门面；② 求解后给出「最优基」或「哪些约束是紧的」，与影子价格呼应，体现理论深度；③ 图1（可行域图解）放在模型建立小节，图2（单纯形路径）放在求解小节，图3、图4（灵敏度与影子价格）放在结果分析小节；④ 若题目要求整数解，先解 LP 松弛并报告上界，再用 08 篇的整数规划方法，对比 gap 大小；⑤ 结论部分把最优方案翻译回「管理语言」：生产多少、买什么资源、赚多少钱，而不是只报数字。

八、延伸阅读

• 对偶单纯形法：单纯形法在「原问题可行、对偶不可行」的顶点间移动，对偶单纯形法反过来在「对偶可行、原不可行」的顶点间移动，逐步恢复原可行性。它的杀手锏场景是灵敏度/再优化：右端项 $b$ 变化后旧最优基不再可行，从旧表出发用对偶单纯形往往一两步就能回到最优，无需从头求解。整数规划的分支定界中每次加约束后也用它快速重优化（08 篇）。
• 内点法（Interior Point Method）：Karmarkar 于 1984 年提出的多项式时间算法，不从顶点走边，而是在可行域内部沿「中心路径」迭代逼近最优（牛顿法 + 对数障碍函数 $\min c^\top x - \mu \sum \ln x_j$）。理论复杂度优于单纯形法（$O(n^{3.5})$ 型），大规模稀疏问题上实践速度也占优；scipy 的 HiGHS 求解器内部正是单纯形法与内点法的混合实现（默认自动选择）。
• 运输问题的表上作业法：产销平衡运输问题（$m$ 个产地、$n$ 个销地、单位运费 $c_{ij}$）是 LP 的特殊结构，可用「西北角法 / 最小元素法 / Vogel 法」给出初始调运方案，再用「位势法」检验、闭回路法调整——全程只在运输表上操作，适合手工演示与竞赛中的小规模算例，能展现扎实的运筹学功底。
• 多目标线性规划：多个目标（利润、污染、就业……）难以同时最优时，常用三种转化：加权法 $\max \sum_k w_k c_k^\top x$（取不同权重画出 Pareto 前沿）、ε-约束法（把一个目标当约束 $\ge \varepsilon$，另一个当目标）、目标规划（设定各目标的期望值，最小化正负偏差量加权和 $\min \sum_k (w_k^+ d_k^+ + w_k^- d_k^-)$）。三者转化后仍可用本文的 LP 求解器。
• 经典文献：Dantzig (1947) 提出单纯形法；Karmarkar (1984) 提出内点法；Chvátal 的 Linear Programming 与《运筹学》（清华版）是竞赛常用的入门教材；HiGHS（Huangfu & Hall, 2018）是当前开源 LP/MILP 求解器的性能标杆，scipy ≥ 1.9 已内置。与本系列衔接：08 篇《整数规划》在本文模型上加取整约束并介绍分支定界与 scipy.optimize.milp。