线性规划

线性规划是运筹学中最基础、应用最广泛的优化方法之一。它研究在一组线性约束条件下，如何使线性目标函数达到最优值的问题。从资源分配到生产调度，从物流运输到投资组合，线性规划为决策者提供了强有力的数学工具。

基本原理

问题定义

线性规划（Linear Programming, LP）是在满足一组线性等式或不等式约束的条件下，求解线性目标函数最大值或最小值的数学方法。

一个线性规划问题包含三个基本要素：

1. 决策变量：需要确定的未知量，通常记为 \( x_1, x_2, \ldots, x_n \)
2. 目标函数：需要最大化或最小化的线性函数
3. 约束条件：决策变量必须满足的线性等式或不等式

数学表述

线性规划问题的一般形式为：

\[ \min \quad z = c_1 x_1 + c_2 x_2 + \cdots + c_n x_n \]

约束条件：

\[ \begin{cases} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n \leq b_1 \\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n \leq b_2 \\ \vdots \\ a_{m1}x_1 + a_{m2}x_2 + \cdots + a_{mn}x_n \leq b_m \\ x_1, x_2, \ldots, x_n \geq 0 \end{cases} \]

基本定理

线性规划的理论基础建立在以下几个关键定理之上：

• 可行域凸性定理：线性规划的可行域是凸多面体（可能无界或为空集）
• 最优解存在性定理：若可行域非空且有界，则最优解一定存在
• 顶点最优性定理：若线性规划有最优解，则一定存在一个顶点（基本可行解）是最优解

这些定理为单纯形法提供了理论基础——我们只需在可行域的顶点中搜索最优解。

标准形式与转化

标准形式

线性规划的标准形式为：

\[ \min \quad \mathbf{c}^T \mathbf{x} \]

\[ \text{s.t.} \quad A\mathbf{x} = \mathbf{b}, \quad \mathbf{x} \geq \mathbf{0} \]

其中 \( \mathbf{c} \in \mathbb{R}^n \)，\( A \in \mathbb{R}^{m \times n} \)，\( \mathbf{b} \in \mathbb{R}^m \)，\( \mathbf{b} \geq \mathbf{0} \)。

转化规则

任何线性规划问题都可以转化为标准形式，转化规则如下：

目标函数转化：最大化问题转为最小化：

\[ \max ; z = \mathbf{c}^T \mathbf{x} \quad \Longleftrightarrow \quad \min ; (-z) = -\mathbf{c}^T \mathbf{x} \]

不等式约束转化：

• “\( \leq \)” 约束加入松弛变量 \( s_i \geq 0 \)：

\[ a_{i1}x_1 + \cdots + a_{in}x_n + s_i = b_i \]

• “\( \geq \)” 约束减去剩余变量 \( s_i \geq 0 \)：

\[ a_{i1}x_1 + \cdots + a_{in}x_n - s_i = b_i \]

无约束变量转化：若 \( x_j \) 无非负约束，令 \( x_j = x_j^+ - x_j^- \)，其中 \( x_j^+, x_j^- \geq 0 \)。

转化示例

将以下问题转化为标准形式：

\[ \max \quad z = 2x_1 + 3x_2 \]

\[ \text{s.t.} \quad x_1 + 2x_2 \leq 8, \quad 4x_1 + 3x_2 \leq 12, \quad x_1, x_2 \geq 0 \]

转化后：

\[ \min \quad z’ = -2x_1 - 3x_2 \]

\[ \text{s.t.} \quad x_1 + 2x_2 + s_1 = 8, \quad 4x_1 + 3x_2 + s_2 = 12, \quad x_1, x_2, s_1, s_2 \geq 0 \]

图解法（二维）

方法原理

当线性规划问题只有两个决策变量时，可以用图解法求解。基本步骤为：

1. 在二维坐标系中画出每个约束条件对应的直线
2. 确定每个约束对应的半平面，取交集得到可行域
3. 画出目标函数的等值线
4. 沿目标函数值增大（或减小）的方向平移等值线
5. 等值线最后接触可行域的点即为最优解

求解示例

考虑如下问题：

