【sklearn笔记】线性规划

1. 引言：什么是线性规划

线性规划（Linear Programming，LP）是运筹学中最基础、应用最广泛的优化模型。它的核心特征是：目标函数和所有约束条件都是决策变量的线性函数。尽管”线性”看似简单，但 LP 所蕴含的对偶理论、单纯形法、内点法等构成了现代数学规划的基石。

1975 年，Kantorovich 和 Koopmans 因资源配置理论（线性规划的核心应用）获得诺贝尔经济学奖。1984 年，Karmarkar 提出的内点法在 AT&T 公司的贝尔实验室引起轰动——它证明了 LP 可以在多项式时间内求解，挑战了统治近 40 年的单纯形法。

LP 的标准应用场景：假设你是一家工厂的生产经理，有多种产品可以使用有限的资源（原料、工时、机器）来生产，每种产品有不同利润，你需要决定各生产多少以最大化总利润。这个问题天然是线性的。

2. 线性规划问题的标准形式

2.1 标准形式（Standard Form）

LP 的标准形式是大多数求解器接受的输入格式：

$$\begin{aligned} \text{minimize} \quad & c^T x \\ \text{subject to} \quad & Ax = b \\ & x \geq 0 \end{aligned}$$

其中：

• x ∈ R^n：决策变量（decision variables）
• c ∈ R^n：成本系数向量（cost vector）
• A ∈ R^{m×n}：约束矩阵（constraint matrix），通常 m < n
• b ∈ R^m：右端常数向量（right-hand side）

重要：标准形式要求所有变量非负（x ≥ 0），且所有约束都是等式（Ax = b）。如果原始问题不满足这些条件，需要进行标准化。

2.2 将任意 LP 转化为标准形式

情形1：最大化问题

max c^T x  →  min -c^T x

情形2：不等式约束（≤）

a_i^T x ≤ b_i  →  a_i^T x + s_i = b_i,  s_i ≥ 0

引入松弛变量（slack variable） s_i。

情形3：不等式约束（≥）

a_i^T x ≥ b_i  →  a_i^T x - s_i = b_i,  s_i ≥ 0

引入剩余变量（surplus variable） s_i。

情形4：无约束变量（free variable）

x_j free  →  x_j = x_j^+ - x_j^-,   x_j^+ ≥ 0, x_j^- ≥ 0

将一个自由变量拆分为两个非负变量之差。

示例：将以下 LP 转化为标准形式：

max   3x_1 - 2x_2
s.t.  x_1 + x_2  ≤ 5
      2x_1 - x_2  ≥ 1
      x_1 ≥ 0, x_2 free

转化过程：

1. max → min：min -3x_1 + 2x_2
2. x_2 = x_2^+ - x_2^- (x_2^+, x_2^- ≥ 0)
3. 第一个约束加松弛变量 s_1 ≥ 0：x_1 + x_2^+ - x_2^- + s_1 = 5
4. 第二个约束加减剩余变量 s_2 ≥ 0：2x_1 - x_2^+ + x_2^- - s_2 = 1
5. 所有变量 ≥ 0：x_1, x_2^+, x_2^-, s_1, s_2 ≥ 0

2.3 规范形式（Canonical Form）

LP 还有另一种常见形式——规范形式（Canonical Form）：

$$\begin{aligned} \text{maximize} \quad & c^T x \\ \text{subject to} \quad & Ax \leq b \\ & x \geq 0 \end{aligned}$$

规范形式中约束全为 ≤ 不等式（针对最大化问题）。这种形式在对称对偶理论中非常方便。

2.4 几何直观

每个 LP 的可行域是一个凸多面体（convex polyhedron）——由有限个半空间的交集定义。对于二维 LP，可行域是一个凸多边形；三维则是一个凸多面体。

LP 基本定理：如果 LP 有最优解，那么至少有一个最优解出现在可行域的顶点（vertex/extreme point）上。这是单纯形法的基础——它只在顶点之间移动。

示例：考虑 LP

max   2x_1 + 3x_2
s.t.  x_1 + 2x_2 ≤ 8
      4x_1      ≤ 16
           4x_2 ≤ 12
      x_1, x_2  ≥ 0

可行域是由 (0,0), (0,3), (2,3), (4,0) 四个顶点围成的直角梯形。目标函数在顶点 (2,3) 处取得最大值 2×2 + 3×3 = 13。

3. 单纯形法（Simplex Method）

单纯形法是 Dantzig 于 1947 年提出的算法，即使在今天仍然是最广泛使用的 LP 求解方法之一。

3.1 核心思想

单纯形法从可行域的一个顶点出发，沿着多面体的边（edge）移动到相邻的顶点，每次移动都改善（或不恶化）目标函数值，直到达到最优解。

为什么只在顶点之间移动？根据 LP 基本定理，如果最优解存在，必在某个顶点取得。因此我们只需要检查顶点。

3.2 基本解与基本可行解

考虑标准形式 Ax = b, x ≥ 0，其中 A 是 m × n 矩阵，rank(A) = m。

将 A 的列分为两部分：

• 基矩阵 B：m 个线性无关的列组成的 m × m 可逆子矩阵
• 非基矩阵 N：剩下的 n - m 列

对应地将 x 分为基变量 x_B 和非基变量 x_N：

$$Ax = [B \quad N] \begin{bmatrix} x_B \\ x_N \end{bmatrix} = Bx_B + Nx_N = b$$

令非基变量 x_N = 0，则 x_B = B^{-1}b。这个解 (x_B, 0) 称为基本解。

如果 x_B ≥ 0，则称其为基本可行解（Basic Feasible Solution, BFS）。BFS 恰好对应可行域的顶点。

3.3 单纯形表（Simplex Tableau）

单纯形表是手工计算时的标准工具。以最大化问题为例，初始表如下：

         x_1    x_2    ...    x_n    |   RHS
------------------------------------+--------
row_0    -c_1   -c_2   ...   -c_n   |    0     (目标函数行)
row_1    a_11   a_12   ...   a_1n   |   b_1
row_2    a_21   a_22   ...   a_2n   |   b_2
...      ...    ...    ...   ...    |   ...
row_m    a_m1   a_m2   ...   a_mn   |   b_m

