运筹优化算法学习实践：从入门到精通

目录

• 第一篇：基础概念与总览
  • 1.1 什么是运筹优化
  • 1.2 优化问题的一般形式
  • 1.3 优化问题的分类体系
  • 1.4 发展历史与技术脉络
• 第二篇：数学基础
  • 2.1 线性代数基础
  • 2.2 微积分与最优化条件
  • 2.3 凸分析基础
  • 2.4 对偶理论
• 第三篇：线性规划（LP）
  • 3.1 线性规划的标准形式
  • 3.2 图解法与几何直觉
  • 3.3 单纯形法
  • 3.4 内点法
  • 3.5 灵敏度分析
  • 3.6 线性规划的Python实现
• 第四篇：整数规划与混合整数规划（IP/MIP）
  • 4.1 整数规划概述
  • 4.2 分支定界法
  • 4.3 割平面法
  • 4.4 分支切割法
  • 4.5 拉格朗日松弛
  • 4.6 Benders分解
  • 4.7 经典建模技巧
  • 4.8 Python实现
• 第五篇：非线性规划（NLP）
  • 5.1 非线性规划概述
  • 5.2 无约束优化方法
  • 5.3 约束优化方法
  • 5.4 二次规划（QP）
  • 5.5 Python实现
• 第六篇：动态规划（DP）
  • 6.1 动态规划基本思想
  • 6.2 最优子结构与重叠子问题
  • 6.3 状态转移方程
  • 6.4 经典问题
  • 6.5 近似动态规划
• 第七篇：图论与网络优化
  • 7.1 图论基础
  • 7.2 最短路径问题
  • 7.3 最小生成树
  • 7.4 网络流问题
  • 7.5 匹配与指派问题
  • 7.6 旅行商问题（TSP）
  • 7.7 车辆路径问题（VRP）
• 第八篇：启发式与元启发式算法
  • 8.1 启发式算法概述
  • 8.2 贪心算法
  • 8.3 局部搜索
  • 8.4 模拟退火（SA）
  • 8.5 遗传算法（GA）
  • 8.6 粒子群优化（PSO）
  • 8.7 蚁群算法（ACO）
  • 8.8 禁忌搜索（TS）
  • 8.9 差分进化（DE）
  • 8.10 自适应大邻域搜索（ALNS）
  • 8.11 超启发式框架
• 第九篇：随机优化与鲁棒优化
  • 9.1 不确定性优化概述
  • 9.2 随机规划
  • 9.3 鲁棒优化
  • 9.4 分布鲁棒优化
• 第十篇：多目标优化
  • 10.1 多目标优化基础
  • 10.2 Pareto最优
  • 10.3 求解方法
  • 10.4 NSGA-II算法
• 第十一篇：约束规划（CP）
  • 11.1 约束规划概述
  • 11.2 约束传播与回溯搜索
  • 11.3 全局约束
  • 11.4 CP与MIP的对比
• 第十二篇：求解器与工具链
  • 12.1 商业求解器
  • 12.2 开源求解器
  • 12.3 建模语言与接口
  • 12.4 求解器选型指南
• 第十三篇：实战案例
  • 13.1 生产调度问题
  • 13.2 物流配送优化
  • 13.3 投资组合优化
  • 13.4 选址问题
  • 13.5 排班问题
  • 13.6 供应链网络设计
• 第十四篇：前沿方向与学习路线
  • 14.1 机器学习+优化
  • 14.2 在线优化与实时决策
  • 14.3 量子优化
  • 14.4 学习路线推荐
  • 14.5 推荐资源

---

第一篇：基础概念与总览

1.1 什么是运筹优化

运筹学（Operations Research, OR） 是一门应用数学学科，起源于二战时期的军事行动研究，旨在通过数学建模、统计分析和优化方法为复杂系统的决策问题提供最优或近优解。

优化（Optimization） 是运筹学的核心，简单来说就是：

在给定约束条件下，从所有可行方案中找到使目标函数最优（最大或最小）的那个方案。

生活中的优化问题举例

场景	决策变量	目标函数	约束条件
物流配送	每辆车的行驶路线	总运输成本最小	车辆容量、时间窗、客户需求
生产计划	各产品生产数量	利润最大	原料、产能、市场需求
投资组合	各资产投资比例	收益最大/风险最小	预算、法规、风险承受
排班调度	员工班次安排	人力成本最小	劳动法、技能匹配、偏好
网络路由	数据包转发路径	延迟最小/吞吐最大	带宽、拓扑约束

运筹优化的核心价值

1. 降本增效：通过数学方法找到最优资源配置方案
2. 科学决策：将主观判断转化为可量化、可复现的决策过程
3. 应对复杂性：处理人类直觉难以解决的大规模组合问题
4. 风险量化：在不确定性环境下做出稳健决策

1.2 优化问题的一般形式

任何一个优化问题都可以用以下三要素来描述：

$$\begin{array}{rlr}\min & f(x)& \text{（目标函数 Objective Function）}\\ \text{s.t.}& g_i(x)\le 0,i=1,...,m& \text{（不等式约束 Inequality Constraints）}\\ & h_j(x)=0,j=1,...,p& \text{（等式约束 Equality Constraints）}\\ & x\in \mathcal{X}& \text{（变量域 Variable Domain）}\end{array}$$

其中：

• $x=(x_1,x_2,...,x_n{)}^T$ 是决策变量（Decision Variables），是我们可以控制和调整的量
• $f(x)$ 是目标函数（Objective Function），是我们要最大化或最小化的指标
• $g_i(x)\le 0$ 是不等式约束
• $h_j(x)=0$ 是等式约束
• $\mathcal{X}$ 是变量的取值范围（如整数、非负等）

可行域（Feasible Region） 是所有满足约束条件的 $x$ 的集合：

F = { x ∈ X ∣ g i ( x ) ≤ 0 , ∀ i ;    h j ( x ) = 0 , ∀ j } \mathcal{F} = \{x \in \mathcal{X} \mid g_i(x) \leq 0, \forall i; \; h_j(x) = 0, \forall j\} F={x∈X∣gi​(x)≤0,∀i;hj​(x)=0,∀j}

最优解（Optimal Solution） $x^{\ast }$ 是在可行域中使目标函数达到最小值的点：

$$f(x^{\ast })\le f(x),\forall x\in \mathcal{F}$$

解的概念层次

┌─────────────────────────────────────────┐
│           全局最优解 (Global Optimal)      │  ← 理想目标
│   ┌─────────────────────────────────┐   │
│   │      局部最优解 (Local Optimal)    │   │  ← 实际常得到
│   │   ┌─────────────────────────┐   │   │
│   │   │    可行解 (Feasible)      │   │   │  ← 满足约束
│   │   │                         │   │   │
│   │   └─────────────────────────┘   │   │
│   └─────────────────────────────────┘   │
└─────────────────────────────────────────┘

1.3 优化问题的分类体系

优化问题可以从多个维度进行分类：

按决策变量类型分类

类型	说明	典型问题
连续优化	变量取实数值	曲线拟合、参数估计
离散/整数优化	变量取整数值	选址、调度、指派
混合优化	部分连续、部分离散	生产排程、投资组合

按目标函数与约束的性质分类

类型	说明	难度
线性规划（LP）	目标和约束均为线性	★☆☆ 多项式可解
二次规划（QP）	目标为二次函数，约束线性	★★☆
凸优化	目标凸、可行域凸	★★☆ 有全局最优保证
非线性规划（NLP）	含非线性函数	★★★ 可能陷入局部最优
整数规划（IP）	变量需取整	★★★★ NP-hard
混合整数规划（MIP）	部分变量整数	★★★★★

按目标数量分类

类型	说明
单目标优化	只有一个优化目标
多目标优化	多个可能冲突的目标，需寻找Pareto最优

按信息确定性分类

类型	说明
确定性优化	所有参数已知确定
随机优化	部分参数为随机变量
鲁棒优化	参数在不确定集内变化
模糊优化	参数具有模糊性

按时间维度分类

类型	说明
静态优化	单次决策
动态优化/动态规划	多阶段序贯决策
在线优化	信息逐步揭示，需实时决策

计算复杂性视角

易解问题 (P)                          难解问题 (NP-hard)
┌──────────────┐                    ┌──────────────────┐
│ 线性规划(LP)  │                    │ 旅行商问题(TSP)   │
│ 最短路径      │    ┌──────────┐   │ 背包问题           │
│ 最小生成树    │    │ NP完全问题│   │ 调度问题           │
│ 网络流        │    └──────────┘   │ 图着色             │
│ 匹配问题      │                    │ 整数规划(IP)       │
└──────────────┘                    │ 集合覆盖           │
                                    └──────────────────┘

1.4 发展历史与技术脉络

运筹优化发展时间线

1939年  ── 二战军事运筹研究（英国雷达部署优化）
    │
1947年  ── George Dantzig 发明单纯形法（线性规划的里程碑）
    │
1950s   ── 线性规划理论成熟，对偶理论建立
    │       动态规划理论（Richard Bellman）
    │
1960s   ── 整数规划的分支定界法
    │       非线性规划的梯度方法
    │
1970s   ── 网络优化算法
    │       拉格朗日松弛方法
    │
1980s   ── 内点法诞生（Karmarkar, 1984）
    │       计算复杂性理论建立
    │       第一批商业求解器出现
    │
1990s   ── 元启发式算法蓬勃发展
    │       (遗传算法、模拟退火、蚁群算法等)
    │       约束规划兴起
    │       随机优化理论
    │
2000s   ── 商业求解器性能飞跃（CPLEX, Gurobi）
    │       开源求解器发展（GLPK, CBC）
    │       鲁棒优化理论
    │
2010s   ── 大规模优化 + 云计算
    │       机器学习 + 优化的交叉
    │       深度学习求解组合优化
    │
2020s   ── AI驱动的求解器（学习增强分支定界）
    │       GPU加速优化
    │       大语言模型辅助建模
    │       量子优化初步探索

关键人物与贡献

年代	人物	贡献
1947	George Dantzig	单纯形法、线性规划
1951	Richard Bellman	动态规划原理
1962	Land & Doig	分支定界法求解整数规划
1979	Leonid Khachiyan	椭球法（首个LP多项式算法）
1984	Narendra Karmarkar	内点法
1985	Fred Glover	禁忌搜索
1989	John Holland	遗传算法的系统理论
1992	Marco Dorigo	蚁群优化算法

---

第二篇：数学基础

2.1 线性代数基础

运筹优化的核心数学工具是线性代数。优化问题的建模和求解都离不开矩阵和向量运算。

向量与矩阵表示

优化问题通常用向量-矩阵形式紧凑表达：

$$\begin{array}{rl}\min & c^{T}x\\ \text{s.t.}& Ax\le b\\ & A_{eq}x=b_{eq}\\ & l\le x\le u\end{array}$$

其中 $c\in {\mathbb{R}}^n$，$A\in {\mathbb{R}}^{m\times n}$，$b\in {\mathbb{R}}^m$。

重要概念

• 仿射组合：$x=\sum {\lambda }_i{x}_i$，其中 $\sum {\lambda }_i=1$
• 凸组合：仿射组合 + ${\lambda }_i\ge 0$
• 凸包：所有凸组合构成的集合， conv ( S ) = { ∑ λ i x i ∣ λ i ≥ 0 , ∑ λ i = 1 , x i ∈ S } \text{conv}(S) = \{\sum \lambda_i x_i \mid \lambda_i \geq 0, \sum \lambda_i = 1, x_i \in S\} conv(S)={∑λi​xi​∣λi​≥0,∑λi​=1,xi​∈S}

2.2 微积分与最优化条件

极值条件

对于无约束优化问题 $\min f(x)$：

一阶必要条件（FONC）：
$$\nabla f(x^{\ast })=0$$

二阶充分条件（SOSC）：
$${\nabla }^{2}f(x^{\ast })\succ 0\text{（Hessian正定）}$$

梯度与Hessian矩阵

$$\nabla f(x)=\left(\begin{array}{c}\frac{\partial f}{\partial x_1}\\ ⋮\\ \frac{\partial f}{\partial x_n}\end{array}\right),{\nabla }^{2}f(x)=\left(\begin{array}{ccc}\frac{{\partial }^{2}f}{\partial x_1^2}& \cdots & \frac{{\partial }^{2}f}{\partial x_1\partial x_n}\\ ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f}{\partial x_n\partial x_1}& \cdots & \frac{{\partial }^{2}f}{\partial x_n^2}\end{array}\right)$$

梯度的含义：梯度指向函数值增长最快的方向，负梯度方向是下降最快的方向。

KKT条件

对于约束优化问题，KKT（Karush-Kuhn-Tucker）条件是最优解的一阶必要条件：

$$\begin{array}{rlr}& \nabla f(x^{\ast })+\sum _{i=1}^m{\mu }_i\nabla g_i(x^{\ast })+\sum _{j=1}^p{\lambda }_j\nabla h_j(x^{\ast })=0& \text{（平稳性）}\\ & g_i(x^{\ast })\le 0,i=1,...,m& \text{（原始可行性）}\\ & h_j(x^{\ast })=0,j=1,...,p& \text{（原始可行性）}\\ & {\mu }_i\ge 0,i=1,...,m& \text{（对偶可行性）}\\ & {\mu }_i{g}_i(x^{\ast })=0,i=1,...,m& \text{（互补松弛性）}\end{array}$$

互补松弛性的含义：如果某个不等式约束在最优解处是严格不等式（$g_i(x^{\ast })<0$），则对应的乘子 ${\mu }_i=0$；反之如果 ${\mu }_i>0$，则 $g_i(x^{\ast })=0$（约束取紧）。

2.3 凸分析基础

凸优化是运筹优化中最重要的子类，因为它保证局部最优即全局最优。

凸集

集合 $C$ 是凸集，当且仅当：对任意 $x,y\in C$ 和 $\lambda \in [0,1]$：
$$\lambda x+(1-\lambda )y\in C$$

常见凸集：超平面、半空间、多面体、椭球、单纯形、锥。

凸函数

函数 $f$ 是凸函数，当且仅当：对任意 $x,y$ 和 $\lambda \in [0,1]$：
$$f(\lambda x+(1-\lambda )y)\le \lambda f(x)+(1-\lambda )f(y)$$

判断方法：

1. 定义法
2. 一阶条件：$f(y)\ge f(x)+\nabla f(x{)}^T(y-x)$
3. 二阶条件：${\nabla }^{2}f(x)⪰0$（Hessian半正定）

凸优化问题

$$\begin{array}{rl}\min & f(x)\text{（}f\text{ 是凸函数）}\\ \text{s.t.}& g_i(x)\le 0,i=1,...,m\text{（}{g}_i\text{ 是凸函数）}\\ & Ax=b\end{array}$$

核心性质：凸优化问题的任何局部最优解都是全局最优解。这是凸优化在工程实践中如此重要的根本原因。

2.4 对偶理论

对偶理论是线性规划和凸优化的核心理论工具。

拉格朗日对偶