\[ \max \quad z = 3x_1 + 5x_2 \]

\[ \text{s.t.} \quad x_1 \leq 4, \quad 2x_2 \leq 12, \quad 3x_1 + 2x_2 \leq 18, \quad x_1, x_2 \geq 0 \]

步骤一：画出约束边界线

• \( x_1 = 4 \)（竖直线）
• \( x_2 = 6 \)（水平线）
• \( 3x_1 + 2x_2 = 18 \)（斜线，过点 \((6,0)\) 和 \((0,9)\)）

步骤二：确定可行域

可行域为满足所有约束的交集区域，其顶点为：

• \( O = (0, 0) \)
• \( A = (4, 0) \)
• \( B = (4, 3) \)（由 \( x_1 = 4 \) 与 \( 3x_1 + 2x_2 = 18 \) 的交点）
• \( C = (2, 6) \)（由 \( 2x_2 = 12 \) 与 \( 3x_1 + 2x_2 = 18 \) 的交点）
• \( D = (0, 6) \)

步骤三：计算各顶点目标函数值

• \( z(O) = 3(0) + 5(0) = 0 \)
• \( z(A) = 3(4) + 5(0) = 12 \)
• \( z(B) = 3(4) + 5(3) = 27 \)
• \( z(C) = 3(2) + 5(6) = 36 \)
• \( z(D) = 3(0) + 5(6) = 30 \)

结论：最优解为 \( x_1 = 2, x_2 = 6 \)，最优值 \( z^* = 36 \)。

特殊情况

图解法还能直观展示线性规划的几种特殊情况：

• 无可行解：各约束的半平面交集为空
• 无界解：可行域在目标函数优化方向上无界
• 多重最优解：目标函数等值线与可行域的一条边重合
• 退化解：某顶点由多于必要数量的约束边界线交汇而成

单纯形法原理与步骤

基本思想

单纯形法是 George Dantzig 于 1947 年提出的求解线性规划的经典算法。其核心思想是：

1. 从一个基本可行解（顶点）出发
2. 沿可行域的边移动到相邻顶点，使目标函数值改善
3. 重复步骤 2 直到无法继续改善，即达到最优解

基本概念

基变量与非基变量：设标准形式中 \( A \) 为 \( m \times n \) 矩阵且 \( \text{rank}(A) = m \)。选取 \( m \) 个线性无关的列组成基矩阵 \( B \)，对应的变量为基变量 \( \mathbf{x}_B \)，其余为非基变量 \( \mathbf{x}_N \)。

基本可行解：令 \( \mathbf{x}_N = \mathbf{0} \)，由 \( B\mathbf{x}_B = \mathbf{b} \) 解得 \( \mathbf{x}_B = B^{-1}\mathbf{b} \)。若 \( \mathbf{x}_B \geq \mathbf{0} \)，则为基本可行解。

单纯形表

单纯形法通常通过单纯形表来组织计算。表的结构如下：

算法步骤

步骤一：初始化

构造初始基本可行解。若标准形式中已有单位矩阵列（如松弛变量），直接作为初始基。

步骤二：最优性检验

计算各非基变量的检验数（reduced cost）：

\[ \sigma_j = c_j - \mathbf{c}_B^T B^{-1} \mathbf{a}_j \]

对于最小化问题，若所有 \( \sigma_j \geq 0 \)，当前解为最优解。

步骤三：确定进基变量

选取检验数最小的非基变量 \( x_k \)（\( \sigma_k < 0 \) 且 \( |\sigma_k| \) 最大）作为进基变量。

步骤四：确定出基变量

计算最小比值（theta 准则）：

\[ \theta = \min_{i: \bar{a}_{ik} > 0} \left{ \frac{\bar{b}i}{\bar{a}{ik}} \right} \]

对应行的基变量为出基变量。若所有 \( \bar{a}_{ik} \leq 0 \)，则问题无界。

步骤五：基变换（旋转操作）

以 \( \bar{a}_{rk} \)（主元素）为中心，进行初等行变换，将主元列化为单位向量。返回步骤二。

计算实例