row_0 存储的是检验数（reduced cost）的相反数。当所有检验数 ≤ 0（最大化）或 ≥ 0（最小化）时，当前解最优。

3.4 单纯形法的步骤

步1 - 最优性检验：计算检验数 σ_j = c_j - z_j = c_j - c_B^T B^{-1} A_j。

• 最大化：若所有 σ_j ≤ 0 → 最优，停止
• 最小化：若所有 σ_j ≥ 0 → 最优，停止

步2 - 选择入基变量（Entering Variable）：

• 最大化：选 σ_j > 0 中最大的（Dantzig 规则）
• 也可以使用 Bland 规则（选择下标最小的正检验数）来避免循环

步3 - 选择出基变量（Leaving Variable）：用最小比值测试（Minimum Ratio Test）：

$$\theta = \min_i \left\{ \frac{x_{B_i}}{y_{ij}} : y_{ij} > 0 \right\}$$

其中 y_{ij} 是入基列中第 i 行的系数。如果不存在 y_{ij} > 0，则问题无界。

步4 - 转轴操作（Pivot）：以 y_{ij} 为主元做 Gauss-Jordan 消元，使入基列成为单位向量。

步5：返回步1。

3.5 手工计算示例

考虑问题：

max   5x_1 + 4x_2
s.t.  6x_1 + 4x_2 ≤ 24
      x_1 + 2x_2 ≤ 6
      x_1, x_2 ≥ 0

引入松弛变量 s_1, s_2：

max   5x_1 + 4x_2 + 0s_1 + 0s_2
s.t.  6x_1 + 4x_2 + s_1      = 24
      x_1 + 2x_2       + s_2  = 6
      x_1, x_2, s_1, s_2 ≥ 0

初始单纯形表：

          x_1   x_2   s_1   s_2  |  RHS
----------------------------------+------
z-row    -5    -4     0     0    |   0
s_1       6     4     1     0    |  24
s_2       1     2     0     1    |   6

当前 BFS：x_1=0, x_2=0, s_1=24, s_2=6, z=0。

迭代1：

• 入基：x_1（检验数 -5 最小，即 5 最大）
• 比值测试：s_1: 24/6=4, s_2: 6/1=6 → s_1 出基
• 主元：6（第一行，第一列）
• 执行转轴操作后：

          x_1   x_2        s_1        s_2  |  RHS
--------------------------------------------+------
z-row      0   -2/3       5/6         0    |  20
x_1        1    2/3       1/6         0    |   4
s_2        0    4/3      -1/6         1    |   2

当前 BFS：x_1=4, x_2=0, s_1=0, s_2=2, z=20。

迭代2：

• 入基：x_2（检验数 -2/3，唯一负值）
• 比值测试：x_1: 4/(2/3)=6, s_2: 2/(4/3)=1.5 → s_2 出基
• 主元：4/3
• 执行转轴操作后：

          x_1   x_2        s_1         s_2  |  RHS
---------------------------------------------+------
z-row      0     0        3/4         1/2  |  21
x_1        1     0        1/4        -1/2  |   3
x_2        0     1       -1/8         3/4  |  1.5

所有检验数 ≥ 0 → 最优。最优解：x_1=3, x_2=1.5, z=21。

3.6 两阶段法（Two-Phase Method）

当初始 BFS 不容易直接获取时（如约束中有 “≥” 或 “=”），需要使用人工变量法。

阶段 I：引入人工变量，构造辅助 LP：

min   Σ a_i  （人工变量之和）
s.t.  Ax + a = b
      x ≥ 0, a ≥ 0

如果最优目标值为 0，说明所有人工变量为 0，原问题有可行解。此时得到原问题的一个 BFS，进入阶段 II。

如果阶段 I 最优目标值 > 0，则原问题无可行解（infeasible）。

阶段 II：使用阶段 I 得出的 BFS 作为初始解，求解原 LP。

3.7 Bland 规则与防止循环

在退化（degenerate）情况下，比值测试可能出现多个相同的最小值，导致算法在选择出基变量时循环（cycling）——即无限地在相同一组基中轮换而不前进。

Bland 规则：

1. 入基：选择下标最小的正（最大化）检验数对应的变量
2. 出基：在比值最小的多个候选中，选择下标最小的基变量

Bland 规则理论上保证了单纯形法的有限终止性。

4. 对偶理论（Duality Theory）

对偶理论是 LP 最优美的部分。每一个 LP（原始问题，Primal）都对应一个对偶问题（Dual），两者之间存在深刻的数学关系。

4.1 对偶问题的构造

原始问题（规范形式，最大化）：

max   c^T x
s.t.  Ax ≤ b
      x ≥ 0

对偶问题：

min   b^T y
s.t.  A^T y ≥ c
      y ≥ 0

重要：

• 原始是最大化 → 对偶是最小化
• 原始约束数 m → 对偶变量数 m（y 的维度）
• 原始变量数 n → 对偶约束数 n
• 约束矩阵转置：A → A^T
• c 和 b 互换角色：原始目标系数 c 变为对偶的右端项，原始右端项 b 变为对偶目标系数

构造法则总结表：

原始（max）	对偶（min）
目标 c^T x	目标 b^T y
约束 i: a_i^T x ≤ b_i	变量 i: y_i ≥ 0
约束 i: a_i^T x ≥ b_i	变量 i: y_i ≤ 0
约束 i: a_i^T x = b_i	变量 i: y_i free
变量 j: x_j ≥ 0	约束 j: A^T_j y ≥ c_j
变量 j: x_j ≤ 0	约束 j: A^T_j y ≤ c_j
变量 j: x_j free	约束 j: A^T_j y = c_j

4.2 弱对偶定理（Weak Duality）

设 x 是原始问题的可行解，y 是对偶问题的可行解。则：

$$c^T x \leq b^T y$$

推论：

• 原始目标函数值的上界由对偶目标值提供
• 如果某个原始解的目标值等于某个对偶解的目标值（c^T x = b^T y），则两者均是最优解