对于原问题（Primal）：
min ⁡ x { f ( x ) ∣ g i ( x ) ≤ 0 , h j ( x ) = 0 } \min_x \{f(x) \mid g_i(x) \leq 0, h_j(x) = 0\} xmin​{f(x)∣gi​(x)≤0,hj​(x)=0}

定义拉格朗日函数：
$$L(x,\mu ,\lambda )=f(x)+\sum _i{\mu }_i{g}_i(x)+\sum _j{\lambda }_j{h}_j(x)$$

对偶函数：
$$q(\mu ,\lambda )=\underset{x}{\inf }L(x,\mu ,\lambda )$$

对偶问题（Dual）：
$$\underset{\mu \ge 0,\lambda }{\max }q(\mu ,\lambda )$$

弱对偶与强对偶

• 弱对偶：$q({\mu }^{\ast },{\lambda }^{\ast })\le f(x^{\ast })$（对偶最优值 ≤ 原问题最优值）
• 强对偶：$q({\mu }^{\ast },{\lambda }^{\ast })=f(x^{\ast })$（两者相等）
  • 对于线性规划，强对偶总是成立（如果原问题有有限最优解）
  • 对于凸优化，在Slater条件下成立

对偶间隙

$$\text{对偶间隙}=f(x^{\ast })-q({\mu }^{\ast },{\lambda }^{\ast })\ge 0$$

对偶间隙为零时，强对偶成立，对偶最优解提供了原问题最优值的一个紧界。

---

第三篇：线性规划（LP）

3.1 线性规划的标准形式

线性规划是最基本、最成熟的优化问题类型：

$$\begin{array}{rl}\min & c^{T}x\\ \text{s.t.}& Ax=b\\ & x\ge 0\end{array}$$

任何LP都可以通过以下变换转化为标准形式：

1. max → min：$\max c^{T}x=-\min (-c^{T}x)$
2. 不等式 → 等式：引入松弛变量 $s_i\ge 0$，使 $a_i^{T}x\le b_i$ 变为 $a_i^{T}x+s_i=b_i$
3. 无约束变量：$x_j=x_j^{+}-x_j^{-}$，其中 $x_j^{+},x_j^{-}\ge 0$
4. 负右端项：等式两边乘以 -1

松弛形式与基

在标准形式 $Ax=b,x\ge 0$ 中（$m$ 个约束，$n$ 个变量，$m\le n$）：

• 基变量（Basic Variables）：$m$ 个，由基矩阵 $B$ 决定
• 非基变量（Non-basic Variables）：$n-m$ 个，设为0
• 基本可行解（BFS）：$x_B=B^{-1}b\ge 0$，$x_N=0$

关键定理：如果线性规划有最优解，则一定存在一个基本可行解是最优解。

3.2 图解法与几何直觉

对于二维LP问题，可以通过画图直观理解：

考虑问题：
min  -x₁ - 2x₂
s.t. x₁ + x₂  ≤ 4     ①
     x₁         ≤ 3     ②
         x₂     ≤ 2     ③
     x₁, x₂    ≥ 0

    x₂
  3 ┤
    │
  2 ┤────────●●●          ③
    │       ╱   ●
    │      ╱     ●
    │     ╱  可行域 ●
    │    ╱         ●
  0 ┤──●───────────●─── x₁
    0  1  2  3  4
              ②

几何直觉：

• 可行域是一个多面体（Polytope）
• 目标函数的等值线是一组平行直线
• 最优解出现在多面体的某个**顶点（Vertex）**上
• 这就是为什么单纯形法沿着顶点搜索

3.3 单纯形法

单纯形法（Simplex Method）由George Dantzig于1947年提出，是求解线性规划最经典的方法。

核心思想

从一个基本可行解（多面体的一个顶点）出发，沿着多面体的边移动到相邻的、目标函数值更好的顶点，直到找到最优解。

算法流程

1. 初始化：找到一个初始基本可行解
2. 最优性检验：计算检验数 σⱼ = cⱼ - cᵦB⁻¹aⱼ
   - 如果所有 σⱼ ≥ 0 → 当前解最优，停止
   - 否则 → 选择 σⱼ < 0 中最负的变量 xⱼ 入基
3. 确定出基变量：最小比值测试
   θ = min{bᵢ/aᵢⱼ | aᵢⱼ > 0}
   对应的基变量出基
4. 更新：旋转（pivot），转到步骤2

单纯形表

以一个具体例子演示：

问题：
$$\min -x_1-2x_2\text{s.t.}{x}_1+x_2+s_1=4,x_1+s_2=3,x_2+s_3=2,x\ge 0$$

初始单纯形表：

基变量	$x_1$	$x_2$	$s_1$	$s_2$	$s_3$	RHS
$s_1$	1	1	1	0	0	4
$s_2$	1	0	0	1	0	3
$s_3$	0	1	0	0	1	2
$z$	1	2	0	0	0	0

检验数 ${\sigma }_{x_2}=2$（最正），选择 $x_2$ 入基。最小比值测试：$4/1=4$，$2/1=2$，选 $s_3$ 出基。经过旋转迭代…

单纯形法的效率

• 最好情况：$O(m)$ 次迭代
• 最坏情况：指数级（Klee-Minty反例）
• 实际表现：通常在 $O(m)$ 到 $O(m^2)$ 次迭代内收敛，对大规模问题非常高效

退化与循环

当基本可行解中某个基变量取值为0时，发生退化。理论上退化可能导致循环（Cycling），但实际中很少发生。Bland规则等方法可以避免循环。

3.4 内点法

内点法（Interior Point Method）由Karmarkar于1984年提出，是从可行域内部逼近最优解的方法。

与单纯形法的对比

特性	单纯形法	内点法
搜索路径	沿多面体边沿顶点	穿过可行域内部
迭代次数	指数最坏，实际很好	$O(\sqrt{n}\log (1/ϵ))$
每次迭代复杂度	$O(mn)$	$O(n^3)$（稠密）
大规模问题	适合稀疏问题	适合大规模问题
热启动	支持好	支持较差

原始-对偶内点法

最常用的内点法形式。引入障碍函数（Barrier）：

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

其中 $\mu >0$ 是障碍参数，随着 $\mu \to 0$，解趋向原问题最优解。

算法框架：

1. 初始化：选择可行内点 x⁰ > 0
2. 对每个障碍参数 μₖ：
   a. 用牛顿法求解障碍问题的KKT系统
   b. 更新 x ← x + αΔx（适当步长）
3. 减小 μ：μₖ₊₁ = σμₖ（0 < σ < 1）
4. 重复直到 μ 足够小，满足精度要求

3.5 灵敏度分析

灵敏度分析研究当问题参数变化时，最优解如何变化。

影子价格（Shadow Price）

对偶变量 ${\lambda }^{\ast }$ 的值反映了约束右端项的影子价格：

• ${\lambda }_i^{\ast }>0$：增加第 $i$ 个约束的右端值（放松约束）可以改善目标值
• ${\lambda }_i^{\ast }=0$：该约束不紧，改变它不影响最优值
• $|{\lambda }_i^{\ast }|$ 表示右端项每变化1个单位，目标值变化 $|{\lambda }_i^{\ast }|$ 个单位

可变范围（Range of Optimality / Feasibility）

• 最优性范围：目标函数系数 $c_j$ 在什么范围内变化时，当前基仍然最优
• 可行性范围：右端项 $b_i$ 在什么范围内变化时，当前基仍可行

3.6 线性规划的Python实现

使用SciPy

from scipy.optimize import linprog

# 问题：min -x₁ - 2x₂
# s.t.  x₁ + x₂ ≤ 4
#       x₁       ≤ 3
#            x₂  ≤ 2
#       x₁, x₂  ≥ 0

c = [-1, -2]                    # 目标函数系数（min c^Tx）
A_ub = [[1, 1], [1, 0], [0, 1]] # 不等式约束 Ax ≤ b
b_ub = [4, 3, 2]
bounds = [(0, None), (0, None)]  # 变量下界

result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')
print(f"最优解: x₁={result.x[0]:.2f}, x₂={result.x[1]:.2f}")
print(f"最优值: {-result.fun:.2f}")  # 注意取反（因为我们用了max变min）

使用PuLP

from pulp import *

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

# 决策变量
x1 = LpVariable("x1", lowBound=0)
x2 = LpVariable("x2", lowBound=0)

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

# 约束条件
prob += x1 + x2 <= 4
prob += x1 <= 3
prob += x2 <= 2

# 求解
prob.solve()

print(f"最优解: x1={value(x1)}, x2={value(x2)}")
print(f"最优值: {value(prob.objective)}")

使用Google OR-Tools

from ortools.linear_solver import pywraplp

solver = pywraplp.Solver.CreateSolver('GLOP')

# 变量
x1 = solver.NumVar(0, solver.infinity(), 'x1')
x2 = solver.NumVar(0, solver.infinity(), 'x2')

# 约束
solver.Add(x1 + x2 <= 4)
solver.Add(x1 <= 3)
solver.Add(x2 <= 2)

# 目标
solver.Maximize(x1 + 2*x2)

status = solver.Solve()
if status == pywraplp.Solver.OPTIMAL:
    print(f"最优解: x1={x1.solution_value()}, x2={x2.solution_value()}")
    print(f"最优值: {solver.Objective().Value()}")

---

第四篇：整数规划与混合整数规划（IP/MIP）

4.1 整数规划概述

当决策变量必须取整数值时，问题变为整数规划：

min ⁡ c T x s.t. A x ≤ b x j ∈ Z （整数约束） x j ∈ { 0 , 1 } （0-1变量，二进制约束） \begin{aligned} \min \quad & c^T x \\ \text{s.t.} \quad & Ax \leq b \\ & x_j \in \mathbb{Z} \quad \text{（整数约束）} \\ & x_j \in \{0, 1\} \quad \text{（0-1变量，二进制约束）} \end{aligned} mins.t.​cTxAx≤bxj​∈Z（整数约束）xj​∈{0,1}（0-1变量，二进制约束）​

为什么整数规划比线性规划难得多？

1. NP-hard：一般整数规划是NP-hard问题
2. 可行域不连续：LP的可行域是连续多面体，IP的可行域是离散点集
3. LP松弛的差距：LP松弛（去掉整数约束）的最优解可能与IP最优解差距很大

整数规划的应用场景

问题类型	说明
选址问题	在哪些候选位置建仓库（0-1变量）
调度问题	任务分配到哪些时间槽（整数变量）
背包问题	选哪些物品装入背包（0-1变量）
网络设计	建设哪些网络连接（0-1变量）
投资决策	选哪些项目投资（0-1变量）

4.2 分支定界法

分支定界法（Branch and Bound, B&B）是求解整数规划最基础也最重要的精确算法。

核心思想

将整数约束松弛为连续约束得到LP松弛，如果松弛问题的最优解不满足整数约束，则选取一个非整数变量进行分支（分成两个子问题），通过定界（计算上下界）来剪枝（排除不可能包含最优解的分支）。

算法流程

分支定界法：

输入：整数规划问题 P
输出：最优解 x*

1. 初始化：
   - 上界 UB = +∞（最小化问题）
   - 下界 LB = -∞
   - 活动节点列表 = {P}
   - 最优解 x* = None

2. 当活动节点列表非空：
   a. 选择一个活动节点 Pₖ
   b. 求解 Pₖ 的LP松弛问题
   c. 如果Pₖ不可行 → 剪枝（Pruned by infeasibility）
   d. 如果松弛最优值 > UB → 剪枝（Pruned by bound）
   e. 如果松弛解满足整数约束 → 更新下界和最优解（如果更优）
      更新 LB = max(LB, 该值)，剪枝（Pruned by optimality）
   f. 否则（分支）：
      - 选一个非整数变量 xⱼ = 5.3
      - 创建两个子问题：
        P_left:  Pₖ + xⱼ ≤ 5
        P_right: Pₖ + xⱼ ≥ 6
      - 将子问题加入活动节点列表

3. 返回 x*

分支变量选择策略

• 最大分数部分：选择 $x_j-\lfloor x_j\rceil$ 最大的变量
• 伪成本法：根据历史信息估计分支后目标值的变化
• 强分支（Strong Branching）：尝试对每个候选变量分支一步，选效果最好的

节点选择策略

• 最佳优先（Best-First）：选择下界最小的节点（最有可能找到改进解）
• 深度优先（Depth-First）：尽快找到可行解，快速收紧上界
• 最佳估计（Best-Estimate）：基于启发式评估最有希望的节点

4.3 割平面法

割平面法（Cutting Plane）通过添加有效的不等式（割）来收紧LP松弛的可行域，排除非整数最优解而不排除任何整数可行解。

Gomory割

对LP松弛最优单纯形表中的某行：
$$x_i+\sum _{j\in N}{\bar{a}}_{ij}{x}_j={\bar{b}}_i$$

其中 $x_i$ 是基变量且 ${\bar{b}}_i$ 非整数，Gomory割为：
$$\sum _{j\in N}{f}_{ij}{x}_j\ge f_{i0}$$

其中 $f_{ij}={\bar{a}}_{ij}-\lfloor {\bar{a}}_{ij}\rfloor$，$f_{i0}={\bar{b}}_i-\lfloor {\bar{b}}_i\rfloor$。

常见割类型

割类型	适用场景
Gomory割	通用整数规划
覆盖割（Cover Cut）	背包约束
流覆盖割	网络设计
Clique割	冲突约束
混合整数舍入（MIR）割	混合整数规划
提升与投影（Lift-and-Project）	0-1规划

4.4 分支切割法

分支切割法（Branch and Cut） 将分支定界和割平面法结合，是现代MIP求解器（如CPLEX、Gurobi）的核心算法：

分支切割法：

在分支定界的每个节点：
1. 求解LP松弛
2. 寻找违反的割平面（分离问题）
3. 如果找到割 → 加入割约束，重新求解LP（不创建子节点）
4. 如果没有有效割 → 按分支定界的规则进行分支

现代MIP求解器的技术栈：

┌──────────────────────────────────────────┐
│          现代MIP求解器架构                  │
├──────────────────────────────────────────┤
│ 预处理（Presolve）                         │
│   ├── 固定变量、消除冗余约束               │
│   ├── 系数简化、对称性检测                 │
│   └── 提取特殊结构                        │
├──────────────────────────────────────────┤
│ 启发式找初始可行解                         │
│   ├── 贪心启发式                          │
│   ├── RINS（松弛诱导邻域搜索）             │
│   └── Feasibility Pump                    │
├──────────────────────────────────────────┤
│ 割平面生成                                │
│   ├── Gomory割、MIR割                     │
│   ├── 覆盖割、流割                         │
│   └── 分离算法（精确/启发式）              │
├──────────────────────────────────────────┤
│ 分支策略                                  │
│   ├── 强分支 / 伪成本                      │
│   ├── 多分支 / 分支外观                    │
│   └── 搜索树管理                          │
├──────────────────────────────────────────┤
│ LP求解器                                  │
│   ├── 单纯形法 / 内点法                    │
│   ├── 热启动                              │
│   └── 并行求解                            │
└──────────────────────────────────────────┘