求解：\( \min ; z = -2x_1 - 3x_2 \)，约束 \( x_1 + 2x_2 + s_1 = 8 \)，\( 4x_1 + 3x_2 + s_2 = 12 \)，所有变量非负。

初始表（基变量为 \( s_1, s_2 \)）：

检验数 \( \sigma_2 = -3 \) 最小，\( x_2 \) 进基。比值：\( 8/2 = 4 \)，\( 12/3 = 4 \)，取第一行，\( s_1 \) 出基。

第一次迭代后：

主元素为第1行第2列的 2，进行行变换：

检验数 \( \sigma_1 = -1/2 < 0 \)，\( x_1 \) 进基。比值：\( 4/(1/2) = 8 \)，\( 0/(5/2) = 0 \)，取第二行，\( s_2 \) 出基。

第二次迭代后：

所有检验数非负，达到最优。最优解 \( x_1 = 0, x_2 = 4 \)，\( z_{\min} = -12 \)（即原最大化问题 \( z_{\max} = 12 \)）。

对偶问题

对偶理论

每个线性规划问题（原问题）都有一个与之对应的对偶问题。对偶关系揭示了线性规划的深层结构。

原问题（P）：

\[ \min \quad \mathbf{c}^T \mathbf{x}, \quad \text{s.t.} \quad A\mathbf{x} \geq \mathbf{b}, ; \mathbf{x} \geq \mathbf{0} \]

对偶问题（D）：

\[ \max \quad \mathbf{b}^T \mathbf{y}, \quad \text{s.t.} \quad A^T\mathbf{y} \leq \mathbf{c}, ; \mathbf{y} \geq \mathbf{0} \]

对偶基本定理

• 弱对偶定理：若 \( \mathbf{x} \) 是原问题的可行解，\( \mathbf{y} \) 是对偶问题的可行解，则 \( \mathbf{c}^T\mathbf{x} \geq \mathbf{b}^T\mathbf{y} \)
• 强对偶定理：若原问题和对偶问题都有可行解，则两者的最优值相等
• 互补松弛定理：设 \( \mathbf{x}^* \) 和 \( \mathbf{y}^* \) 分别是原问题和对偶问题的最优解，则：

\[ y_i^* (A_i \mathbf{x}^* - b_i) = 0, \quad i = 1, \ldots, m \]

\[ x_j^* (c_j - \mathbf{y}^{*T} \mathbf{a}_j) = 0, \quad j = 1, \ldots, n \]

对偶问题的经济解释

对偶变量 \( y_i \) 可以解释为第 \( i \) 个约束对应资源的“影子价格“（Shadow Price），即该资源增加一个单位时目标函数的改善量。这为管理决策提供了重要的经济信息。

灵敏度分析

分析目的

灵敏度分析研究当模型参数发生变化时，最优解如何变化。这对于：

• 评估最优解的稳定性
• 确定关键参数
• 为决策提供更丰富的信息

具有重要意义。

主要内容

目标函数系数变化：

当目标函数系数 \( c_j \) 在某个范围内变化时，最优基不变（即最优解的结构不变）。对于非基变量 \( x_j \)，其系数可以在使检验数保持非负的范围内变化。

右端项变化：

当约束右端项 \( b_i \) 变化时，基保持可行的条件为 \( B^{-1}\mathbf{b} \geq \mathbf{0} \)。在此范围内，最优目标函数值的变化为：

\[ \Delta z = y_i^* \cdot \Delta b_i \]

其中 \( y_i^* \) 就是对偶变量值（影子价格）。

新增变量分析：

若引入新变量 \( x_{n+1} \)，计算其检验数 \( \sigma_{n+1} = c_{n+1} - \mathbf{c}B^T B^{-1} \mathbf{a}{n+1} \)。若 \( \sigma_{n+1} < 0 \)（最小化问题），则该变量应进入基中，即当前解不再最优。

新增约束分析：

检验当前最优解是否满足新约束。若满足，则最优解不变；否则需要重新求解。

实际案例分析：生产计划问题

问题描述

某工厂生产甲、乙两种产品，需要使用 A、B、C 三种原材料。已知信息如下：