4.3 强对偶定理（Strong Duality）

如果原始问题有最优解 x*，则对偶问题也有最优解 y*，且：

$$c^T x^* = b^T y^*$$

这意味着在最优状态下，原始和对偶的目标值完全相等——对偶间隙（duality gap）为零。

4.4 互补松弛性（Complementary Slackness）

设 x* 和 y* 分别是原始和对偶的最优解，则：

1. 对于原始约束 i：y_i* · (b_i - a_i^T x*) = 0
2. 对于原始变量 j：x_j* · (A_j^T y* - c_j) = 0

直观含义：如果一个原始约束不紧（有松弛，即 a_i^T x* < b_i），则对应的对偶变量必须为零。反之，如果对偶变量 y_i* > 0，则对应的原始约束必为等式（紧约束）。

互补松弛性在经济学中表示影子价格（shadow price）：y_i* 度量了约束 i 的边际价值——增加一单位资源 b_i 带来的目标函数增量。

示例：在上面单纯形法求解的最终表中，检验数行对应松弛变量的系数 (3/4, 1/2) 恰好就是对偶变量的最优值 y* = (0.75, 0.5)。这意味着：

• 资源1（第一个约束的 RHS）每增加 1 单位，目标函数增加 0.75
• 资源2 每增加 1 单位，目标函数增加 0.5

4.5 对偶单纯形法（Dual Simplex）

对偶单纯形法与原始单纯形法对称：

• 原始单纯形：始终保持原始可行性（x_B ≥ 0），逐步满足对偶可行性（检验数 ≥ 0）
• 对偶单纯形：始终保持对偶可行性（检验数 ≥ 0），逐步满足原始可行性（x_B ≥ 0）

应用场景：已有最优解后 RHS b 发生变化（如敏感度分析），或整数规划的分支定界中——此时原始可行性可能被破坏但最优性仍保持，用对偶单纯形法可以更快地重新优化。

5. 内点法（Interior Point Methods）

1984 年，Karmarkar 提出的内点法在理论上具有多项式时间复杂度，在实践中对大规模稀疏 LP 问题通常比单纯形法更快。

5.1 核心思想

与单纯形法在可行域边界顶点移动不同，内点法从可行域内部出发，沿着中心路径（central path） 朝最优解移动，并在接近边界时通过障碍项减速（被”推回”内部）。

5.2 障碍法（Barrier Method）

将 LP 的标准形式转化为带障碍项的等价问题：

$$\text{minimize} \quad c^T x - \mu \sum_{j=1}^{n} \ln(x_j)$$
$$\text{subject to} \quad Ax = b$$

这里 -Σ ln(x_j) 是对数障碍函数（logarithmic barrier），当任何 x_j → 0+ 时，障碍函数 → ∞，从而强制保持 x > 0。

μ > 0 是障碍参数。当 μ → 0 时，障碍问题的最优解收敛到原始 LP 的最优解。

5.3 中心路径

对于每个 μ > 0，障碍问题有唯一最优解 x(μ)，当 μ 从 ∞ 逐渐减小到 0 时，x(μ) 形成的路径就是中心路径。

实际算法不会为每个 μ 精确求解，而是采用路径跟随（path-following）策略：用一个递减序列 μ_k → 0，对每个 μ_k 只做一步（或几步）Newton 迭代来更新解。

5.4 KKT 条件与 Primal-Dual 方法

对于 LP：

min   c^T x
s.t.  Ax = b,  x ≥ 0

其 Lagrange 函数为：
$$L(x, y, s) = c^T x - y^T(Ax - b) - s^T x$$

其中 y 对应等式约束的乘子，s ≥ 0 对应 x ≥ 0 约束的乘子。

KKT 条件（最优性的必要条件，对 LP 也是充分的）：

1. 驻点条件：c - A^T y - s = 0
2. 原始可行性：Ax = b, x ≥ 0
3. 对偶可行性：s ≥ 0
4. 互补松弛：x_j · s_j = 0 (∀j)

Primal-Dual 内点法同时处理原始变量 x 和对偶变量 (y, s)。互补松弛条件 x_j s_j = 0 被松弛为 x_j s_j = μ（中心路径条件），然后逐步让 μ ↓ 0。

5.5 单纯形法 vs 内点法

特性	单纯形法	内点法
复杂度（最坏情况）	指数级（O(2^n)），Klee-Minty 反例	多项式级 O(√n log(1/ε)) 或 O(n^{3.5})
实际运行速度	对小/中型问题极快	对大规模稀疏问题更快
解的性质	精确顶点解（BFS）	近似内部解（需舍入到顶点）
热启动（warm start）	容易	困难
敏感度分析	自然支持（有最优基矩阵信息）	需要额外计算
整数规划中的应用	分支定界中大量使用	较少使用

实践中，商业求解器（CPLEX、Gurobi）同时实现了两种方法，并在运行时自动选择或切换。

5.6 Klee-Minty 反例：单纯形法的最坏情况

Klee 和 Minty 于 1972 年构造了一个 n 维 LP，使用 Dantzig 规则（最大检验数入基）时，单纯形法需要遍历 2^n - 1 个顶点才能到达最优解，证明其最坏情况是指数复杂度。

2 维的 Klee-Minty 问题：

max   2x_1 + x_2
s.t.  x_1       ≤ 5
      4x_1 + x_2 ≤ 25
      x_1, x_2   ≥ 0

对于更高维的推广，单纯形法确实需要指数步数。然而，在绝大多数实际问题中，单纯形法的迭代次数通常仅为约束数的 2~3 倍，这被称为”单纯形法的谜题”——理论和实践的显著差异。Spielman 和 Teng 的”平滑分析”（smoothed analysis，2001）为这一现象提供了理论解释：在轻微随机扰动下，单纯形法的期望复杂度是多项式级的。

6. 整数线性规划（Integer Linear Programming, ILP）

当部分或全部决策变量必须取整数值时，LP 升级为整数规划。ILP 通常是 NP-hard 的（相比之下，LP 是 P 问题）。