4.5 拉格朗日松弛

当约束使得问题难以求解时，可以将"难"约束放入目标函数中，用拉格朗日乘子惩罚违反的程度。

基本思想

原问题：
z ∗ = min ⁡ { c x ∣ A x ≤ b , D x ≤ d , x ∈ X } z^* = \min\{cx \mid Ax \leq b, Dx \leq d, x \in X\} z∗=min{cx∣Ax≤b,Dx≤d,x∈X}

假设去掉约束 $Ax\le b$ 后问题变得容易求解（ X ′ = { x ∣ D x ≤ d , x ∈ X } X' = \{x \mid Dx \leq d, x \in X\} X′={x∣Dx≤d,x∈X} 容易处理），则拉格朗日松弛问题为：

z L R ( λ ) = min ⁡ { c x + λ ( A x − b ) ∣ x ∈ X ′ } z_{LR}(\lambda) = \min\{cx + \lambda(Ax - b) \mid x \in X'\} zLR​(λ)=min{cx+λ(Ax−b)∣x∈X′}

对偶问题：
$$z_{LD}=\underset{\lambda \ge 0}{\max }{z}_{LR}(\lambda )$$

次梯度优化

求解对偶问题的常用方法是次梯度优化：

λ k + 1 = max ⁡ { 0 , λ k + t k ( A x k − b ) } \lambda^{k+1} = \max\{0, \lambda^k + t_k(Ax^k - b)\} λk+1=max{0,λk+tk​(Axk−b)}

其中步长 $t_k=\frac{{\alpha }_k(z_{LR}({\lambda }^k)-z_{UB})}{\Vert A{x}^k-b{\Vert }^2}$，${\alpha }_k\in (0,2]$逐步递减。

4.6 Benders分解

Benders分解适用于含有"难"变量和"易"变量的问题，将大问题分解为主问题（Master Problem）和子问题（Subproblem）。

适用结构

$$\begin{array}{rl}\min & c^{T}x+d^{T}y\\ \text{s.t.}& Ax+By\le b\\ & x\in X,y\ge 0\end{array}$$

其中固定 $x$ 后，关于 $y$ 的子问题是LP。

算法流程

Benders分解：

1. 初始化：主问题（只有x变量）加上下界估计
2. 重复：
   a. 求解主问题，得到 x*
   b. 固定 x=x*，求解子问题（关于y的LP）
   c. 如果子问题不可行 → 添加可行性割（排除 x*）
   d. 如果子问题最优 → 添加最优性割（收紧下界）
   e. 直到上下界足够接近

4.7 经典建模技巧

逻辑约束建模

大M法是最常用的技巧，引入一个足够大的常数 $M$ 和二进制变量 δ ∈ { 0 , 1 } \delta \in \{0,1\} δ∈{0,1} 来建模逻辑条件：

逻辑条件	建模方式
如果-则：若 $\delta =1$ 则 $f(x)\le 0$	$f(x)\le M(1-\delta )$
或约束：$g_1(x)\le 0$ 或 $g_2(x)\le 0$	$g_1(x)\le M\delta ,g_2(x)\le M(1-\delta )$
选k个：$n$ 个约束中满足至少 $k$ 个	引入 $n$ 个二进制变量 ${\delta }_i$
SOS约束	特殊有序集

固定成本建模

$$\text{总成本}=\sum _j{f}_j{\delta }_j+\sum _j{c}_j{x}_j$$

其中 ${\delta }_j=1$ 当且仅当 $x_j>0$，需要约束 $x_j\le M_j{\delta }_j$。

分段线性函数

用SOS2（Special Ordered Set Type 2）或二进制变量建模分段线性凸/凹函数。

4.8 Python实现

使用PuLP

from pulp import *

# 背包问题示例
# 物品重量、价值
weights = [2, 3, 4, 5]
values = [3, 4, 5, 6]
capacity = 8
n = len(weights)

prob = LpProblem("Knapsack", LpMaximize)

# 二进制变量：是否选择物品i
x = [LpVariable(f"x{i}", cat='Binary') for i in range(n)]

# 目标：最大化总价值
prob += lpSum(values[i] * x[i] for i in range(n))

# 约束：总重量不超过容量
prob += lpSum(weights[i] * x[i] for i in range(n)) <= capacity

prob.solve()

print(f"选中物品: {[i for i in range(n) if value(x[i]) > 0.5]}")
print(f"最大价值: {value(prob.objective)}")

使用Google OR-Tools

from ortools.sat.python import cp_model

model = cp_model.CpModel()

# 选址问题：5个候选仓库，6个客户
num_warehouses = 5
num_customers = 6
fixed_cost = [100, 120, 150, 130, 110]  # 开设成本
# 运输成本矩阵（仓库→客户）
transport_cost = [
    [10, 12, 15, 20, 18, 25],
    [15, 10, 12, 18, 22, 20],
    [20, 18, 10, 15, 12, 18],
    [18, 15, 20, 10, 15, 12],
    [12, 20, 18, 15, 10, 15],
]
demand = [50, 40, 60, 30, 70, 45]
capacity_w = [100, 120, 150, 130, 110]

# 决策变量
open_w = [model.NewBoolVar(f'open_{w}') for w in range(num_warehouses)]
assign = [[model.NewBoolVar(f'assign_{w}_{c}') 
           for c in range(num_customers)] for w in range(num_warehouses)]

# 约束：每个客户恰好被一个仓库服务
for c in range(num_customers):
    model.Add(sum(assign[w][c] for w in range(num_warehouses)) == 1)

# 约束：只有开放的仓库才能服务客户
for w in range(num_warehouses):
    for c in range(num_customers):
        model.Add(assign[w][c] <= open_w[w])

# 约束：仓库容量限制
for w in range(num_warehouses):
    model.Add(sum(demand[c] * assign[w][c] for c in range(num_customers)) 
              <= capacity_w[w] * open_w[w])

# 目标：最小化总成本
total_cost = (sum(fixed_cost[w] * open_w[w] for w in range(num_warehouses)) +
              sum(transport_cost[w][c] * assign[w][c] 
                  for w in range(num_warehouses) 
                  for c in range(num_customers)))
model.Minimize(total_cost)

# 求解
solver = cp_model.CpSolver()
solver.parameters.max_time_in_seconds = 60
status = solver.Solve(model)

if status in (cp_model.OPTIMAL, cp_model.FEASIBLE):
    print(f"总成本: {solver.ObjectiveValue()}")
    for w in range(num_warehouses):
        if solver.Value(open_w[w]):
            clients = [c for c in range(num_customers) if solver.Value(assign[w][c])]
            print(f"  仓库{w} 开放，服务客户: {clients}")

---

第五篇：非线性规划（NLP）

5.1 非线性规划概述

当目标函数或约束中包含非线性函数时，问题变为非线性规划：

$$\begin{array}{rl}\min & f(x)\\ \text{s.t.}& g_i(x)\le 0,i=1,...,m\\ & h_j(x)=0,j=1,...,p\end{array}$$

非线性规划的挑战

1. 局部最优：非凸NLP可能有多个局部最优解，找到全局最优是NP-hard
2. 收敛速度：许多算法只能保证局部收敛
3. 初始化敏感：算法的性能高度依赖初始点
4. 约束处理复杂：需要专门的约束处理机制

5.2 无约束优化方法

最速下降法（梯度下降法）

最简单的迭代方法，沿负梯度方向搜索：

$$x^{k+1}=x^k-{\alpha }_k\nabla f(x^k)$$

其中 ${\alpha }_k$ 是步长，通过线搜索确定。

线搜索方法：

• 精确线搜索：${\alpha }_k=\arg {\min }_{\alpha >0}f(x^k-\alpha \nabla f(x^k))$
• Armijo回溯线搜索：满足 $f(x^k+\alpha d^k)\le f(x^k)+c_1\alpha \nabla f(x^k{)}^T{d}^k$
• Wolfe条件：Armijo条件 + 曲率条件

收敛性：

• 收敛速度：线性（$O((1-\frac{m}{M}{)}^k)$，其中 $m,M$ 是Hessian的最小和最大特征值）
• 对条件数敏感：当 $\kappa =M/m$ 很大时收敛很慢
• 优点：简单、每步计算便宜
• 缺点：锯齿形路径，收敛慢

牛顿法

利用二阶信息（Hessian矩阵）加速收敛：

$$x^{k+1}=x^k-[{\nabla }^{2}f(x^k){]}^{-1}\nabla f(x^k)$$

收敛性：

• 局部二次收敛：$\Vert x^{k+1}-x^{\ast }\Vert \le C\Vert x^k-x^{\ast }{\Vert }^2$
• 缺点：需要计算和存储Hessian矩阵，$O(n^3)$ 求逆
• 可能Hessian不正定

拟牛顿法

用近似矩阵 $B_k\approx {\nabla }^{2}f(x^k)$ 替代Hessian，避免二阶导数计算：

BFGS公式（最流行的拟牛顿法）：
$$B_{k+1}=B_k+\frac{y_k{y}_k^T}{y_k^T{s}_k}-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}$$

其中 $s_k=x^{k+1}-x^k$，$y_k=\nabla f(x^{k+1})-\nabla f(x^k)$。

L-BFGS（Limited-memory BFGS）：只保存最近 $m$ 步的 { s k , y k } \{s_k, y_k\} {sk​,yk​} 信息，适合大规模问题。

共轭梯度法

不需要存储Hessian或其近似，适合大规模问题：

$$d^k=-\nabla f(x^k)+{\beta }_k{d}^{k-1}$$

不同 ${\beta }_k$ 的选择：

• Fletcher-Reeves：${\beta }_k=\frac{\Vert \nabla f(x^k){\Vert }^2}{\Vert \nabla f(x^{k-1}){\Vert }^2}$
• Polak-Ribière：${\beta }_k=\frac{\nabla f(x^k{)}^T(\nabla f(x^k)-\nabla f(x^{k-1}))}{\Vert \nabla f(x^{k-1}){\Vert }^2}$

方法对比

方法	收敛速度	每步计算量	存储需求	适用规模
最速下降	线性	$O(n)$	$O(n)$	大规模
牛顿法	超线性/二次	$O(n^3)$	$O(n^2)$	小规模
BFGS	超线性	$O(n^2)$	$O(n^2)$	中等规模
L-BFGS	超线性	$O(mn)$	$O(mn)$	大规模
共轭梯度	线性/超线性	$O(n)$	$O(n)$	大规模

5.3 约束优化方法

序列二次规划（SQP）

SQP是求解约束NLP最有效的方法之一。在每一步迭代中，求解一个二次规划子问题：

$$\begin{array}{rl}\underset{d}{\min }& \nabla f(x^k{)}^{T}d+\frac{1}{2}{d}^T{H}_{k}d\\ \text{s.t.}& g_i(x^k)+\nabla g_i(x^k{)}^{T}d\le 0\\ & h_j(x^k)+\nabla h_j(x^k{)}^{T}d=0\end{array}$$

其中 $H_k$ 是拉格朗日函数Hessian的BFGS近似。

内点法（Barrier Method）

将约束通过障碍函数纳入目标：

$$\min f(x)-\mu \sum _i\ln (-g_i(x))-\mu \sum _j(h_j(x){)}^2$$

增广拉格朗日法（ALM）

结合拉格朗日乘子法和罚函数法：

$$L_A(x,\lambda ,\rho )=f(x)+\sum _i{\lambda }_i{g}_i(x)+\frac{\rho }{2}\sum _i[\max (0,g_i(x)){]}^2$$

罚函数法

将约束违反度加入目标函数：

$$\min P(x,\rho )=f(x)+\rho \left[\sum _i[\max (0,g_i(x)){]}^2+\sum _j{h}_j(x{)}^2\right]$$

随着 $\rho \to \infty$，罚函数的最优解趋向原问题最优解。

5.4 二次规划（QP）

二次规划是一类特殊的非线性规划，目标函数为二次，约束为线性：

$$\begin{array}{rl}\min & \frac{1}{2}{x}^{T}Qx+c^{T}x\\ \text{s.t.}& Ax\le b\\ & A_{eq}x=b_{eq}\end{array}$$

凸QP（$Q⪰0$）有全局最优解，可以用有效算法求解。

求解方法：

• 积极集法：维护有效约束集合，迭代调整
• 内点法：适合大规模问题
• 对偶法：转化为对偶问题

应用：投资组合优化（Markowitz模型）、SVM训练、模型预测控制。

5.5 Python实现

使用SciPy

import numpy as np
from scipy.optimize import minimize

# 无约束优化：Rosenbrock函数
def rosenbrock(x):
    return 100*(x[1]-x[0]**2)**2 + (1-x[0])**2

x0 = [-1.0, 1.0]
result = minimize(rosenbrock, x0, method='BFGS')
print(f"最优解: {result.x}")
print(f"最优值: {result.fun}")

from scipy.optimize import minimize

# 约束优化
# min (x₁ - 1)² + (x₂ - 2.5)²
# s.t. x₁ - 2x₂ + 2 ≥ 0
#      -x₁ - 2x₂ + 6 ≥ 0
#      -x₁ + 2x₂ + 2 ≥ 0
#      0 ≤ x₁ ≤ 5, 0 ≤ x₂ ≤ 5

def objective(x):
    return (x[0]-1)**2 + (x[1]-2.5)**2

constraints = [
    {'type': 'ineq', 'fun': lambda x: x[0] - 2*x[1] + 2},
    {'type': 'ineq', 'fun': lambda x: -x[0] - 2*x[1] + 6},
    {'type': 'ineq', 'fun': lambda x: -x[0] + 2*x[1] + 2},
]
bounds = [(0, 5), (0, 5)]

result = minimize(objective, [0, 0], method='SLSQP',
                  bounds=bounds, constraints=constraints)
print(f"最优解: x₁={result.x[0]:.4f}, x₂={result.x[1]:.4f}")
print(f"最优值: {result.fun:.4f}")

使用CVXPY（凸优化建模）

import cvxpy as cp
import numpy as np

# 二次规划：投资组合优化
n = 5  # 5只股票
np.random.seed(42)
returns = np.random.randn(100, n)  # 历史收益率
mu = returns.mean(axis=0)          # 期望收益率
Sigma = np.cov(returns.T)          # 协方差矩阵

# 决策变量：投资比例
w = cp.Variable(n)

# 目标：最小化风险（方差）
risk = cp.quad_form(w, Sigma)

# 约束
constraints = [
    cp.sum(w) == 1,          # 比例之和为1
    mu @ w >= 0.05,          # 最低收益5%
    w >= 0,                  # 不做空
]

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

print(f"最优投资比例: {np.round(w.value, 4)}")
print(f"最小风险(方差): {prob.value:.6f}")
print(f"期望收益: {mu @ w.value:.4f}")

---

第六篇：动态规划（DP）

6.1 动态规划基本思想