甲产品单位利润 50 元，乙产品单位利润 100 元。问：如何安排生产计划使总利润最大？

数学建模

设 \( x_1 \) 为甲产品产量，\( x_2 \) 为乙产品产量。

\[ \max \quad z = 50x_1 + 100x_2 \]

\[ \text{s.t.} \begin{cases} x_1 + 2x_2 \leq 8 \quad \text{(原材料 A)} \\ 4x_1 \leq 16 \quad \text{(原材料 B)} \\ 4x_2 \leq 12 \quad \text{(原材料 C)} \\ x_1, x_2 \geq 0 \end{cases} \]

图解法求解

简化约束：\( x_1 + 2x_2 \leq 8 \)，\( x_1 \leq 4 \)，\( x_2 \leq 3 \)。

可行域顶点：

• \( O = (0, 0) \)：\( z = 0 \)
• \( A = (4, 0) \)：\( z = 200 \)
• \( B = (4, 2) \)：\( z = 50(4) + 100(2) = 400 \)
• \( C = (2, 3) \)：\( z = 50(2) + 100(3) = 400 \)
• \( D = (0, 3) \)：\( z = 300 \)

顶点 B 和 C 的目标函数值相同，均为 400。这意味着线段 BC 上所有点都是最优解（多重最优解）。

最优解为 \( x_1 = 4, x_2 = 2 \) 或 \( x_1 = 2, x_2 = 3 \)（以及它们的凸组合），最大利润为 400 元。

单纯形法求解

转化为标准形式（最小化）：

\[ \min \quad z’ = -50x_1 - 100x_2 \]

添加松弛变量 \( s_1, s_2, s_3 \)：

\[ x_1 + 2x_2 + s_1 = 8, \quad 4x_1 + s_2 = 16, \quad 4x_2 + s_3 = 12 \]

初始单纯形表：

\( x_2 \) 进基（\( \sigma_2 = -100 \)），比值：\( 8/2=4, -, 12/4=3 \)，第三行 \( s_3 \) 出基。

第一次迭代：

\( x_1 \) 进基（\( \sigma_1 = -50 \)），比值：\( 2/1=2, 16/4=4, - \)，第一行 \( s_1 \) 出基。

第二次迭代：

所有检验数 \( \geq 0 \)，达到最优。注意 \( \sigma_5 = 0 \)（对应 \( s_3 \)），说明存在多重最优解。

最优解：\( x_1 = 2, x_2 = 3 \)，\( z’_{\min} = -400 \)，即最大利润 \( z = 400 \) 元。

对偶分析与灵敏度

对偶变量（影子价格）从最终单纯形表中读出：

• 原材料 A：\( y_1 = 50 \)（对应 \( s_1 \) 列的检验数）
• 原材料 B：\( y_2 = 0 \)（对应 \( s_2 \) 列的检验数）
• 原材料 C：\( y_3 = 0 \)（对应 \( s_3 \) 列的检验数）

经济解释：

• 原材料 A 的影子价格为 50，表示每增加 1 kg 原材料 A 的供应，最大利润可增加 50 元
• 原材料 B 的影子价格为 0，表示原材料 B 有剩余（\( s_2 = 8 > 0 \)），增加供应不会改善利润
• 原材料 C 的影子价格为 0，但 \( s_3 \) 的检验数为 0 意味着约束处于临界状态

Python 代码实现

使用 scipy.optimize.linprog

import numpy as np
from scipy.optimize import linprog

# ============================================================
# 生产计划问题求解
# max z = 50*x1 + 100*x2
# s.t. x1 + 2*x2 <= 8
#      4*x1       <= 16
#           4*x2  <= 12
#      x1, x2 >= 0
# ============================================================

# linprog 求解的是最小化问题，因此取目标函数系数的负值
c = [-50, -100]

# 不等式约束 A_ub @ x <= b_ub
A_ub = [
    [1, 2],    # 原材料A约束
    [4, 0],    # 原材料B约束
    [0, 4],    # 原材料C约束
]
b_ub = [8, 16, 12]