6.1 分类

• 纯整数规划（Pure ILP）：所有变量都是整数
• 混合整数规划（Mixed ILP, MILP）：部分变量连续，部分整数
• 0-1 整数规划（Binary ILP）：变量只能取 0 或 1

0-1 规划的强大之处在于可以刻画”是/否”决策：

• 项目 j 是否投资：x_j ∈ {0, 1}
• 任务 i 是否分配给机器 j：x_{ij} ∈ {0, 1}
• 设施 j 是否建设：x_j ∈ {0, 1}

6.2 LP 松弛（LP Relaxation）

求解 ILP 最基本的方法是先忽略整数约束，求解对应的 LP（称为 LP 松弛）。LP 松弛的最优目标值提供 ILP 最优值的界：

• 最大化 ILP：LP 松弛值 ≥ ILP 最优值（上界）
• 最小化 ILP：LP 松弛值 ≤ ILP 最优值（下界）

示例：整数背包问题

max   16x_1 + 19x_2 + 23x_3
s.t.  2x_1 + 3x_2 + 4x_3 ≤ 8
      x_1, x_2, x_3 ∈ {0, 1}

LP 松弛（0 ≤ x_j ≤ 1）的最优解可能是 x = (1, 1, 0.75)，目标值 52.25。但真实的整数最优解是 x = (0, 0, 1) 或 (1, 0, 1) 等等——需要进一步探索。

6.3 分支定界法（Branch and Bound）

分支定界是求解 ILP 最基础的精确算法：

定界（Bound）：求解当前子问题的 LP 松弛。如果 LP 松弛的最优值已经比当前已知的最好整数解差（劣于已找到的解），则剪枝（prune）——该分支不可能含有更好的整数解。

分支（Branch）：在 LP 松弛解中选择一个取值不是整数的变量 x_j，生成两个子问题：

1. 加入约束 x_j ≤ ⌊x_j⌋
2. 加入约束 x_j ≥ ⌈x_j⌉

然后在子问题上递归地求解。

搜索策略：

• 深度优先（DFS）：快速找到一个可行整数解
• 最佳界优先（Best-First）：优先探索最优 LP 松弛值的节点（通常更快找到最优解）

# 分支定界伪代码
def branch_and_bound(problem, best_obj=float('-inf'), best_sol=None):
    # 求解 LP 松弛
    lp_sol, lp_obj = solve_lp(problem)
    
    # 不可行 → 剪枝
    if not lp_sol:
        return best_sol, best_obj
    
    # LP 目标值 ≤ 当前最好解 → 剪枝（最大化问题）
    if lp_obj <= best_obj:
        return best_sol, best_obj
    
    # 检查是否全整数
    if all(is_integer(x_j) for x_j in lp_sol):
        if lp_obj > best_obj:
            return lp_sol, lp_obj
        else:
            return best_sol, best_obj
    
    # 分支：选择一个非整数变量
    j = choose_branching_var(lp_sol)
    floor_val = int(np.floor(lp_sol[j]))
    
    # 左分支：x_j ≤ floor_val
    left_sol, left_obj = branch_and_bound(
        problem + constraint(x_j <= floor_val),
        best_obj, best_sol
    )
    if left_obj > best_obj:
        best_obj, best_sol = left_obj, left_sol
    
    # 右分支：x_j ≥ floor_val + 1
    right_sol, right_obj = branch_and_bound(
        problem + constraint(x_j >= floor_val + 1),
        best_obj, best_sol
    )
    if right_obj > best_obj:
        best_obj, best_sol = right_obj, right_sol
    
    return best_sol, best_obj

6.4 割平面法（Cutting Plane）

割平面法的思想是：求解 LP 松弛，如果解不是整数，则添加一个线性约束（割平面）来割掉当前的非整数 LP 解，但不割掉任何可行的整数解。反复进行直到得到整数解。

Gomory 割：经典方法，从单纯形表的一行生成：

若 x_{B_i} + Σ f_{ij} x_{N_j} = f_{i0}
其中 f_{i0} 不是整数，则添加约束：
x_{B_i} + Σ ⌊f_{ij}⌋ x_{N_j} ≤ ⌊f_{i0}⌋

Gomory 证明了这种割的有限收敛性。实践中，割平面通常与分支定界结合——称为分支-切割法（Branch and Cut）。

6.5 求解 ILP 的常用策略总结

方法	精确/近似	适用规模	特点
穷举	精确	n < 20	仅限极小问题
分支定界 (B&B)	精确	n < 200	通用框架，效率依赖 LP 松弛质量
割平面	精确	n < 200	通常与 B&B 结合
分支-切割 (B&C)	精确	n < 500	现代 MILP 求解器标准方法
LP 松弛 + 舍入	近似	n > 1000	快速但可能较差
启发式/元启发式	近似	n > 1000	遗传算法、禁忌搜索等

7. 敏感度分析（Sensitivity Analysis）

敏感度分析回答的问题是：当模型参数发生变化时，最优解如何改变？这在现实中至关重要——数据往往是不精确的，我们需要知道最优决策对参数变化的稳健性。

7.1 目标系数的变化范围

对于基变量 x_{B_j} 的目标系数 c_{B_j}，其允许的变化范围是：

$$\max \left\{ \frac{-\sigma_k}{\bar{a}_{jk}} : \bar{a}_{jk} > 0 \right\} \leq \Delta c_{B_j} \leq \min \left\{ \frac{-\sigma_k}{\bar{a}_{jk}} : \bar{a}_{jk} < 0 \right\}$$

其中 σ_k 是当前表中非基变量的检验数，\bar{a}_{jk} 是转轴后的约束矩阵元素。

如果 Δc 在此范围内，最优基不变（检验数符号不变）；但如果 Δc 超出这个范围，需要使用对偶单纯形法重新优化。

7.2 右端项（RHS）的变化与影子价格

对于约束 i 的 RHS b_i 的变化：