动态规划（Dynamic Programming, DP）由Richard Bellman于1950年代提出，用于求解具有多阶段决策过程特征的优化问题。

核心原理

Bellman最优性原理：

一个最优策略具有这样的性质：不论初始状态和初始决策如何，余下的决策序列对于由初始决策导致的状态而言，必须构成最优策略。

用Bellman方程表示：
V ( s ) = max ⁡ a ∈ A ( s ) { R ( s , a ) + γ ∑ s ′ P ( s ′ ∣ s , a ) V ( s ′ ) } V(s) = \max_{a \in A(s)} \{R(s,a) + \gamma \sum_{s'} P(s'|s,a) V(s')\} V(s)=a∈A(s)max​{R(s,a)+γs′∑​P(s′∣s,a)V(s′)}

（这是随机DP的形式，确定性DP中 $P(s^{\prime }|s,a)$ 是确定性的）

适用条件

1. 最优子结构：问题的最优解包含子问题的最优解
2. 重叠子问题：不同的子问题可能共享更小的子问题
3. 无后效性：未来决策只依赖当前状态，与到达该状态的路径无关

6.2 最优子结构与重叠子问题

最优子结构示例

最短路径问题：从A到D的最短路径经过B，则A→B的子路径也是A到B的最短路径。

重叠子问题示例

斐波那契数列：

f(5) = f(4) + f(3)
f(4) = f(3) + f(2)    ← f(3)被重复计算
f(3) = f(2) + f(1)    ← f(2)被重复计算

DP通过记忆化（自顶向下）或制表法（自底向上）避免重复计算。

6.3 状态转移方程

设计DP算法的关键是定义：

1. 状态：描述问题子阶段的变量
2. 决策：在每个状态下的可选操作
3. 状态转移方程：状态之间的递推关系
4. 边界条件：最小子问题的解

自顶向下 vs 自底向上

# 自顶向下（记忆化递归）
def fib_memo(n, memo={}):
    if n in memo: return memo[n]
    if n <= 1: return n
    memo[n] = fib_memo(n-1, memo) + fib_memo(n-2, memo)
    return memo[n]

# 自底向上（制表法）
def fib_tab(n):
    dp = [0] * (n + 1)
    dp[1] = 1
    for i in range(2, n + 1):
        dp[i] = dp[i-1] + dp[i-2]
    return dp[n]

6.4 经典问题

0-1背包问题

问题：$n$ 个物品，每个物品有重量 $w_i$ 和价值 $v_i$，背包容量 $W$，最大化总价值。

状态定义：$dp[i][w]$ = 前 $i$ 个物品在容量为 $w$ 时的最大价值

状态转移方程：
d p [ i ] [ w ] = max ⁡ { d p [ i − 1 ] [ w ] ,    d p [ i − 1 ] [ w − w i ] + v i } dp[i][w] = \max\{dp[i-1][w], \; dp[i-1][w-w_i] + v_i\} dp[i][w]=max{dp[i−1][w],dp[i−1][w−wi​]+vi​}

（不选第 $i$ 个物品 vs 选第 $i$ 个物品）

Python实现：

def knapsack_01(weights, values, capacity):
    n = len(weights)
    dp = [[0] * (capacity + 1) for _ in range(n + 1)]
    
    for i in range(1, n + 1):
        for w in range(capacity + 1):
            dp[i][w] = dp[i-1][w]  # 不选物品i
            if w >= weights[i-1]:
                dp[i][w] = max(dp[i][w], 
                               dp[i-1][w-weights[i-1]] + values[i-1])
    
    # 回溯找出选择了哪些物品
    selected = []
    w = capacity
    for i in range(n, 0, -1):
        if dp[i][w] != dp[i-1][w]:
            selected.append(i-1)
            w -= weights[i-1]
    
    return dp[n][capacity], selected

weights = [2, 3, 4, 5]
values = [3, 4, 5, 6]
capacity = 8
max_val, items = knapsack_01(weights, values, capacity)
print(f"最大价值: {max_val}, 选中物品: {items}")

最长公共子序列（LCS）

$$dp[i][j]=\{\begin{array}{ll}dp[i-1][j-1]+1& \text{if }{s}_1[i]=s_2[j]\\ \max (dp[i-1][j],dp[i][j-1])& \text{otherwise}\end{array}$$

最长递增子序列（LIS）

d p [ i ] = max ⁡ j < i , a j < a i { d p [ j ] } + 1 dp[i] = \max_{j<i, a_j<a_i}\{dp[j]\} + 1 dp[i]=j<i,aj​<ai​max​{dp[j]}+1

编辑距离

$$dp[i][j]=\min \{\begin{array}{ll}dp[i-1][j]+1& \text{（删除）}\\ dp[i][j-1]+1& \text{（插入）}\\ dp[i-1][j-1]+[s_1[i]\ne s_2[j]]& \text{（替换）}\end{array}$$

矩阵链乘法

$n$ 个矩阵 $A_1,A_2,...,A_n$（$A_i$ 的维度为 $p_{i-1}\times p_i$），求乘法顺序使标量乘法次数最少。

d p [ i ] [ j ] = min ⁡ i ≤ k < j { d p [ i ] [ k ] + d p [ k + 1 ] [ j ] + p i − 1 p k p j } dp[i][j] = \min_{i \leq k < j}\{dp[i][k] + dp[k+1][j] + p_{i-1} p_k p_j\} dp[i][j]=i≤k<jmin​{dp[i][k]+dp[k+1][j]+pi−1​pk​pj​}

经典DP问题分类

类别	问题	状态定义
线性DP	最长递增子序列、编辑距离	$dp[i]$ 或 $dp[i][j]$
区间DP	矩阵链乘法、石子合并	$dp[i][j]$ = 区间 $[i,j]$ 的最优值
背包DP	0-1背包、完全背包	$dp[i][w]$
树形DP	树的最大独立集	$dp[u][0/1]$
状态压缩DP	TSP、集合覆盖	$dp[S][i]$ = 已访问集合 $S$ 且在 $i$
数位DP	数字中不含某模式的计数	$dp[pos][state][...]$
概率DP	马尔可夫决策过程	$dp[i]$ = 期望值

6.5 近似动态规划

当状态空间太大时，精确DP不可行（维度灾难），需要近似方法：

函数近似

用参数化函数 $V(s;\theta )\approx V^{\ast }(s)$ 替代精确值函数：

• 线性函数近似：$V(s;\theta )={\theta }^T\varphi (s)$
• 神经网络近似：$V(s;\theta )=\text{NN}(s;\theta )$

Rollout算法

从当前状态向前模拟（rollout）若干步，用启发式策略评估，选择最优行动。

近似动态规划（ADP）/ 强化学习

ADP是动态规划与机器学习的交叉领域，核心思想是通过采样和函数近似来克服维度灾难。这与强化学习（尤其是Q-learning和策略梯度方法）有密切联系。

---

第七篇：图论与网络优化

7.1 图论基础

图 $G=(V,E)$ 由顶点集 $V$ 和边集 $E$ 组成。

图的分类

分类维度	类型	说明
边方向	无向图 / 有向图	边是否带方向
边权重	无权图 / 加权图	边是否带数值
连通性	连通图 / 非连通图	任意两点间是否可达
特殊结构	树、二分图、完全图	特定拓扑结构

图的存储

• 邻接矩阵：$A[i][j]$ 表示 $i$ 到 $j$ 的边权，$O(n^2)$ 空间
• 邻接表：每个顶点维护邻居列表，$O(n+m)$ 空间
• 边列表：直接存储所有边，适合稀疏图

7.2 最短路径问题

Dijkstra算法

适用于非负权重的单源最短路径问题。

时间复杂度：$O((n+m)\log n)$（使用优先队列）

import heapq

def dijkstra(graph, source):
    """
    graph: 邻接表 {u: [(v, w), ...]}
    返回从source到所有点的最短距离
    """
    dist = {v: float('inf') for v in graph}
    dist[source] = 0
    pq = [(0, source)]  # (距离, 顶点)
    
    while pq:
        d, u = heapq.heappop(pq)
        if d > dist[u]:
            continue  # 跳过已处理的旧条目
        for v, w in graph[u]:
            if dist[u] + w < dist[v]:
                dist[v] = dist[u] + w
                heapq.heappush(pq, (dist[v], v))
    
    return dist

Bellman-Ford算法

适用于含负权边的单源最短路径，可检测负权环。

时间复杂度：$O(nm)$

def bellman_ford(edges, n, source):
    """
    edges: [(u, v, w), ...]
    n: 顶点数
    """
    dist = [float('inf')] * n
    dist[source] = 0
    
    for _ in range(n - 1):
        for u, v, w in edges:
            if dist[u] + w < dist[v]:
                dist[v] = dist[u] + w
    
    # 检测负权环
    for u, v, w in edges:
        if dist[u] + w < dist[v]:
            return None  # 存在负权环
    
    return dist

Floyd-Warshall算法

全源最短路径，允许负权边（无负权环）。

时间复杂度：$O(n^3)$

def floyd_warshall(n, dist):
    """dist: n×n 邻接矩阵，dist[i][j]为i到j的直接距离"""
    for k in range(n):
        for i in range(n):
            for j in range(n):
                dist[i][j] = min(dist[i][j], dist[i][k] + dist[k][j])
    return dist

A*算法

带启发式的最短路径搜索，常用于游戏和导航：

$$f(n)=g(n)+h(n)$$

• $g(n)$：从起点到 $n$ 的实际距离
• $h(n)$：从 $n$ 到终点的启发式估计（如欧几里得距离、曼哈顿距离）
• 要求 $h(n)$ 是可容许的（不高估实际距离）

7.3 最小生成树

Kruskal算法

按边权从小到大排序，逐条加入，如果形成环则跳过。使用**并查集（Union-Find）**高效判断环。

时间复杂度：$O(m\log m)$

Prim算法

从一个顶点出发，每次加入连接已选顶点集和未选顶点集的最小权边。

时间复杂度：$O(m\log n)$（使用优先队列）

7.4 网络流问题

最大流问题

问题：给定有向图，每条边有容量上限，求从源点 $s$ 到汇点 $t$ 的最大流量。

核心定理：最大流 = 最小割（Max-Flow Min-Cut Theorem）

Ford-Fulkerson方法

重复：
1. 在残余图中找到s→t的增广路径
2. 计算瓶颈容量 = 路径上最小残余容量
3. 沿增广路径推送流量，更新残余图
直到没有增广路径

Edmonds-Karp算法

Ford-Fulkerson的BFS实现，总复杂度 $O(n{m}^2)$。

Dinic算法

使用层次图和阻塞流，复杂度 $O(n^{2}m)$，实践中非常高效。

最小费用最大流

在满足最大流的前提下，使总运输费用最小。常用连续最短路算法（Successive Shortest Path）或圈消除算法。

7.5 匹配与指派问题

二分图最大匹配

匈牙利算法：$O(n^3)$ 时间求解二分图最大匹配。

指派问题

$n$ 个工人分配 $n$ 个任务，每人一个任务，最小化总成本。这是特殊的LP和特殊的网络流问题。

min ⁡ ∑ i ∑ j c i j x i j s.t. ∑ j x i j = 1 , ∀ i ∑ i x i j = 1 , ∀ j x i j ∈ { 0 , 1 } \begin{aligned} \min \quad & \sum_i \sum_j c_{ij} x_{ij} \\ \text{s.t.} \quad & \sum_j x_{ij} = 1, \quad \forall i \\ & \sum_i x_{ij} = 1, \quad \forall j \\ & x_{ij} \in \{0, 1\} \end{aligned} mins.t.​i∑​j∑​cij​xij​j∑​xij​=1,∀ii∑​xij​=1,∀jxij​∈{0,1}​

匈牙利算法可在 $O(n^3)$ 时间精确求解。

7.6 旅行商问题（TSP）

问题定义

旅行商要访问 $n$ 个城市各一次并返回出发城市，求最短总路程。这是经典的NP-hard组合优化问题。

min ⁡ ∑ i , j c i j x i j s.t. ∑ j x i j = 1 , ∀ i （每个城市恰好一条出边） ∑ i x i j = 1 , ∀ j （每个城市恰好一条入边） u i − u j + n x i j ≤ n − 1 , ∀ i ≠ j , i , j ≥ 2 （子回路消除） x i j ∈ { 0 , 1 } \begin{aligned} \min \quad & \sum_{i,j} c_{ij} x_{ij} \\ \text{s.t.} \quad & \sum_j x_{ij} = 1, \quad \forall i \quad \text{（每个城市恰好一条出边）} \\ & \sum_i x_{ij} = 1, \quad \forall j \quad \text{（每个城市恰好一条入边）} \\ & u_i - u_j + n x_{ij} \leq n-1, \quad \forall i \neq j, i,j \geq 2 \quad \text{（子回路消除）} \\ & x_{ij} \in \{0, 1\} \end{aligned} mins.t.​i,j∑​cij​xij​j∑​xij​=1,∀i（每个城市恰好一条出边）i∑​xij​=1,∀j（每个城市恰好一条入边）ui​−uj​+nxij​≤n−1,∀i=j,i,j≥2（子回路消除）xij​∈{0,1}​

求解方法

方法	类型	适用规模
MIP精确求解	精确	≤ 50 城市
动态规划（Held-Karp）	精确	≤ 25 城市
最近邻启发式	构造式启发式	大规模
Christofides算法	近似算法（比值1.5）	大规模
Lin-Kernighan	局部搜索	大规模
遗传算法	元启发式	大规模
蚁群算法	元启发式	大规模
模拟退火	元启发式	大规模

Python实现（MIP求解）

from ortools.sat.python import cp_model
import numpy as np

def solve_tsp(distance_matrix):
    n = len(distance_matrix)
    model = cp_model.CpModel()
    
    # 边变量
    x = {}
    for i in range(n):
        for j in range(n):
            if i != j:
                x[i,j] = model.NewBoolVar(f'x_{i}_{j}')
    
    # 度约束：每个节点恰好一进一出
    for i in range(n):
        model.Add(sum(x[i,j] for j in range(n) if j != i) == 1)
        model.Add(sum(x[j,i] for j in range(n) if j != i) == 1)
    
    # MTZ子回路消除
    u = [model.NewIntVar(0, n-1, f'u_{i}') for i in range(n)]
    model.Add(u[0] == 0)
    for i in range(1, n):
        for j in range(1, n):
            if i != j:
                model.Add(u[i] - u[j] + n * x[i,j] <= n - 1)
    
    # 目标
    model.Minimize(sum(distance_matrix[i][j] * x[i,j] 
                       for i in range(n) for j in range(n) if i != j))
    
    solver = cp_model.CpSolver()
    solver.parameters.max_time_in_seconds = 30
    solver.Solve(model)
    
    # 提取路径
    path = [0]
    current = 0
    for _ in range(n-1):
        for j in range(n):
            if j != current and solver.Value(x[current,j]):
                path.append(j)
                current = j
                break
    path.append(0)
    
    return path, solver.ObjectiveValue()

7.7 车辆路径问题（VRP）

问题定义

VRP是TSP的推广：有 $K$ 辆车从配送中心出发服务 $n$ 个客户，每辆车有容量限制，需要决定每辆车的路线使总成本最小。

VRP变体

变体	说明
CVRP	带容量约束的VRP
VRPTW	带时间窗的VRP
VRPPD	带取送货的VRP
MDVRP	多配送中心VRP
DVRP	动态VRP（需求实时到来）

求解框架

VRP求解方法：
├── 精确方法
│   ├── 分支定界/分支切割
│   ├── 分支定价（Branch and Price）
│   └── 集合覆盖/集合划分
├── 构造式启发式
│   ├── 节约算法（Clarke-Wright）
│   ├── 最近邻
│   └── 扫描算法
├── 改进式启发式
│   ├── 2-opt, 3-opt
│   ├── Or-opt
│   └── Cross exchange
└── 元启发式
    ├── ALNS（自适应大邻域搜索）⭐ 最有效
    ├── 遗传算法
    ├── 模拟退火
    └── 禁忌搜索

---

第八篇：启发式与元启发式算法

8.1 启发式算法概述

对于NP-hard问题，精确算法在问题规模大时计算时间不可接受。启发式算法牺牲最优性保证，换取合理的计算时间，通常能找到质量可接受的近优解。

分类

启发式算法
├── 构造式启发式（Constructive）
│   ├── 贪心算法
│   ├── 最近邻
│   └── 节约算法
├── 改进式启发式（Improvement）
│   ├── 局部搜索
│   ├── k-opt
│   └── 邻域搜索
└── 元启发式（Metaheuristic）  ← 框架级别
    ├── 单解元启发式
    │   ├── 模拟退火
    │   ├── 禁忌搜索
    │   ├── 变邻域搜索
    │   └── 自适应大邻域搜索
    └── 种群元启发式
        ├── 遗传算法
        ├── 粒子群优化
        ├── 蚁群算法
        └── 差分进化

8.2 贪心算法

核心思想：每一步做出当前看起来最好的选择，不考虑未来影响。

应用场景

• Huffman编码：频率低的字符用长编码
• 活动选择问题：选最早结束的活动
• Kruskal算法：每次选最小权边

局限性

贪心不一定得到最优解，但对某些问题（拟阵结构）可以证明最优性。

8.3 局部搜索

基本框架

局部搜索：
1. 生成初始解 s
2. 重复：
   a. 在邻域 N(s) 中找到更好的解 s'
   b. 如果存在 s' → s = s'
   c. 否则 → 停止（到达局部最优）
3. 返回 s

邻域结构

问题	常用邻域
TSP	2-opt（交换两条边）、3-opt、Or-opt
调度	交换两个任务的顺序、移动任务
背包	翻转一个物品的选择
VRP	relocate（移动客户）、swap（交换客户）、2-opt*

局部最优的困境

目标值
  ↑
  │    ╱╲
  │   ╱  ╲      ╱╲      ← 全局最优
  │  ╱    ╲    ╱  ╲
  │ ╱  局部 ╲  ╱    ╲
  │╱   最优  ╲╱      ╲
  └────────────────────→ 解空间
       ↑        ↑
    局部最优  全局最优

元启发式算法的核心就是如何跳出局部最优。

8.4 模拟退火（SA）

核心思想

借鉴金属退火过程：高温时原子自由运动（接受差解），逐步降温使系统达到最低能量状态（最优解）。

算法

模拟退火算法：

1. 初始化：初始解 s，初始温度 T₀，终止温度 T_min，降温系数 α
2. 当 T > T_min：
   a. 重复 L 次（每温度的迭代次数）：
      - 在邻域随机产生新解 s'
      - 计算 ΔE = f(s') - f(s)
      - 如果 ΔE < 0 → 接受 s'（改进解）
      - 否则以概率 exp(-ΔE/T) 接受 s'（Metropolis准则）
   b. 降温：T ← αT（通常 α ∈ [0.95, 0.99]）
3. 返回找到的最好解

接受概率：
$$P(\text{accept})=\{\begin{array}{ll}1& \text{if }\Delta E\le 0\\ e^{-\Delta E/T}& \text{if }\Delta E>0\end{array}$$

关键参数

参数	影响	调参建议
$T_0$	初始接受率	初始接受率约80%
$\alpha$	降温速度	0.95-0.99
$L$	每温度搜索深度	与问题规模成正比
$T_{min}$	终止条件	足够小或迭代次数

Python实现

import math
import random

def simulated_annealing(problem, initial_solution, 
                         T0=100, T_min=0.01, alpha=0.99, L=100):
    """
    problem: 问题对象，需提供 evaluate(sol) 和 neighbor(sol) 方法
    """
    current = initial_solution
    best = current
    best_cost = problem.evaluate(current)
    T = T0
    
    while T > T_min:
        for _ in range(L):
            # 生成邻域解
            new_sol = problem.neighbor(current)
            new_cost = problem.evaluate(new_sol)
            current_cost = problem.evaluate(current)
            
            delta = new_cost - current_cost
            
            # Metropolis准则
            if delta < 0 or random.random() < math.exp(-delta / T):
                current = new_sol
                
                # 更新最优解
                if new_cost < best_cost:
                    best = new_sol
                    best_cost = new_cost
        
        T *= alpha
    
    return best, best_cost

# TSP示例
class TSPProblem:
    def __init__(self, cities):
        self.cities = cities
        self.n = len(cities)
    
    def evaluate(self, route):
        """计算路线总长度"""
        total = 0
        for i in range(len(route)):
            c1 = self.cities[route[i]]
            c2 = self.cities[route[(i+1) % self.n]]
            total += math.sqrt((c1[0]-c2[0])**2 + (c1[1]-c2[1])**2)
        return total
    
    def neighbor(self, route):
        """2-opt邻域"""
        new_route = route[:]
        i, j = sorted(random.sample(range(self.n), 2))
        new_route[i:j+1] = reversed(new_route[i:j+1])
        return new_route

8.5 遗传算法（GA）

核心思想

借鉴自然选择和遗传机制：种群中的个体（解）通过选择、交叉、变异进化，适应度高的个体更容易生存和繁殖。

算法流程

遗传算法：

1. 初始化：随机生成初始种群 P₀（大小为N）
2. 对每一代 g = 1, 2, ...:
   a. 适应度评估：计算每个个体的适应度 f(sᵢ)
   b. 选择（Selection）：按适应度选择父代
      - 轮盘赌选择（适应度比例选择）
      - 锦标赛选择
      - 精英保留
   c. 交叉（Crossover）：父代配对，交换基因生成子代
      - 单点交叉、两点交叉、均匀交叉
      - 交叉概率 p_c（通常0.6-0.9）
   d. 变异（Mutation）：以小概率改变个体的某些基因
      - 变异概率 p_m（通常0.01-0.1）
   e. 更新种群：P_{g+1} = 选择(P_g ∪ 子代)
3. 返回种群中适应度最高的个体

编码方式

问题	编码	交叉	变异
TSP	排列编码	PMX/OX/CX	交换/逆序/插入
背包	二进制串	均匀交叉	位翻转
调度	优先级/排列	顺序交叉	交换/移动
连续优化	实数向量	SBX/BLX	高斯扰动

TSP专用交叉算子

**顺序交叉（OX）**示例：

父代1:  1 2 | 3 4 5 | 6 7 8
父代2:  3 7 | 5 1 6 | 8 2 4

子代1:  5 1 | 3 4 5 | 6 8 2  → 去重后调整 → 5 1 6 3 4 5 7 8 2
         (保留父代1片段，按父代2顺序填入)

Python实现

import random
import numpy as np

def genetic_algorithm(problem, pop_size=100, generations=500,
                      crossover_rate=0.8, mutation_rate=0.1, 
                      elitism_count=2):
    """
    problem: 需提供 create_individual(), fitness(ind), 
             crossover(p1, p2), mutate(ind) 方法
    """
    # 初始化种群
    population = [problem.create_individual() for _ in range(pop_size)]
    
    for gen in range(generations):
        # 评估适应度
        fitness_scores = [(ind, problem.fitness(ind)) for ind in population]
        fitness_scores.sort(key=lambda x: x[1], reverse=True)
        
        # 精英保留
        new_population = [fs[0] for fs in fitness_scores[:elitism_count]]
        
        # 生成新个体
        while len(new_population) < pop_size:
            # 锦标赛选择
            parent1 = tournament_select(fitness_scores, k=3)
            parent2 = tournament_select(fitness_scores, k=3)
            
            # 交叉
            if random.random() < crossover_rate:
                child1, child2 = problem.crossover(parent1, parent2)
            else:
                child1, child2 = parent1[:], parent2[:]
            
            # 变异
            if random.random() < mutation_rate:
                child1 = problem.mutate(child1)
            if random.random() < mutation_rate:
                child2 = problem.mutate(child2)
            
            new_population.extend([child1, child2])
        
        population = new_population[:pop_size]
        
        if gen % 50 == 0:
            best_fitness = fitness_scores[0][1]
            print(f"代 {gen}: 最佳适应度 = {best_fitness:.4f}")
    
    return fitness_scores[0]

def tournament_select(fitness_scores, k=3):
    """锦标赛选择"""
    candidates = random.sample(fitness_scores, k)
    return max(candidates, key=lambda x: x[1])[0]

8.6 粒子群优化（PSO）

核心思想

模拟鸟群觅食行为。每个粒子代表一个解，具有位置和速度，通过跟踪个体最优和全局最优来更新运动方向。

算法

PSO算法：

初始化：N个粒子，随机位置和速度
重复直到满足终止条件：
  对每个粒子 i：
    1. 评估当前位置的适应度
    2. 更新个体最优 pbest_i
    3. 更新全局最优 gbest
    4. 更新速度：
       v_i = ω·v_i + c₁·r₁·(pbest_i - x_i) + c₂·r₂·(gbest - x_i)
    5. 更新位置：
       x_i = x_i + v_i

参数说明：
  ω  = 惯性权重（平衡全局和局部搜索）
  c₁ = 认知系数（向个体最优学习的程度）
  c₂ = 社会系数（向全局最优学习的程度）
  r₁, r₂ = [0,1]随机数

速度更新的三个分量：

1. 惯性项 $\omega v_i$：保持当前运动趋势
2. 认知项 $c_1{r}_1(p_{best}-x_i)$：向自身最优位置移动
3. 社会项 $c_2{r}_2(g_{best}-x_i)$：向全局最优位置移动

Python实现

import numpy as np

class PSO:
    def __init__(self, objective, n_dim, n_particles=50,
                 w=0.7, c1=1.5, c2=1.5, bounds=None):
        self.objective = objective
        self.n_dim = n_dim
        self.n_particles = n_particles
        self.w = w
        self.c1 = c1
        self.c2 = c2
        self.bounds = bounds  # [(low, high), ...] for each dim
        
        # 初始化粒子
        if bounds:
            low = np.array([b[0] for b in bounds])
            high = np.array([b[1] for b in bounds])
            self.positions = np.random.uniform(low, high, (n_particles, n_dim))
            self.velocities = np.random.uniform(-(high-low)/2, 
                                                 (high-low)/2, 
                                                 (n_particles, n_dim))
        else:
            self.positions = np.random.randn(n_particles, n_dim)
            self.velocities = np.random.randn(n_particles, n_dim)
        
        self.pbest_positions = self.positions.copy()
        self.pbest_scores = np.full(n_particles, np.inf)
        self.gbest_position = None
        self.gbest_score = np.inf
    
    def optimize(self, max_iter=1000):
        for iteration in range(max_iter):
            # 评估
            scores = np.array([self.objective(p) for p in self.positions])
            
            # 更新个体最优
            better = scores < self.pbest_scores
            self.pbest_scores[better] = scores[better]
            self.pbest_positions[better] = self.positions[better]
            
            # 更新全局最优
            best_idx = np.argmin(scores)
            if scores[best_idx] < self.gbest_score:
                self.gbest_score = scores[best_idx]
                self.gbest_position = self.positions[best_idx].copy()
            
            # 更新速度和位置
            r1 = np.random.random((self.n_particles, self.n_dim))
            r2 = np.random.random((self.n_particles, self.n_dim))
            
            self.velocities = (self.w * self.velocities +
                              self.c1 * r1 * (self.pbest_positions - self.positions) +
                              self.c2 * r2 * (self.gbest_position - self.positions))
            self.positions += self.velocities
            
            # 边界处理
            if self.bounds:
                for d in range(self.n_dim):
                    self.positions[:, d] = np.clip(self.positions[:, d],
                                                   self.bounds[d][0],
                                                   self.bounds[d][1])
            
            if iteration % 100 == 0:
                print(f"Iter {iteration}: Best = {self.gbest_score:.6f}")
        
        return self.gbest_position, self.gbest_score

# 示例：求Rosenbrock函数最小值
pso = PSO(lambda x: 100*(x[1]-x[0]**2)**2 + (1-x[0])**2,
          n_dim=2, bounds=[(-5,5), (-5,5)])
best_pos, best_val = pso.optimize(max_iter=1000)
print(f"最优解: {best_pos}, 最优值: {best_val}")

8.7 蚁群算法（ACO）

核心思想

模拟蚂蚁觅食行为。蚂蚁在路径上释放信息素（pheromone），短路径上信息素积累更快，吸引更多蚂蚁，形成正反馈。

算法（以TSP为例）

蚁群算法：

1. 初始化：所有边上信息素 τᵢⱼ = τ₀
2. 重复直到收敛：
   a. 对每只蚂蚁 k：
      - 从随机城市出发
      - 重复选择下一个城市，概率为：
        p_{ij}^k = [τᵢⱼ]^α · [ηᵢⱼ]^β / Σ [τᵢₗ]^α · [ηᵢₗ]^β
        其中 ηᵢⱼ = 1/dᵢⱼ（启发式信息），α和β是参数
      - 完成一条回路，记录总长度
   b. 信息素更新：
      - 蒸发：τᵢⱼ ← (1-ρ)τᵢⱼ
      - 增强：τᵢⱼ += Σ Δτᵢⱼ^k
        Δτᵢⱼ^k = Q/L_k（L_k是蚂蚁k的路径长度）

参数：

参数	含义	建议值
$\alpha$	信息素重要程度	1-5
$\beta$	启发信息重要程度	2-5
$\rho$	信息素蒸发率	0.1-0.5
$Q$	信息素强度常数	1

8.8 禁忌搜索（TS）

核心思想

在局部搜索的基础上引入禁忌表（Tabu List），记录最近做过的移动，禁止在一定步数内做逆向移动，从而跳出局部最优。

算法

禁忌搜索：

1. 初始解 s，禁忌表 T = ∅
2. 重复：
   a. 生成邻域 N(s)
   b. 从中选择不在禁忌表中的最优解 s'
      （除非s'比当前最优解更好——渴望准则/Aspiration）
   c. s ← s'
   d. 更新禁忌表（FIFO策略）
3. 直到满足终止条件

关键设计

• 禁忌表长度：通常 5-20，太短容易循环，太长限制搜索
• 渴望准则：如果候选解比历史最优更好，即使被禁忌也接受
• 多样化策略：当搜索停滞时，跳到新的搜索区域

8.9 差分进化（DE）

核心思想

差分进化是一种简单但强大的实数优化算法，通过种群中个体之间的差分向量来产生变异，再通过交叉和选择来进化。

算法（DE/rand/1/bin）

差分进化：

初始化种群 {x₁, x₂, ..., xN}（随机生成）
重复：
  对每个个体 xᵢ：
    1. 变异：随机选择三个不同个体 xᵣ₁, xᵣ₂, xᵣ₃
       v = xᵣ₁ + F · (xᵣ₂ - xᵣ₃)  （F ∈ [0,2]为缩放因子）
    2. 交叉：
       uⱼ = vⱼ  if rand() < CR or j = j_rand
       uⱼ = xᵢⱼ otherwise  （CR ∈ [0,1]为交叉概率）
    3. 选择：
       如果 f(u) ≤ f(xᵢ) → xᵢ ← u
       否则 → 保持 xᵢ 不变
直到满足终止条件

8.10 自适应大邻域搜索（ALNS）

ALNS是求解VRP、调度等复杂组合优化问题的最有效框架之一。

核心思想

维护一组**破坏（Destroy）和修复（Repair）**算子，根据算子的历史表现自适应地选择算子组合。

算法

ALNS：

1. 生成初始解 s₀
2. 初始化所有算子的权重 wᵢ = 1
3. 重复：
   a. 根据权重选择破坏算子 d 和修复算子 r
      （轮盘赌选择）
   b. 应用 d 破坏当前解（移除部分元素）
   c. 应用 r 修复解（重新插入元素）
   d. 接受准则（如模拟退火接受准则）
   e. 更新算子权重：
      - 如果找到新最优 → wᵢ += σ₁（大奖）
      - 如果改进当前解 → wᵢ += σ₂（中奖）
      - 如果接受但未改进 → wᵢ += σ₃（小奖）
4. 直到满足终止条件

常用算子

破坏算子：

• Random Removal：随机移除 $q$ 个客户
• Worst Removal：移除成本增加最多的客户
• Shaw Removal：移除相似的客户（位置/时间相近）
• Route Removal：移除整条路线

修复算子：

• Greedy Insertion：将客户插入成本增加最少的位置
• Regret Insertion：选择后悔值最大的客户（不插最佳位置的损失）
• Random Insertion：随机插入

8.11 超启发式框架

自适应算子选择

除了ALNS的权重自适应，还有：

• **多臂赌博机（MAB）**模型：将算子选择建模为赌博机问题
• 强化学习驱动：用RL学习算子选择策略

变邻域搜索（VNS）

VNS：
1. 定义邻域结构 N₁, N₂, ..., Nₖ
2. 初始解 s
3. 重复：
   k = 1
   while k ≤ K:
     a. 扰动：在 Nₖ(s) 中随机选 s'
     b. 局部搜索：对 s' 做局部搜索得到 s''
     c. 如果 s'' 优于 s → s = s'', k = 1
     d. 否则 → k++

启发式方法选择指南

问题特征	推荐方法
小规模精确解	MIP求解器
VRP类路径优化	ALNS
连续函数优化	PSO / DE
排列组合问题	GA + 专用编码
需要快速得到解	贪心 + 局部搜索
调度问题	禁忌搜索 / GA
多峰函数	模拟退火 / DE

---

第九篇：随机优化与鲁棒优化

9.1 不确定性优化概述

现实中的优化问题通常面临不确定性：

• 未来需求不确定
• 供应能力波动
• 价格变化
• 设备故障
• 天气影响

不确定性建模方法

方法	不确定性描述	决策风格
随机规划	概率分布已知	最小化期望成本/风险
鲁棒优化	不确定集内任意取值	最坏情况下的最优
分布鲁棒	属于某个分布族	最坏分布下的最优
模糊优化	模糊隶属度	模糊决策

9.2 随机规划

两阶段随机规划

第一阶段（Here-and-Now）：在不确定参数揭示之前做出决策
第二阶段（Recourse）：不确定参数揭示后，根据实际情况做补救决策

$$\begin{array}{rl}\min & c^{T}x+{\mathbb{E}}_{\xi }[Q(x,\xi )]\\ \text{s.t.}& Ax\le b\\ & x\ge 0\end{array}$$

其中 Q ( x , ξ ) = min ⁡ { q T y ∣ W y ≥ h ( ξ ) − T ( ξ ) x ,    y ≥ 0 } Q(x, \xi) = \min\{q^T y \mid Wy \geq h(\xi) - T(\xi)x, \; y \geq 0\} Q(x,ξ)=min{qTy∣Wy≥h(ξ)−T(ξ)x,y≥0} 是第二阶段问题。

情景树方法

用有限个**情景（Scenarios）**近似不确定参数的分布：

$${\mathbb{E}}_{\xi }[Q(x,\xi )]\approx \sum _{s=1}^S{p}_{s}Q(x,{\xi }^s)$$

多阶段随机规划

阶段0     阶段1     阶段2     阶段3
  │         │         │         │
  x₀ ────── x₁ ────── x₂ ────── x₃
  │    ξ₁   │    ξ₂   │    ξ₃   │
  └─ 不确定 ─┘─ 不确定 ─┘─ 不确定 ─┘

9.3 鲁棒优化

核心思想

不假设不确定参数的概率分布，而是定义一个不确定集 $\mathcal{U}$，优化最坏情况下的表现：

$$\underset{x}{\min }\underset{\xi \in \mathcal{U}}{\max }f(x,\xi )$$

不确定集类型

类型	形式	特点
盒式集	$	\xi_i - \bar{\xi}_i
椭球集	$|\xi -\bar{\xi }{|}_2\le \Omega$	考虑相关性
多面体集	$|\xi -\bar{\xi }{|}_1\le \Gamma$	可调保守度

Bertsimas-Sim鲁棒优化

引入预算参数 $\Gamma$ 控制不确定性的保守程度：

$$\max \left\{\sum _i{a}_{ij}{\xi }_i\mid \sum _i\frac{|{\xi }_i-{\bar{a}}_{ij}|}{{\hat{a}}_{ij}}\le \Gamma ,|{\xi }_i-{\bar{a}}_{ij}|\le {\hat{a}}_{ij}\right\}$$

• $\Gamma =0$：确定性优化（不考虑不确定性）
• $\Gamma =n$：完全鲁棒（考虑所有不确定参数同时取最坏值）
• 中间值：平衡鲁棒性和最优性

9.4 分布鲁棒优化

分布鲁棒优化（DRO）假设不确定参数的概率分布属于某个模糊集（Ambiguity Set），优化最坏分布下的期望：

$$\underset{x}{\min }\underset{P\in \mathcal{P}}{\max }{\mathbb{E}}_P[f(x,\xi )]$$

模糊集的构建方法：

• 矩模糊集：已知均值和协方差
• Wasserstein球：基于Wasserstein距离
• KL散度球：基于KL散度

---

第十篇：多目标优化

10.1 多目标优化基础

现实决策中经常需要同时优化多个可能冲突的目标：

$$\begin{array}{rl}\min & F(x)=(f_1(x),f_2(x),...,f_k(x))\\ \text{s.t.}& x\in \mathcal{F}\end{array}$$

例子：

• 投资组合：收益 vs 风险
• 产品设计：性能 vs 成本
• 物流：时间 vs 成本 vs 碳排放

10.2 Pareto最优

Pareto支配

解 $x$ 支配 解 $y$（记为 $x\prec y$），当且仅当：
$$\forall i:f_i(x)\le f_i(y)\text{且}\exists j:f_j(x)<f_j(y)$$

即 $x$ 在所有目标上不差于 $y$，且至少在一个目标上严格更优。

Pareto最优

解 $x^{\ast }$ 是Pareto最优，当且仅当不存在任何解支配 $x^{\ast }$。

Pareto前沿

所有Pareto最优解对应的目标向量构成Pareto前沿（Pareto Front）：

f₂ ↑
   │
   │  ●         ← Pareto前沿
   │   ●●
   │     ●●●
   │        ●●●●
   │            ●●●●●
   │ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ → f₁
   │ (最小化f₁和f₂)

10.3 求解方法

加权和法

$$\min \sum _{i=1}^k{w}_i{f}_i(x),w_i>0,\sum w_i=1$$

优点：简单。缺点：无法找到非凸Pareto前沿上的解。

ε-约束法

选择一个目标最小化，其他目标转化为约束：

$$\min f_1(x)\text{s.t.}{f}_i(x)\le {ϵ}_i,i=2,...,k$$

通过改变 $ϵ$ 值得到Pareto前沿上的不同解。

目标规划法

设定每个目标的期望值 $T_i$，最小化与目标的偏差：

$$\min \sum _i(d_i^{+}+d_i^{-})\text{s.t.}{f}_i(x)+d_i^{-}-d_i^{+}=T_i$$

Tchebycheff标量化

min ⁡ max ⁡ i { w i ∣ f i ( x ) − z i ∗ ∣ } \min \max_i \{w_i |f_i(x) - z_i^*|\} minimax​{wi​∣fi​(x)−zi∗​∣}

其中 $z_i^{\ast }$ 是第 $i$ 个目标的最优值。可以找到非凸前沿上的解。

10.4 NSGA-II算法

**NSGA-II（Non-dominated Sorting Genetic Algorithm II）**是最流行的多目标进化算法。

核心机制

1. 非支配排序：将种群按支配关系分层（Front 1, Front 2, …）
2. 拥挤距离：同一前沿内，密度大的区域个体被淘汰，保持多样性
3. 精英保留：父代和子代合并后选择

算法流程

NSGA-II：

1. 随机初始化种群 P₀（大小N）
2. 对每一代 t：
   a. 通过选择、交叉、变异生成子代 Q_t（大小N）
   b. 合并 R_t = P_t ∪ Q_t（大小2N）
   c. 对 R_t 进行非支配排序：F₁, F₂, ...
   d. 从 F₁ 开始逐层加入新种群：
      - 如果 Σ|Fᵢ| ≤ N → 全部加入
      - 如果加入 Fₖ 后超过 N → 按拥挤距离排序取前N个
   e. P_{t+1} = 选出的N个个体

Python实现

import numpy as np

def nsga2(problem, n_pop=100, n_gen=200, crossover_rate=0.9, 
          mutation_rate=0.1):
    """
    problem: 需提供 create(), evaluate(ind) -> [f1, f2, ...],
             crossover(p1, p2), mutate(ind)
    """
    # 初始化
    pop = [problem.create() for _ in range(n_pop)]
    objectives = [problem.evaluate(ind) for ind in pop]
    
    for gen in range(n_gen):
        # 生成子代
        offspring = []
        for _ in range(n_pop // 2):
            p1, p2 = tournament_select(pop, objectives, 2)
            if np.random.random() < crossover_rate:
                c1, c2 = problem.crossover(p1, p2)
            else:
                c1, c2 = p1[:], p2[:]
            if np.random.random() < mutation_rate:
                c1 = problem.mutate(c1)
            if np.random.random() < mutation_rate:
                c2 = problem.mutate(c2)
            offspring.extend([c1, c2])
        
        # 合并父代和子代
        combined = pop + offspring
        combined_obj = [problem.evaluate(ind) for ind in combined]
        
        # 非支配排序
        fronts = fast_non_dominated_sort(combined_obj)
        
        # 选择下一代
        new_pop = []
        new_obj = []
        for front in fronts:
            if len(new_pop) + len(front) <= n_pop:
                for idx in front:
                    new_pop.append(combined[idx])
                    new_obj.append(combined_obj[idx])
            else:
                # 按拥挤距离排序
                remaining = n_pop - len(new_pop)
                cd = crowding_distance(combined_obj, front)
                sorted_front = sorted(front, key=lambda i: cd[i], reverse=True)
                for idx in sorted_front[:remaining]:
                    new_pop.append(combined[idx])
                    new_obj.append(combined_obj[idx])
                break
        
        pop = new_pop
        objectives = new_obj
    
    return pop, objectives


def fast_non_dominated_sort(objectives):
    """快速非支配排序"""
    n = len(objectives)
    domination_count = [0] * n  # 被支配次数
    dominated_set = [[] for _ in range(n)]
    fronts = [[]]
    
    for i in range(n):
        for j in range(i+1, n):
            if dominates(objectives[i], objectives[j]):
                dominated_set[i].append(j)
                domination_count[j] += 1
            elif dominates(objectives[j], objectives[i]):
                dominated_set[j].append(i)
                domination_count[i] += 1
        
        if domination_count[i] == 0:
            fronts[0].append(i)
    
    current_front = 0
    while fronts[current_front]:
        next_front = []
        for i in fronts[current_front]:
            for j in dominated_set[i]:
                domination_count[j] -= 1
                if domination_count[j] == 0:
                    next_front.append(j)
        current_front += 1
        fronts.append(next_front)
    
    return fronts[:-1]


def dominates(obj1, obj2):
    """判断obj1是否支配obj2"""
    better_in_any = False
    for a, b in zip(obj1, obj2):
        if a > b:  # 假设最小化
            return False
        if a < b:
            better_in_any = True
    return better_in_any


def crowding_distance(objectives, front):
    """计算拥挤距离"""
    n = len(front)
    distances = {i: 0 for i in front}
    
    for m in range(len(objectives[0])):
        sorted_indices = sorted(front, key=lambda i: objectives[i][m])
        distances[sorted_indices[0]] = float('inf')
        distances[sorted_indices[-1]] = float('inf')
        
        obj_range = (objectives[sorted_indices[-1]][m] - 
                     objectives[sorted_indices[0]][m])
        if obj_range == 0:
            continue
        
        for i in range(1, n-1):
            distances[sorted_indices[i]] += (
                (objectives[sorted_indices[i+1]][m] - 
                 objectives[sorted_indices[i-1]][m]) / obj_range
            )
    
    return distances

---

第十一篇：约束规划（CP）

11.1 约束规划概述

约束规划（Constraint Programming, CP） 是一种基于约束满足的建模和求解范式，特别擅长处理逻辑约束、全局约束和调度问题。

与数学规划的区别

维度	数学规划（MIP）	约束规划（CP）
建模语言	数学公式（线性/整数）	约束+变量域
求解方法	LP松弛+分支定界	约束传播+回溯搜索
擅长领域	资源分配、网络流	调度、排班、配置
约束表达	需线性化	直接表达复杂约束
全局约束	不直接支持	原生支持
最优性证明	强（LP松弛提供界）	弱（通常需要与B&B结合）

11.2 约束传播与回溯搜索

约束传播

核心思想：通过约束推理，逐步缩小每个变量的取值域。

弧相容（Arc Consistency）：对于约束 $C(X_i,X_j)$，对 $X_i$ 域中的每个值，$X_j$ 域中至少存在一个值满足约束，反之亦然。

初始：X₁ ∈ {1,2,3,4,5}, X₂ ∈ {1,2,3,4,5}
约束：X₁ < X₂

传播后：
X₁ ∈ {1,2,3,4}  （X₁不能为5，因为没有X₂>5）
X₂ ∈ {2,3,4,5}  （X₂不能为1，因为没有X₁<1）

回溯搜索（Backtracking）

深度优先搜索 + 约束传播：

1. 选择一个未赋值变量 Xᵢ
2. 从域中选一个值 v（使用值选择启发式）
3. 赋值 Xᵢ = v
4. 约束传播：缩小其他变量的域
5. 如果某个变量域为空 → 回溯
6. 递归搜索下一个变量
7. 如果所有变量都赋值 → 找到解

关键启发式

变量选择：

• MRV（Minimum Remaining Values）：选域最小的变量（最少选择，最快失败）
• Degree启发式：选约束最多的变量

值选择：

• LCV（Least Constraining Value）：选对其他变量域影响最小的值

11.3 全局约束

全局约束是CP的强大之处，能对一组变量的整体关系进行高效推理。

全局约束	含义	传播算法
AllDifferent	所有变量取不同值	匹配算法
Global Cardinality (GCC)	每个值出现指定次数	网络流
Cumulative	资源累积使用量约束	时间表传播
Circuit	变量构成一条回路	基于匹配
NoOverlap / Disjunctive	任务不重叠	时间推理
Regular	序列满足自动机	基于自动机传播
Element	$X[Y]=Z$（下标访问）	域传播

11.4 CP与MIP的对比

何时使用CP

• 调度问题（尤其是带复杂约束的）
• 排班问题
• 配置问题
• 约束逻辑编程
• 需要快速找到可行解

何时使用MIP

• 资源分配问题
• 需要最优性证明
• 有明确的线性结构
• 需要灵敏度分析

CP求解器

• Google OR-Tools CP-SAT：最先进的开源CP求解器
• IBM CP Optimizer：商业CP求解器
• MiniZinc：CP建模语言

---

第十二篇：求解器与工具链

12.1 商业求解器

求解器	开发商	擅长	特点
Gurobi	Gurobi Optimization	LP/MIP/QP/QCQP	性能最强之一，API友好
CPLEX	IBM	LP/MIP/QP/CP	行业标杆，功能全面
MOSEK	MOSEK ApS	锥优化/SDP	凸优化专精
FICO Xpress	FICO	LP/MIP	在金融领域应用广
CP Optimizer	IBM	CP/调度	调度问题特别强

12.2 开源求解器

求解器	类型	特点
HiGHS	LP/MIP	性能接近商业，目前最活跃的开源LP/MIP求解器
CBC	MIP	COIN-OR项目，成熟的开源MIP求解器
GLPK	LP/MIP	老牌开源求解器
CLP	LP	COIN-OR的LP求解器
SCIP	MIP/CP	学术界广泛使用，支持约束整数规划
OR-Tools	混合	Google开发，支持CP-SAT、LP、路由
IPOPT	NLP	内点法NLP求解器
BONMIN	MINLP	混合整数非线性规划
Couenne	MINLP	凸MINLP
ECOS	SOCP	嵌入式锥优化
SLSQP	NLP	SciPy集成
OSQP	QP	一阶方法，适合嵌入式
GEKKO	MINLP	Python友好，支持DAE

12.3 建模语言与接口

Python建模库

库	特点	适用场景
PuLP	简单易用	LP/MIP入门
CVXPY	凸优化DSL	凸优化、锥优化
Pyomo	功能全面	LP/MIP/NLP/MINLP
OR-Tools	Google出品	CP、路由、LP
Gurobipy	Gurobi官方API	高性能MIP
Python-MIP	简洁的MIP接口	MIP
DOcplex	IBM CPLEX API	LP/MIP/CP
GEKKO	在线优化	NLP/MINLP/DAE
Picat	逻辑+约束	CP/SAT

专用建模语言

语言	特点
AMPL	代数建模语言，最经典
GAMS	经济学/工程优化常用
MiniZinc	约束优化建模语言
OPL	IBM的建模语言
ZIMPL	开源，配合SCIP使用
JuMP	Julia语言的优化建模

12.4 求解器选型指南

你遇到的是什么类型的问题？
│
├── 线性规划（LP）
│   └── 推荐：HiGHS（开源）/ Gurobi（商业）
│
├── 混合整数规划（MIP）
│   ├── 小规模（<10000整数变量）
│   │   └── PuLP + CBC / HiGHS
│   ├── 中等规模
│   │   └── Gurobi / CPLEX
│   └── 大规模
│       └── Gurobi（需要调参和分解技巧）
│
├── 凸优化（QP/SOCP/SDP）
│   └── CVXPY + MOSEK / ECOS / SCS
│
├── 非线性规划（NLP）
│   ├── 光滑问题 → IPOPT / SciPy
│   └── 非光滑问题 → 次梯度 / 近端方法
│
├── 混合整数非线性（MINLP）
│   └── BONMIN / Couenne / BARON
│
├── 约束规划/调度
│   └── OR-Tools CP-SAT / IBM CP Optimizer
│
└── 组合优化/TSP/VRP
    └── ALNS / 元启发式 / OR-Tools 路由

实际建议

1. 从建模开始：先用PuLP/CVXPY快速建模，用默认求解器验证
2. 逐步升级：如果求解慢，换更好的求解器
3. 调参：预处理、启发式强度、割平面策略等
4. 分解：如果仍然太大，考虑Benders分解、列生成等
5. 启发式：对于NP-hard问题，好的启发式可能比精确算法更实用

---

第十三篇：实战案例

13.1 生产调度问题

问题描述

工厂有 $n$ 个工件需要在 $m$ 台机器上加工，每个工件的每道工序有确定的加工时间和机器要求，目标是最小化最大完工时间（Makespan）。

数学建模

$$\begin{array}{rl}\min & C_{\max }\\ \text{s.t.}& C_{\max }\ge C_{ij},\forall i,j\\ & C_{ij}\ge s_{ij}+p_{ij},\forall i,j\text{（完成时间）}\\ & s_{ij}\ge C_{i,j-1},\forall i,j>1\text{（工序顺序）}\\ & \text{机器不冲突约束（大M法或位置变量）}\end{array}$$

Python实现（OR-Tools）

from ortools.sat.python import cp_model

def job_shop_scheduling(jobs_data):
    """
    jobs_data: [[(machine, duration), ...], ...]
    jobs_data[i][j] = (machine, duration) 表示工件i的第j道工序
    """
    model = cp_model.CpModel()
    
    num_jobs = len(jobs_data)
    all_machines = set()
    for job in jobs_data:
        for machine, _ in job:
            all_machines.add(machine)
    num_machines = len(all_machines)
    
    horizon = sum(dur for job in jobs_data for _, dur in job)
    
    # 创建变量
    all_tasks = {}
    for i, job in enumerate(jobs_data):
        for j, (machine, duration) in enumerate(job):
            suffix = f'_{i}_{j}'
            start = model.NewIntVar(0, horizon, f'start{suffix}')
            end = model.NewIntVar(0, horizon, f'end{suffix}')
            interval = model.NewIntervalVar(start, duration, end, f'interval{suffix}')
            all_tasks[i, j] = {
                'start': start, 'end': end, 
                'interval': interval, 'machine': machine
            }
    
    # 工序顺序约束
    for i, job in enumerate(jobs_data):
        for j in range(len(job) - 1):
            model.Add(all_tasks[i, j+1]['start'] >= all_tasks[i, j]['end'])
    
    # 机器不冲突约束
    for machine in all_machines:
        intervals = []
        for i, job in enumerate(jobs_data):
            for j, (m, _) in enumerate(job):
                if m == machine:
                    intervals.append(all_tasks[i, j]['interval'])
        model.AddNoOverlap(intervals)
    
    # 最小化Makespan
    makespan = model.NewIntVar(0, horizon, 'makespan')
    model.AddMaxEquality(makespan, [all_tasks[i, len(job)-1]['end'] 
                                     for i, job in enumerate(jobs_data)])
    model.Minimize(makespan)
    
    # 求解
    solver = cp_model.CpSolver()
    solver.parameters.max_time_in_seconds = 60
    status = solver.Solve(model)
    
    if status in (cp_model.OPTIMAL, cp_model.FEASIBLE):
        print(f"Makespan: {solver.Value(makespan)}")
        for i, job in enumerate(jobs_data):
            for j in range(len(job)):
                task = all_tasks[i, j]
                print(f"  工件{i} 工序{j}: 开始={solver.Value(task['start'])}, "
                      f"结束={solver.Value(task['end'])}, "
                      f"机器={task['machine']}")

# 示例：3个工件，3台机器
jobs = [
    [(0, 3), (1, 2), (2, 2)],  # 工件0：机器0加工3→机器1加工2→机器2加工2
    [(0, 2), (2, 1), (1, 4)],  # 工件1
    [(1, 4), (2, 3)],          # 工件2
]
job_shop_scheduling(jobs)

13.2 物流配送优化

车辆路径问题（VRP）实现

from ortools.constraint_solver import routing_enums_pb2, pywrapcp

def vrp_solver(distance_matrix, num_vehicles, depot):
    """求解VRP问题"""
    manager = pywrapcp.RoutingIndexManager(
        len(distance_matrix), num_vehicles, depot)
    routing = pywrapcp.RoutingModel(manager)
    
    # 距离回调
    def distance_callback(from_index, to_index):
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return distance_matrix[from_node][to_node]
    
    transit_callback_index = routing.RegisterTransitCallback(distance_callback)
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
    
    # 搜索参数
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
    search_parameters.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
    search_parameters.time_limit.seconds = 30
    
    # 求解
    solution = routing.SolveWithParameters(search_parameters)
    
    if solution:
        total_distance = 0
        for vehicle_id in range(num_vehicles):
            index = routing.Start(vehicle_id)
            route = []
            route_distance = 0
            while not routing.IsEnd(index):
                route.append(manager.IndexToNode(index))
                previous_index = index
                index = solution.Value(routing.NextVar(index))
                route_distance += routing.GetArcCostForVehicle(
                    previous_index, index, vehicle_id)
            route.append(manager.IndexToNode(index))
            print(f"车辆{vehicle_id}: {' → '.join(map(str, route))}, "
                  f"距离={route_distance}")
            total_distance += route_distance
        print(f"总距离: {total_distance}")

# 示例
dist = [
    [0, 10, 15, 20, 25, 30],
    [10, 0, 35, 25, 20, 15],
    [15, 35, 0, 30, 20, 25],
    [20, 25, 30, 0, 15, 10],
    [25, 20, 20, 15, 0, 12],
    [30, 15, 25, 10, 12, 0],
]
vrp_solver(dist, num_vehicles=2, depot=0)

13.3 投资组合优化

Markowitz均值-方差模型

import cvxpy as cp
import numpy as np

def portfolio_optimization(expected_returns, cov_matrix, 
                           min_return=None, risk_free=0.02):
    """
    expected_returns: 预期收益率向量
    cov_matrix: 协方差矩阵
    """
    n = len(expected_returns)
    w = cp.Variable(n)
    
    # 组合收益和风险
    portfolio_return = expected_returns @ w
    portfolio_risk = cp.quad_form(w, cov_matrix)
    
    # 约束
    constraints = [
        cp.sum(w) == 1,    # 权重和为1
        w >= 0,            # 不做空
        w <= 0.3,          # 单只股票不超过30%
    ]
    
    if min_return is not None:
        constraints.append(portfolio_return >= min_return)
    
    # 最小化风险
    prob = cp.Problem(cp.Minimize(portfolio_risk), constraints)
    prob.solve(solver=cp.SCS)
    
    return w.value, prob.value

# 有效前沿计算
def efficient_frontier(expected_returns, cov_matrix, n_points=50):
    """计算有效前沿"""
    min_ret = expected_returns.min()
    max_ret = expected_returns.max()
    target_returns = np.linspace(min_ret, max_ret, n_points)
    
    results = []
    for target in target_returns:
        try:
            w, risk = portfolio_optimization(expected_returns, cov_matrix, 
                                             min_return=target)
            results.append({
                'return': target,
                'risk': np.sqrt(risk),
                'weights': w
            })
        except:
            pass
    
    return results

13.4 选址问题

设施选址（Facility Location）

from pulp import *

def facility_location(fixed_costs, transport_costs, demands, capacities):
    """
    fixed_costs: [f_j] 各设施开设成本
    transport_costs: [[c_ij]] 运输成本矩阵
    demands: [d_i] 各客户需求
    capacities: [K_j] 各设施容量
    """
    n_facilities = len(fixed_costs)
    n_customers = len(demands)
    
    prob = LpProblem("Facility_Location", LpMinimize)
    
    # 决策变量
    open_fac = [LpVariable(f"open_{j}", cat='Binary') 
                for j in range(n_facilities)]
    assign = [[LpVariable(f"assign_{i}_{j}", lowBound=0, upBound=1) 
               for j in range(n_facilities)] for i in range(n_customers)]
    
    # 目标函数
    total_cost = (lpSum(fixed_costs[j] * open_fac[j] 
                        for j in range(n_facilities)) +
                  lpSum(demands[i] * transport_costs[i][j] * assign[i][j]
                        for i in range(n_customers) 
                        for j in range(n_facilities)))
    prob += total_cost
    
    # 每个客户的需求被满足
    for i in range(n_customers):
        prob += lpSum(assign[i][j] for j in range(n_facilities)) == 1
    
    # 只有开放的设施才能服务客户
    for i in range(n_customers):
        for j in range(n_facilities):
            prob += assign[i][j] <= open_fac[j]
    
    # 容量约束
    for j in range(n_facilities):
        prob += (lpSum(demands[i] * assign[i][j] 
                       for i in range(n_customers)) 
                 <= capacities[j] * open_fac[j])
    
    prob.solve(PULP_CBC_CMD(msg=0))
    
    print(f"总成本: {value(prob.objective)}")
    for j in range(n_facilities):
        if value(open_fac[j]) > 0.5:
            served = [i for i in range(n_customers) 
                      if value(assign[i][j]) > 0.5]
            print(f"  设施{j} 开放，服务客户: {served}")

13.5 排班问题

护士排班

from ortools.sat.python import cp_model

def nurse_scheduling(num_nurses, num_days, shifts_per_day, 
                     min_nurses_per_shift, max_consecutive_days):
    """
    求解护士排班问题
    """
    model = cp_model.CpModel()
    
    # 变量：shift[n][d][s] = 1 表示护士n在第d天的第s个班次工作
    shift = {}
    for n in range(num_nurses):
        for d in range(num_days):
            for s in range(shifts_per_day):
                shift[n, d, s] = model.NewBoolVar(f'shift_{n}_{d}_{s}')
    
    # 每个班次需要最少人数
    for d in range(num_days):
        for s in range(shifts_per_day):
            model.Add(sum(shift[n, d, s] for n in range(num_nurses)) 
                      >= min_nurses_per_shift[s])
    
    # 每个护士每天最多一个班次
    for n in range(num_nurses):
        for d in range(num_days):
            model.Add(sum(shift[n, d, s] for s in range(shifts_per_day)) <= 1)
    
    # 连续工作天数限制
    for n in range(num_nurses):
        for d in range(num_days - max_consecutive_days):
            model.Add(
                sum(shift[n, d+k, s] 
                    for k in range(max_consecutive_days + 1)
                    for s in range(shifts_per_day)) <= max_consecutive_days)
    
    # 公平性：尽量均匀分配
    # 每个护士的总班次数
    total_shifts = []
    for n in range(num_nurses):
        total = model.NewIntVar(0, num_days * shifts_per_day, f'total_{n}')
        model.Add(total == sum(shift[n, d, s] 
                               for d in range(num_days) 
                               for s in range(shifts_per_day)))
        total_shifts.append(total)
    
    # 最小化最大与最小班次数之差
    max_shifts = model.NewIntVar(0, num_days * shifts_per_day, 'max_shifts')
    min_shifts = model.NewIntVar(0, num_days * shifts_per_day, 'min_shifts')
    model.AddMaxEquality(max_shifts, total_shifts)
    model.AddMinEquality(min_shifts, total_shifts)
    model.Minimize(max_shifts - min_shifts)
    
    # 求解
    solver = cp_model.CpSolver()
    solver.Solve(model)
    
    # 输出排班表
    days = ['周一', '周二', '周三', '周四', '周五', '周六', '周日']
    shift_names = ['早班', '中班', '晚班']
    
    for n in range(num_nurses):
        schedule = []
        for d in range(num_days):
            for s in range(shifts_per_day):
                if solver.Value(shift[n, d, s]):
                    schedule.append(f"{days[d]}{shift_names[s]}")
        print(f"护士{n}: {', '.join(schedule) if schedule else '休息'}")

nurse_scheduling(num_nurses=8, num_days=7, shifts_per_day=3,
                 min_nurses_per_shift=[2, 2, 1], max_consecutive_days=5)

13.6 供应链网络设计

多级供应链优化

from pulp import *

def supply_chain_design():
    """多级供应链网络设计"""
    
    # 数据
    suppliers = ['S1', 'S2']           # 供应商
    factories = ['F1', 'F2', 'F3']    # 工厂候选
    warehouses = ['W1', 'W2']          # 仓库
    customers = ['C1', 'C2', 'C3', 'C4']  # 客户
    
    # 开设成本
    factory_open_cost = {'F1': 100, 'F2': 120, 'F3': 90}
    warehouse_open_cost = {'W1': 50, 'W2': 60}
    
    # 供应能力
    supply_capacity = {'S1': 200, 'S2': 150}
    
    # 工厂产能
    factory_capacity = {'F1': 180, 'F2': 200, 'F3': 160}
    
    # 仓库容量
    warehouse_capacity = {'W1': 150, 'W2': 170}
    
    # 客户需求
    demand = {'C1': 40, 'C2': 50, 'C3': 60, 'C4': 45}
    
    prob = LpProblem("Supply_Chain", LpMinimize)
    
    # 决策变量
    open_f = {f: LpVariable(f"open_{f}", cat='Binary') for f in factories}
    open_w = {w: LpVariable(f"open_{w}", cat='Binary') for w in warehouses}
    
    # 运输量
    x_sf = {(s,f): LpVariable(f"x_{s}_{f}", lowBound=0) 
            for s in suppliers for f in factories}
    x_fw = {(f,w): LpVariable(f"x_{f}_{w}", lowBound=0) 
            for f in factories for w in warehouses}
    x_wc = {(w,c): LpVariable(f"x_{w}_{c}", lowBound=0) 
            for w in warehouses for c in customers}
    
    # 简化的运输成本
    c_sf = {(s,f): 5 for s in suppliers for f in factories}
    c_fw = {(f,w): 3 for f in factories for w in warehouses}
    c_wc = {(w,c): 4 for w in warehouses for c in customers}
    
    # 目标：最小化总成本
    prob += (lpSum(factory_open_cost[f] * open_f[f] for f in factories) +
             lpSum(warehouse_open_cost[w] * open_w[w] for w in warehouses) +
             lpSum(c_sf[s,f] * x_sf[s,f] for s in suppliers for f in factories) +
             lpSum(c_fw[f,w] * x_fw[f,w] for f in factories for w in warehouses) +
             lpSum(c_wc[w,c] * x_wc[w,c] for w in warehouses for c in customers))
    
    # 供应能力约束
    for s in suppliers:
        prob += lpSum(x_sf[s,f] for f in factories) <= supply_capacity[s]
    
    # 流量平衡（工厂）
    for f in factories:
        prob += (lpSum(x_sf[s,f] for s in suppliers) == 
                 lpSum(x_fw[f,w] for w in warehouses))
        prob += (lpSum(x_fw[f,w] for w in warehouses) 
                 <= factory_capacity[f] * open_f[f])
    
    # 流量平衡（仓库）
    for w in warehouses:
        prob += (lpSum(x_fw[f,w] for f in factories) == 
                 lpSum(x_wc[w,c] for c in customers))
        prob += (lpSum(x_wc[w,c] for c in customers) 
                 <= warehouse_capacity[w] * open_w[w])
    
    # 客户需求满足
    for c in customers:
        prob += lpSum(x_wc[w,c] for w in warehouses) == demand[c]
    
    prob.solve(PULP_CBC_CMD(msg=0))
    
    print(f"总成本: {value(prob.objective)}")
    print("\n开放的工厂:", [f for f in factories if value(open_f[f]) > 0.5])
    print("开放的仓库:", [w for w in warehouses if value(open_w[w]) > 0.5])

supply_chain_design()

---

第十四篇：前沿方向与学习路线

14.1 机器学习+优化

学习增强分支定界

• 分支变量选择：用GNN预测哪个变量分支效果最好
• 割平面选择：学习选择最有效的割
• 节点选择：预测哪个子树最有希望

神经网络求解组合优化

• Pointer Network：用注意力机制直接输出TSP路径
• Transformer模型：求解VRP、调度等
• 图神经网络（GNN）：利用问题的图结构

优化辅助机器学习

• 超参数优化：贝叶斯优化
• 神经架构搜索（NAS）：搜索最优网络结构
• 模型压缩：通过优化剪枝和量化

可微优化

将优化问题嵌入神经网络中，通过隐式微分端到端训练：

$$\underset{\theta }{\min }\sum _i\mathcal{L}(y_i,{\text{OPT}}_{\theta }(x_i))$$

14.2 在线优化与实时决策

• 在线凸优化：对抗性环境下的在线决策
• Bandit算法：探索-利用权衡
• 模型预测控制（MPC）：实时滚动优化

14.3 量子优化

• 量子退火：D-Wave量子计算机
• QAOA：量子近似优化算法
• 变分量子本征求解器（VQE）
• 目前仍在探索阶段，实际优势尚不明确

14.4 学习路线推荐

阶段一：入门（1-2个月）

目标：理解优化基本概念，能用工具求解简单问题

学习内容：
├── 线性规划基础（标准形式、图解法）
├── 整数规划基础（0-1变量建模）
├── Python工具入门（PuLP / SciPy linprog）
└── 简单案例（背包、选址、运输问题）

推荐资源：
├── 教材：《运筹学》（清华版，第1-4章）
├── 在线课程：Coursera "Discrete Optimization"（墨尔本大学）
└── 练习：PuLP官方教程

阶段二：进阶（2-3个月）

目标：掌握主流算法原理，能建模和求解中等规模问题

学习内容：
├── 单纯形法原理
├── 分支定界法
├── 动态规划
├── 图论与网络优化
├── 非线性规划基础
├── MIP建模技巧（大M法、SOS等）
└── 求解器使用（Gurobi/CPLEX，OR-Tools）

推荐资源：
├── 教材：《Introduction to Linear Optimization》(Bertsimas & Tsitsiklis)
├── 教材：《Combinatorial Optimization》(Papadimitriou & Steiglitz)
├── Gurobi官方文档和示例
└── OR-Tools官方教程

阶段三：高级（3-6个月）

目标：深入算法原理，能处理大规模/复杂结构问题

学习内容：
├── 对偶理论深入
├── Benders分解、列生成
├── 拉格朗日松弛
├── 元启发式算法（GA, SA, ALNS等）
├── 约束规划
├── 随机优化/鲁棒优化
├── 多目标优化
└── 实际项目实践

推荐资源：
├── 教材：《Integer Programming》(Wolsey)
├── 教材：《Convex Optimization》(Boyd & Vandenberghe)
├── ALNS论文和代码
└── 学术会议：INFORMS, CPAIOR

阶段四：精通（持续学习）

目标：前沿研究，优化+AI结合

学习内容：
├── 学习增强的优化算法
├── 大规模分布式优化
├── 实时/在线优化
├── 领域特定建模技巧
└── 阅读最新论文，参与开源社区

14.5 推荐资源

经典教材

教材	内容侧重	适合人群
《运筹学》清华大学	全面入门	零基础
Introduction to Linear Optimization (Bertsimas)	LP理论与算法	入门-进阶
Integer Programming (Wolsey)	整数规划理论	进阶
Convex Optimization (Boyd)	凸优化理论	进阶-高级
Combinatorial Optimization (Papadimitriou)	组合优化	进阶
Metaheuristics (Gendreau & Potvin)	元启发式	进阶
Introduction to Stochastic Programming (Birge & Louveaux)	随机规划	高级
Robust Optimization (Ben-Tal et al.)	鲁棒优化	高级

在线课程

课程	平台	特点
Discrete Optimization	Coursera（墨尔本大学）	最好的优化入门课，有编程作业
Linear Programming	edX（MIT）	MIT开放课程
Optimization Methods	Coursera（斯坦福）	凸优化入门
Gurobi Webinars	Gurobi官网	实用建模技巧
COIN-OR教程	COIN-OR	开源优化工具

开源项目与社区

资源	说明
OR-Tools	Google优化工具包，文档丰富
COIN-OR	开源运筹优化软件集合
NEOS Server	在线优化求解平台
OR-Library	经典优化测试问题集合
awesome-decision-optimization	GitHub上优化资源汇总

学术期刊与会议

名称	说明
Operations Research	顶刊
Mathematical Programming	理论顶刊
INFORMS Journal on Computing	计算方向
European Journal of Operational Research	应用广泛
INFORMS年会	最大的运筹学会议
CPAIOR	约束规划+整数规划
ISMP	数学规划国际会议

Python生态速查

# ============ 快速参考 ============

# 1. 线性规划 (LP)
from scipy.optimize import linprog
from pulp import LpProblem, LpVariable, lpSum

# 2. 凸优化
import cvxpy as cp

# 3. 非线性规划
from scipy.optimize import minimize

# 4. 混合整数规划 (MIP)
from pulp import *  # 简单MIP
import gurobipy as gp  # Gurobi（商业，最快）
from ortools.sat.python import cp_model  # CP-SAT（开源，强）

# 5. 约束规划/调度
from ortools.sat.python import cp_model
# 或用 MiniZinc

# 6. 路径优化 (TSP/VRP)
from ortools.constraint_solver import routing_enums_pb2, pywrapcp

# 7. 多目标优化
from pymoo.algorithms.moo.nsga2 import NSGA2
import pymoo

# 8. 元启发式
import pygmo  # 多种算法
# 或自行实现 SA/GA/PSO/ACO

# 9. 图算法
import networkx as nx  # 最短路径、最小生成树、网络流等

# 10. 高级建模
import pyomo.environ as pyo  # 通用代数建模语言

---

附录A：常见优化问题速查表

问题	类型	复杂度	推荐方法
线性规划	LP	P	单纯形/内点法
背包问题	IP	NP-hard	DP/MIP
指派问题	IP	P	匈牙利算法
最短路径	网络	P	Dijkstra/Bellman-Ford
最大流	网络	P	Dinic/Edmonds-Karp
最小费用流	网络	P	SSP算法
最小生成树	网络	P	Prim/Kruskal
TSP	组合	NP-hard	B&C/ALNS/LKH
VRP	组合	NP-hard	ALNS/OR-Tools
调度(JSSP)	组合	NP-hard	CP/B&B
设施选址	IP	NP-hard	B&C/Benders
集合覆盖	IP	NP-hard	B&C/LP松弛
投资组合	QP/QCQP	P(凸)	CVXPY
短期调度	MIP	NP-hard	B&C/CP
网络设计	MIP	NP-hard	Benders分解

附录B：术语对照表

英文	中文	说明
Optimization	优化
Objective Function	目标函数
Decision Variable	决策变量
Constraint	约束
Feasible Region	可行域
Optimal Solution	最优解
Local Optimum	局部最优
Global Optimum	全局最优
Linear Programming (LP)	线性规划
Integer Programming (IP)	整数规划
Mixed-Integer Programming (MIP)	混合整数规划
Nonlinear Programming (NLP)	非线性规划
Dynamic Programming (DP)	动态规划
Simplex Method	单纯形法
Branch and Bound (B&B)	分支定界法
Cutting Plane	割平面法
Branch and Cut	分支切割法
Lagrangian Relaxation	拉格朗日松弛
Benders Decomposition	Benders分解
Duality	对偶理论
KKT Conditions	KKT条件
Convex Optimization	凸优化
Heuristic	启发式
Metaheuristic	元启发式
Simulated Annealing (SA)	模拟退火
Genetic Algorithm (GA)	遗传算法
Particle Swarm Optimization (PSO)	粒子群优化
Ant Colony Optimization (ACO)	蚁群算法
Tabu Search (TS)	禁忌搜索
Differential Evolution (DE)	差分进化
ALNS	自适应大邻域搜索
Constraint Programming (CP)	约束规划
Robust Optimization	鲁棒优化
Stochastic Programming	随机规划
Pareto Optimality	Pareto最优
Multi-objective Optimization	多目标优化
Makespan	最大完工时间	调度问题常用指标
Shadow Price	影子价格	对偶变量的经济学解释
Presolve	预处理
Warm Start	热启动

---

附录C：算法复杂度速查

算法	时间复杂度	适用问题
单纯形法	指数最坏，实际多项式	LP
内点法	$O(n^{3.5}\log (1/ϵ))$	LP/QP/凸优化
分支定界	指数	IP/MIP
Dijkstra	$O((n+m)\log n)$	单源最短路径（非负权）
Bellman-Ford	$O(nm)$	单源最短路径（含负权）
Floyd-Warshall	$O(n^3)$	全源最短路径
匈牙利算法	$O(n^3)$	指派问题/二分匹配
最大流(Dinic)	$O(n^{2}m)$	网络流
动态规划	看具体问题	多种问题
遗传算法	$O(GNP)$	$G$=代数，$N$=种群，$P$=评估代价
模拟退火	无理论保证	通用组合优化