# 变量范围
x1_bounds = (0, None)
x2_bounds = (0, None)
bounds = [x1_bounds, x2_bounds]

# 求解
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

# 输出结果
print("=" * 50)
print("生产计划线性规划求解结果")
print("=" * 50)
print(f"求解状态: {result.message}")
print(f"甲产品产量 x1 = {result.x[0]:.4f}")
print(f"乙产品产量 x2 = {result.x[1]:.4f}")
print(f"最大利润 z = {-result.fun:.4f} 元")
print("=" * 50)

获取对偶信息与灵敏度分析

from scipy.optimize import linprog
import numpy as np

# 问题数据
c = [-50, -100]
A_ub = [[1, 2], [4, 0], [0, 4]]
b_ub = [8, 16, 12]
bounds = [(0, None), (0, None)]

# 使用 HiGHS 求解器（scipy >= 1.9）
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

# 对偶信息
# result.ineqlin 包含不等式约束的对偶信息
print("对偶变量（影子价格）:")
print(f"  原材料A: {-result.ineqlin.marginals[0]:.4f}")
print(f"  原材料B: {-result.ineqlin.marginals[1]:.4f}")
print(f"  原材料C: {-result.ineqlin.marginals[2]:.4f}")

# 约束松弛量
slack = b_ub - np.array(A_ub) @ result.x
print("\n约束松弛量（资源剩余）:")
print(f"  原材料A剩余: {slack[0]:.4f} kg")
print(f"  原材料B剩余: {slack[1]:.4f} kg")
print(f"  原材料C剩余: {slack[2]:.4f} kg")

更复杂的示例：饮食问题

from scipy.optimize import linprog
import numpy as np

# ============================================================
# 饮食问题（Diet Problem）
# 目标：以最低成本满足营养需求
#
# 食物：面包、牛奶、鸡蛋、牛肉
# 营养素：蛋白质、钙、热量
# ============================================================

# 食物单价（元/单位）
# 面包  牛奶  鸡蛋  牛肉
cost = [2, 3.5, 2.5, 12]

# 每单位食物含营养素量
# 行：蛋白质(g)、钙(mg)、热量(kcal)
# 列：面包、牛奶、鸡蛋、牛肉
nutrition = np.array([
    [3,   8,   6,  20],   # 蛋白质
    [20, 120,  25,  10],  # 钙
    [250, 150, 75, 180],  # 热量
])

# 每日最低需求
min_requirement = [50, 800, 2000]  # 蛋白质50g, 钙800mg, 热量2000kcal

# 每日最高摄入（上限约束）
max_intake = [4000]  # 热量不超过4000kcal

# 构建约束
# 下限约束：nutrition @ x >= min_requirement
# 转化为：-nutrition @ x <= -min_requirement
A_ub_lower = -nutrition
b_ub_lower = -np.array(min_requirement)

# 上限约束：热量 <= 4000
A_ub_upper = [nutrition[2]]  # 热量行
b_ub_upper = max_intake

# 合并不等式约束
A_ub = np.vstack([A_ub_lower, A_ub_upper])
b_ub = np.concatenate([b_ub_lower, b_ub_upper])

# 变量范围：每种食物消费 0-10 单位
bounds = [(0, 10)] * 4