$$\max \left\{ \frac{-\bar{b}_r}{\bar{B}^{-1}_{ri}} : \bar{B}^{-1}_{ri} > 0 \right\} \leq \Delta b_i \leq \min \left\{ \frac{-\bar{b}_r}{\bar{B}^{-1}_{ri}} : \bar{B}^{-1}_{ri} < 0 \right\}$$

在此范围内，最优基不变。b_i 每增加 1 单位，目标函数增加 y_i*（影子价格）。

影子价格的经济解释：如果约束 i 代表”原材料 A 的可用量不超过 b_i 吨”，那么 y_i*（影子价格）就是额外 1 吨原材料的边际价值——你愿意为它支付多少。如果市场价格低于影子价格，你应该购买更多；反之则不应。

7.3 100% 规则

当多个系数同时变化时，可以使用 100% 规则做近似检查：

• 对每个变化计算 ratio = (实际变化量) / (允许的最大变化量)
• 如果 Σ ratio ≤ 1，最优基不变
• 如果 Σ ratio > 1，无法确定，需要实际重新求解

8. Python 实现：LP 求解实操

8.1 scipy.optimize.linprog

SciPy 提供了 linprog 函数，支持多种求解算法：

from scipy.optimize import linprog
import numpy as np

# 问题：min c^T x  s.t. A_ub x ≤ b_ub, A_eq x = b_eq, bounds
# 示例：min -2x1 - 3x2  s.t. x1 + 2x2 ≤ 8, 4x1 ≤ 16, x1,x2 ≥ 0

c = [-2, -3]  # 注意是 -2, -3 因为我们最小化（等价于 max 2x1+3x2）

A_ub = [[1, 2],
        [4, 0],
        [0, 4]]
b_ub = [8, 16, 12]

bounds = [(0, None), (0, None)]  # x1 ≥ 0, x2 ≥ 0

# 使用 HiGHS 求解器（SciPy 1.6+）
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

print(f"状态: {result.message}")
print(f"最优值: {-result.fun:.2f}")  # 注意符号翻转
print(f"最优解: x1 = {result.x[0]:.2f}, x2 = {result.x[1]:.2f}")
print(f"迭代次数: {result.nit}")

SciPy 支持的方法：

• method='highs'（推荐）：HiGHS 实现，性能优秀
• method='simplex'（已弃用）：经典两阶段单纯形法
• method='interior-point'（已弃用）：Primal-Dual 内点法

8.2 PuLP（教学首选）

PuLP 提供优雅的建模语法，适合教学和中小规模问题：

from pulp import *

# 创建问题
prob = LpProblem("Production_Planning", LpMaximize)

# 变量
x1 = LpVariable("x1", lowBound=0)   # 产品1产量
x2 = LpVariable("x2", lowBound=0)   # 产品2产量

# 目标函数
prob += 2*x1 + 3*x2

# 约束
prob += x1 + 2*x2 <= 8
prob += 4*x1 <= 16
prob += 4*x2 <= 12

# 求解
prob.solve()

print(f"状态: {LpStatus[prob.status]}")
print(f"最优值: {value(prob.objective)}")
print(f"x1 = {value(x1)}, x2 = {value(x2)}")

# 影子价格（对偶变量）
for name, constraint in prob.constraints.items():
    print(f"{name}: 影子价格 = {constraint.pi}, 松弛 = {constraint.slack}")

8.3 CVXPY（凸优化框架）

CVXPY 定位为领域特定语言（DSL），可求解 LP、QP、SDP 等多种凸优化问题：

import cvxpy as cp

# 变量
x = cp.Variable(2, nonneg=True)

# 目标
objective = cp.Maximize(2*x[0] + 3*x[1])

# 约束
constraints = [
    x[0] + 2*x[1] <= 8,
    4*x[0] <= 16,
    4*x[1] <= 12
]

# 求解
prob = cp.Problem(objective, constraints)
prob.solve()

print(f"最优值: {prob.value:.2f}")
print(f"最优解: {x.value}")

CVXPY 的优势：如果问题是凸的（包括 LP），CVXPY 可以自动验证凸性并使用最合适的求解器。代码几乎就是数学公式的直接翻译。

8.4 Gurobi / CPLEX（商业级求解器）

对于工业级 MILP 问题，Gurobi 和 CPLEX 是行业标准：

# Gurobi 示例
import gurobipy as gp
from gurobipy import GRB

# 创建模型
m = gp.Model("production")

# 变量
x1 = m.addVar(name="x1")
x2 = m.addVar(name="x2")

# 目标
m.setObjective(2*x1 + 3*x2, GRB.MAXIMIZE)

# 约束
m.addConstr(x1 + 2*x2 <= 8)
m.addConstr(4*x1 <= 16)
m.addConstr(4*x2 <= 12)

# 求解
m.optimize()

print(f"最优值: {m.objVal}")
for v in m.getVars():
    print(f"{v.varName} = {v.x}")

Gurobi/CPLEX 的优势：

• 对 MILP 的启发式算法极强（预求解、割生成、强分支策略）
• 支持并行计算
• 学术/免费许可可用

8.5 工具对比

工具	许可证	适用规模	特点
scipy.optimize.linprog	BSD	中小型 LP	轻量，NumPy 生态，HiGHS 后端
PuLP	MIT	中小型 LP/MILP	语法友好，可调用多种求解器
CVXPY	Apache 2.0	中小型 LP/QP/SDP	DSL 风格，自动凸性验证
OR-Tools (Google)	Apache 2.0	中大型 MILP	Google 出品，约束编程强
HiGHS	MIT	大型 LP	目前最快的开源 LP 求解器
GLPK	GPL	中型 LP/MILP	GNU 项目，教学使用
Gurobi	商业/Academic	工业级	顶级 MILP 性能
CPLEX (IBM)	商业/Academic	工业级	与 Gurobi 齐名

推荐选择：

• 学习/教学：PuLP + GLPK
• 纯 LP 研究/原型：scipy.linprog (HiGHS)
• 生产环境 MILP：Gurobi 或 CPLEX
• 凸优化研究：CVXPY

9. 经典应用场景

9.1 资源分配问题

场景：工厂生产 n 种产品，有 m 种资源约束（原料、工时、机器容量）。每种产品 j 需要 a_{ij} 单位资源 i，利润为 c_j。求最大利润的生产方案。

LP 模型：

max   Σ c_j x_j
s.t.  Σ a_{ij} x_j ≤ b_i  (∀i)
      x_j ≥ 0

9.2 运输问题（Transportation Problem）

场景：有 p 个供应地（工厂），q 个需求地（零售店）。工厂 i 产能为 s_i，零售店 j 需求为 d_j，从 i 到 j 的单位运输成本为 c_{ij}。如何安排运输使总成本最小？

LP 模型：

min   Σ_i Σ_j c_{ij} x_{ij}
s.t.  Σ_j x_{ij} = s_i   (∀i)   — 每个工厂运出量等于产能
      Σ_i x_{ij} = d_j   (∀j)   — 每个零售店收到量等于需求
      x_{ij} ≥ 0

运输问题是一个高度结构化的问题：约束矩阵是全幺模（totally unimodular）的，因此 LP 松弛自动给出整数解（当 s_i 和 d_j 为整数时）。这意味着运输问题可以在多项式时间内精确求解。

9.3 指派问题（Assignment Problem）

场景：n 个工人，n 个任务。工人 i 完成任务 j 的成本为 c_{ij}。每人分配恰好一个任务，每个任务恰好分配一人。

LP 模型：

min   Σ_i Σ_j c_{ij} x_{ij}
s.t.  Σ_j x_{ij} = 1   (∀i)   — 每人一个任务
      Σ_i x_{ij} = 1   (∀j)   — 每个任务一人
      x_{ij} ≥ 0  （或 x_{ij} ∈ {0,1}）

由于约束矩阵的全幺模性，LP 松弛自动产生 0-1 解。匈牙利算法可以在 O(n^3) 时间内求解，比单纯形法更高效。

9.4 网络流问题（Network Flow）

最大流问题：有向图 G=(V,E)，源点 s，汇点 t，每条边 (u,v) 有容量 c_{uv}。求从 s 到 t 的最大流量。

LP 模型：

max   Σ_v f_{sv}   (源点总流出)
s.t.  0 ≤ f_{uv} ≤ c_{uv}        (∀(u,v)∈E)  — 容量约束
      Σ_{v:(u,v)∈E} f_{uv} - Σ_{v:(v,u)∈E} f_{vu} = 0  (∀u∈V\{s,t})  — 流量守恒

虽然可以写成 LP，但实际中从不这么求解——Ford-Fulkerson 算法（O(E·max_flow)）或 Dinic 算法（O(V^2·E)）高效得多。这说明了解问题结构的重要性：LP 是通用工具，但专用算法几乎总是更快。

9.5 投资组合优化

场景（简化版）：n 个投资标的，每个有预期收益率 r_j。总预算为 B。需要在给定风险约束下最大化预期收益。

LP 模型（线性版本，使用绝对偏差代替方差）：

max   Σ r_j x_j
s.t.  Σ x_j = B             — 预算约束
      Σ |r_j - r̄| x_j ≤ R  — 风险约束（线性化 MAD）
      x_j ≥ 0

更精确的模型是二次规划（QP）——Markowitz 均值-方差模型。

9.6 生产计划与库存管理

场景：T 个时间周期，每个周期 t 有需求 d_t、生产成本 c_t、库存持有成本 h_t。决策变量：产量 p_t 和库存量 inv_t。

LP 模型：

min   Σ_t (c_t · p_t + h_t · inv_t)
s.t.  inv_{t-1} + p_t = d_t + inv_t   (∀t)  — 库存平衡
      p_t ≥ 0, inv_t ≥ 0

10. 多目标线性规划（Multi-Objective LP）

现实中的决策问题往往涉及多个相互冲突的目标。例如”利润最大化”与”污染最小化”通常不可兼得。

10.1 目标规划（Goal Programming）

目标规划将每个目标设定一个期望值（目标水平），然后最小化与目标的偏差：

形式：

min   Σ w_k (d_k^+ + d_k^-)
s.t.  c_k^T x - d_k^+ + d_k^- = g_k   (∀k)  — 目标约束
      Ax ≤ b                              — 硬约束
      d_k^+, d_k^- ≥ 0

其中：

• g_k 是目标 k 的期望水平
• d_k^+ 是超出的偏差（≥ 0）
• d_k^- 是不足的偏差（≥ 0）
• w_k 是该目标的权重

扩展：

• 优先级（Preemptive）目标规划：目标有严格优先级 P_1 >> P_2 >> …，必须先满足高优先级的再考虑低优先级的
• 加权目标规划：所有目标同时考虑，通过权重 w 反映相对重要性

10.2 Pareto 前沿

在多目标优化中，一个解称为Pareto 最优，如果不存在另一个解在所有目标上都不差于此解且至少在一个目标上严格优于它。

获取 Pareto 前沿的常用方法之一是ε-约束法：将 k-1 个目标转化为约束，优化剩余一个目标，通过改变约束阈值遍历 Pareto 前沿。

min   f_1(x)
s.t.  f_2(x) ≤ ε_2
      f_3(x) ≤ ε_3
      ...
      x ∈ X

变化 ε_2, ε_3, … 即可获得 Pareto 前沿的不同点。

11. 线性规划与其他优化问题的比较

11.1 LP vs QP（二次规划）

特性	LP	QP
目标函数	线性：c^T x	二次：1/2 x^T Q x + c^T x
约束	线性	线性
凸性	自动凸（线性即凸）	需要 Q 半正定才凸
最优解位置	总是在顶点	可能在约束的内部或边界上
复杂度	P 问题	凸时 P 问题
经典应用	资源分配、运输	投资组合、SVM、MPC

11.2 LP vs SDP（半定规划）

特性	LP	SDP
决策变量	向量 x ∈ R^n	对称矩阵 X ∈ S^n
非负约束	x ≥ 0	X ≽ 0（半正定性）
应用	运筹优化	组合优化松弛、控制系统、量子信息
复杂度	P 问题	P 问题（但更”昂贵”，≈ O(n^6)）