# 求解
result = linprog(cost, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

# 输出
food_names = ['面包', '牛奶', '鸡蛋', '牛肉']
nutrient_names = ['蛋白质(g)', '钙(mg)', '热量(kcal)']

print("=" * 50)
print("饮食问题线性规划求解结果")
print("=" * 50)
print(f"求解状态: {result.message}")
print(f"\n最低日花费: {result.fun:.2f} 元\n")

print("每日食物消费计划:")
for i, name in enumerate(food_names):
    print(f"  {name}: {result.x[i]:.4f} 单位 (花费 {cost[i]*result.x[i]:.2f} 元)")

print("\n营养素摄入验证:")
actual_nutrition = nutrition @ result.x
for i, name in enumerate(nutrient_names):
    status = "满足" if actual_nutrition[i] >= min_requirement[i] - 1e-6 else "不足"
    print(f"  {name}: {actual_nutrition[i]:.2f} (需求: {min_requirement[i]}) [{status}]")

封装为通用求解函数

from scipy.optimize import linprog
import numpy as np


def solve_lp(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None,
             bounds=None, maximize=False, verbose=True):
    """
    通用线性规划求解函数

    Parameters
    ----------
    c : array_like
        目标函数系数向量
    A_ub : array_like, optional
        不等式约束矩阵（<= 形式）
    b_ub : array_like, optional
        不等式约束右端项
    A_eq : array_like, optional
        等式约束矩阵
    b_eq : array_like, optional
        等式约束右端项
    bounds : list of tuples, optional
        变量上下界
    maximize : bool
        是否为最大化问题
    verbose : bool
        是否打印详细结果

    Returns
    -------
    dict : 包含最优解、最优值、对偶信息等
    """
    c = np.array(c, dtype=float)
    obj_c = -c if maximize else c

    result = linprog(obj_c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
                     bounds=bounds, method='highs')

    if not result.success:
        if verbose:
            print(f"求解失败: {result.message}")
        return None

    optimal_value = -result.fun if maximize else result.fun

    output = {
        'x': result.x,
        'optimal_value': optimal_value,
        'status': result.message,
        'slack': result.slack if hasattr(result, 'slack') else None,
    }

    if verbose:
        print(f"最优值: {optimal_value:.6f}")
        print(f"最优解: {result.x}")
        if hasattr(result, 'ineqlin') and result.ineqlin is not None:
            print(f"影子价格: {result.ineqlin.marginals}")

    return output


# 使用示例
if __name__ == '__main__':
    # 生产计划问题
    result = solve_lp(
        c=[50, 100],
        A_ub=[[1, 2], [4, 0], [0, 4]],
        b_ub=[8, 16, 12],
        bounds=[(0, None), (0, None)],
        maximize=True
    )

应用注意事项与局限性

建模注意事项

1. 线性假设的合理性：线性规划要求目标函数和约束条件都是线性的。在建模前必须验证：

   • 成本/收益是否与产量成正比
   • 资源消耗是否与生产量成线性关系
   • 是否存在规模效应或边际递减

2. 变量连续性假设：线性规划的解通常是连续值。若决策变量必须取整数（如设备台数、人员数量），应使用整数规划。

3. 确定性假设：标准线性规划假设所有参数已知且确定。若参数存在不确定性，需考虑鲁棒优化或随机规划。

4. 约束条件的完整性：确保所有实际限制都被纳入模型，遗漏关键约束会导致不切实际的解。

求解中的常见问题

1. 退化问题：当基本可行解中某些基变量为零时，单纯形法可能出现循环。解决方法包括 Bland 规则或字典序法。

2. 数值稳定性：系数矩阵的条件数过大会导致数值误差。建议：

   • 对变量和约束进行适当缩放
   • 使用高精度求解器
   • 检验解的可行性

3. 大规模问题：变量和约束数量庞大时，标准单纯形法效率下降。可采用：

   • 修正单纯形法
   • 内点法（如 method='highs' 中的 IPM 选项）
   • 列生成或分解技术

局限性

1. 线性近似的局限：现实中许多关系是非线性的（如规模经济、折扣定价），线性规划只能给出近似结果。

2. 静态模型的局限：标准线性规划是单期模型，不能直接处理多阶段动态决策问题。

3. 单目标的局限：线性规划只能优化单个目标函数。多目标情况需使用多目标规划方法（如加权法、约束法）。

4. 解的敏感性：参数的微小变化可能导致最优解跳变（顶点之间切换），需要结合灵敏度分析来评估解的稳健性。

与其他方法的关系

实用建议

• 先用小规模算例验证模型的正确性
• 善用对偶理论和灵敏度分析理解解的经济含义
• 关注约束的紧与松，识别瓶颈资源
• 在实际应用中，结合灵敏度分析进行决策，而非仅依赖单一最优解
• 对于大规模实际问题，考虑使用专业求解器（如 Gurobi、CPLEX）以获得更好的性能