SDP 可以视为 LP 的推广——将”元素非负”推广为”矩阵半正定”。许多 NP-hard 组合优化问题（如 MAX-CUT）有高质量的 SDP 松弛。

11.3 什么时候 LP 够用？（何时不需要更复杂的模型）

• LP 适用：问题天然具有比例性（线性）+ 可加性
• 需要 IP：存在”是/否”决策、数量必须为整数
• 需要 QP：目标函数存在二次项（如风险最小化、力最小化）
• 需要 NLP：存在非线性物理/化学规律（如化学平衡、电力潮流）
• 需要 SDP：需要矩阵变量或组合优化的紧松弛

12. 计算复杂度与实用性能

12.1 LP 的理论复杂度

• 椭球法（Khachiyan，1979）：首次证明 LP ∈ P，复杂度 O(n^6 L^2)。理论上重要，但实践上极慢。
• 内点法（Karmarkar，1984）：O(n^{3.5} L) 级，既有理论保证又有实用性能。
• 单纯形法：最坏情况指数级，但平均/实用时为 O(m log n) 到 O(n) 次迭代。

12.2 为什么单纯形法在大多数实际问题上这么快？

这个问题困扰了优化社区数十年。Spielman 和 Teng（2001）的平滑分析给出了部分答案：通过对输入系数施加小的随机扰动，单纯形法的期望运行时间是多项式级的。这意味着”病态”的 Klee-Minty 型问题在实际数据中极为罕见——实际数据本身就有天然的”噪声”。

直观地说，单纯形法的每次迭代大约需要 O(mn) 的工作量（对单纯形表做 Gauss-Jordan 消元），而迭代次数通常在 O(m) 到 O(2m) 之间。因此总复杂度约为 O(m^2 n)，在 m << n 时非常高效。

12.3 选择求解策略的实用决策树

问题规模？
├── n < 1000, m < 1000 → 单纯形法（简单、有顶点解）
│   └── 需要敏感度分析 → 单纯形法（有基矩阵信息）
├── n > 10000, m > 10000 → 内点法或 HiGHS
│   └── 稀疏约束矩阵 → 内点法（利用稀疏线性代数）
└── 含整数变量？
    ├── n < 200 → B&B / Branch-and-Cut
    ├── n > 500 → 考虑启发式（时间有限的场景）
    └── n > 1000 → LP 松弛 + 舍入，或使用商业 MILP 求解器

13. 数值稳定性与预求解（Presolve）

13.1 数值问题

LP 求解中常见的数值问题：

• 缩放不佳：约束系数量级差异巨大（如 10^{-6} 和 10^6 同时出现），导致舍入误差
• 退化：多个约束交汇于同一点，导致基矩阵接近奇异
• 接近不可行：约束几乎矛盾

解决：商业求解器自动进行预求解（Presolve），包括缩放（scaling）、冗余约束移除、变量固定（固定能推出的变量值）等。

13.2 预求解示例

# 原始问题
min   x1 + 2x2
s.t.  x1 + x2 ≥ 100
      x1 + x2 ≤ 100
      1000x1 + 1000x2 = 100000  # 冗余！
      x2 ≥ 50
      x1 ≥ 0

# 预求解后：
# 约束1 + 约束2 → x1 + x2 = 100（等式化）
# 约束3 除以 1000 → x1 + x2 = 100（与上面相同，删除）
# x1 + x2 = 100, x2 ≥ 50 → x1 ≤ 50
# 消元 x1 = 100 - x2
# 问题简化为：min (100 - x2) + 2x2 = 100 + x2, s.t. 50 ≤ x2 ≤ 100
# → x2 = 50 最优

现代求解器的预求解步骤可能将问题规模缩减 50% 以上，某些极端情况甚至直接求解出最优解。

14. 面试高频问答

Q1: 请解释单纯形法为什么总是在顶点之间移动，而不是在可行域内部搜索？内点法与单纯形法有什么根本区别？

A: 单纯形法依赖 LP 的几何性质——LP 的可行域是凸多面体，目标函数是线性的。线性函数在凸多面体上的极值一定在所有顶点上取得（如果最优解存在）。用反证法：如果一个内部点 x 是最优解，那么存在一个可移动的方向使目标函数值不变（否则就不是最优点），顺着该方向走到可行域边界，不断重复最终到达顶点，目标值不会变差。因此只需检查顶点即可找到最优解。

单纯形法沿着边界从顶点走到相邻顶点，本质是一种组合搜索。内点法则从内部沿”中心路径”逼近最优解——它使用 Newton 步在约束定义的流形上移动。两者的根本区别：

• 单纯形法：沿可行域边界的离散跳跃，每次迭代走一条边
• 内点法：在可行域内部的连续路径上移动，避开边界直到最后一刻

从复杂度看，单纯形法最坏情况是指数级（Klee-Minty），但通常只要 O(m) 次迭代；内点法理论上是多项式级，但每次迭代求解 Newton 系统的代价更高。实践中，单纯形法对中小规模或需要热启动的场景更好，内点法对大规模稀疏问题更优。

---

Q2: 什么是影子价格？如何从单纯形法的最终表中获取影子价格？它在管理决策中如何使用？

A: 影子价格（shadow price）是约束的对偶变量最优值 y_i*，表示约束 i 右端项 b_i 增加一单位时目标函数最优值的变化率。

在单纯形法的最终表中，影子价格就是目标函数行（检验数行）中初始基变量（如松弛变量）对应的系数。以上面的手工计算例子：最终表中松弛变量 s_1 的检验数为 3/4，s_2 为 1/2，说明：

• 资源 1 每增加 1 单位，最优利润增加 0.75
• 资源 2 每增加 1 单位，最优利润增加 0.5

在管理决策中：

• 如果某种资源的市场价格 < 影子价格 → 购买更多是有利可图的
• 如果某种资源的市场价格 > 影子价格 → 甚至可以卖出部分资源获利
• 影子价格为 0 的约束 → 该资源有剩余（非紧约束），增加它不会改善目标（但也不会恶化）

需要注意的是，影子价格只在 b_i 的变化范围内有效。超过这个范围，最优基改变，影子价格也会改变。这体现了敏感度分析的重要性——不能无限外推影子价格的适用性。

---

Q3: 请解释整数规划为什么比线性规划难得多？为什么不能简单地求解 LP 松弛然后四舍五入？

A: LP 是 P 问题（多项式时间可解），而 ILP 是 NP-hard 问题。根本原因在于整数约束破坏了凸性——可行域从凸多面体变成一组孤立的整数点。

简单地四舍五入有多个严重问题：

1. 可能违反约束：LP 松弛解 x_1=2.3, x_2=1.7 舍入为 (2, 2)，但在给定约束下 (2, 2) 可能不可行。

2. 可能远离最优：以下面的一维例子说明：

   max  x
   s.t. x ≤ 100.4
        x ≥ 0, x 整数

   LP 松弛给出 x=100.4。四舍五入得 x=100 是好的。但考虑：

   max  10x
   s.t. x ≤ 0.49
        x ≥ 0, x 整数

   LP 松弛得 x=0.49，四舍五入为 0。但真正的整数最优解就是 0。虽然这个例子刚好正确，但换个系数：

   max  c · x
   s.t. a · x ≤ b, x ∈ {0, 1}

   令 c/a 很大时，LP 解可能是 x=1，但实际约束可能在整数点下精确阻止了 x=1。

3. 整数约束改变了可行域的形状：整数可行域不是凸的，LP 松弛的最优解和 ILP 的最优解可能相距甚远，无论你怎么舍入都无法到达。

更具体的例子（0-1 背包）：

max  15x_1 + 12x_2 + 10x_3
s.t. 8x_1 + 6x_2 + 5x_3 ≤ 10
     x_j ∈ {0, 1}

LP 松弛（0 ≤ x_j ≤ 1）：x = (1, 1/3, 0), 目标 = 15 + 4 = 19。四舍五入得 (1, 0, 0) 目标 = 15。但真正的整数最优解是 (0, 0, 1) 目标 = 10？不对… 让我重新算。x_1=1 占 8 单位重量，剩余 2 单位，放不下 x_2 (6) 或 x_3 (5)。所以只有 x_1=1，目标 15。LP 松弛的 (1, 0.33, 0) 四舍五入为 (1, 0, 0)，恰好是这个例子中最优的 15。不过设计一个使舍入失败的例子也很容易，此处不赘述。关键是：ILP 必须使用分支定界、割平面等专用方法，而非简单舍入。

---

Q4: 如何判断一个 LP 是否有可行解？如果不可行，怎么找到”最接近”可行的解？

A: 判断 可行性有两种方法：

方法一：两阶段法中的阶段 I。引入人工变量 a_i ≥ 0，求解：

min  Σ a_i
s.t. Ax + a = b
     x ≥ 0, a ≥ 0

如果此问题的最优值为 0，则存在可行解（a_i = 0, ∀i）。如果 > 0，则原 LP 不可行。

方法二：不可行性证明（Farkas 引理）。Farkas 引理指出，LP 不可行当且仅当存在向量 y 使得：

A^T y ≥ 0,  b^T y < 0

找到了这样的 y 即提供了一个不可行性证明（certificate of infeasibility）。商业求解器在报告不可行时通常会提供这样的证明。

找到”最接近”可行的解：如果问题不可行，可以通过以下方式松弛：

1. 引入人工变量到每个约束，允许约束被”弹性化”（elastic programming）
2. 求解 min Σ（违反量），得到最小违反的解
3. 或者使用目标规划将硬约束转化为带惩罚项的软约束：

   min  c^T x + M · Σ s_i
   s.t. Ax - s ≤ b  （或类似形式）
        s_i ≥ 0  （s_i 为约束违反量）

   其中 M 是一个大的惩罚系数。这种方法给出了”在约束几乎满足的情况下，最优决策是什么”的答案。

实际工程中，当模型报告不可行时，通常需要回到问题本身，检查是否有矛盾的业务约束——如”预算 ≤ 100 万”和”最低支出 ≥ 200 万”同时存在。这是建模错误，而非求解问题。

---

Q5: 内点法中的 barrier method 是什么？为什么参数 μ 要从大逐步减小？这和模拟退火中的温度参数有何相似之处？

A: Barrier method 的基本思想是将带不等式约束 x ≥ 0 的 LP 转化为只含等式约束但目标函数带障碍项的问题：

$$\min \quad c^T x - \mu \sum \ln(x_j) \quad \text{s.t.} \quad Ax = b$$

对数障碍函数 -ln(x_j) 在 x_j → 0^+ 时 → ∞，这确保了解始终保持在内部（x_j > 0），不会”碰到”边界。

μ 的角色：μ 控制障碍项的强度。

• μ 很大 → 障碍项主导，解被强力推向可行域中心（远离边界），此时优化的是障碍函数，而非原始目标 c^T x
• μ → 0 → 障碍项影响减弱，解逐渐接近原始 LP 的最优解（通常在边界上）

为什么逐步减小：如果一开始就用非常小的 μ，障碍函数曲面会非常”陡峭”——在边界附近急剧变化，使 Newton 迭代很难收敛。从大的 μ 开始，曲面平缓，容易找到解；然后逐步减小 μ，每次用上一步的解作为热启动，继续追踪中心路径，直到逼近真实最优解。

与模拟退火的类比：

• 退火中的温度 T：高温允许在能量面上大幅跳跃（探索），低温则精确逼近局部极小
• 内点法中的 μ：大 μ 使解保持在内部”安全区”（探索），小 μ 则精确逼近边界上的最优解

但有一个关键区别：退火通常用于非凸优化，需要在探索（exploration）和利用（exploitation）之间权衡；而 LP 是凸的，中心路径是唯一的，只要遵循中心路径就会必然收敛到全局最优。因此内点法的 μ 减小策略不像退火的降温策略那样需要精心设计——简单的几何递减 μ_{k+1} = σμ_k（σ ∈ [0.01, 0.5]）通常就足够好了。