组合优化深度指南

想象一下,你是一家快递公司的调度员。每天早上睁开眼,就要面对一个让人头秃的问题:怎么安排路线才能让所有包裹在最短时间内送完?怎么分配车辆才能省油钱?怎么安排快递员的工作时间才能让所有人都满意?这些问题听起来像是拍脑袋就能解决的,但实际操作起来,你会发现它们都有一个共同的特点——选项多得吓人,而且每个选择都影响着其他选择。

这就是组合优化要研究的问题。

组合优化是运筹学的核心领域,专门对付这种”离散”的决策问题。跟连续变量不同,组合优化处理的是那些必须”一个一个选”的问题——选或不选、去这家店还是那家店、先干这个还是先干那个。正因为选项是离散的,所以解空间往往是天文数字,暴力枚举根本不可能。

比如说旅行商问题——你要去10个城市推销,每个城市只能去一次,最后回到出发点。听起来挺简单对吧?但实际上你有9!=362880种走法。要是20个城市呢?那就是2432902008176640000种可能,宇宙诞生到现在都没这么长时间给你枚举。

所以组合优化研究的核心就是:怎么在这些海量的可能性中,快速找到最优解(或者至少是足够好的解)。


组合优化入门:离散世界的优化问题

什么叫做”组合”?

“组合”这两个字在数学里通常意味着”排列组合”——从一堆东西里挑出一些,按照某种顺序排起来。组合优化就是研究这类”挑挑选选排排序”问题的学科。

你可能觉得这不就是枚举吗?有那么复杂吗?

问题在于,组合问题的解空间增长得太快了。假设你有n个东西要安排,最简单的排列问题都有n!种可能。当n=10的时候,10!=3628800,这个数字已经需要好好想想了。当n=20的时候,20!≈2.4×10^18,这个数字比银河系的星星还多。你用最快的超级计算机,枚举完所有排列也需要几百年。

组合优化要做的,就是在不解枚举的情况下,找到最好的那一个(或者至少是一个很不错的)。

组合优化问题的特点

组合优化问题通常有三个特征:

第一个,决策是离散的。 比如选择哪条路线、分配哪个工人、接受哪个订单。变量通常只能取整数或布尔值(0或1)。这跟线性规划不一样,线性规划的变量可以是任何实数。

第二个,解空间是有限的但规模指数级。 再多的选项也是有限的,但架不住增长太快。n个城市的旅行商问题是有限的,但你永远不可能枚举完。

第三个,存在约束条件。 现实问题总是有各种限制——时间不够用、预算有限、每个工人每天只能工作8小时……这些约束把解空间切成一块一块的,我们只能在满足约束的解里面挑最优的。

常见组合优化问题一览

组合优化领域有很多经典问题,它们之所以经典,是因为:

  1. 本身有重要的实际应用价值
  2. 解决它们的技术可以推广到其他问题
  3. 它们是检验算法能力的”试金石”

背包问题:给你一堆东西,每个东西有价值有重量,你的背包只能装一定重量。怎么选才能让价值最大?这是资源分配的抽象。

指派问题:n个人n项工作,每个人做每项工作有不同的效率。怎么分配才能让总效率最高?这就是经典的0-1规划问题。

旅行商问题:刚才说过了,是路径优化的标杆。

调度问题:一堆任务要安排到时间轴上,怎么排才能让总完成时间最短?这关系到工厂排产、项目管理。

装箱问题:怎么把一堆物品装进最少数量的箱子?这是仓储物流的经典问题。

这些问题的共同特点是:看起来简单,但实际求解极其困难。它们大多属于NP难问题——目前没有已知的多项式时间算法能保证找到最优解。

但这不意味着我们束手无策。组合优化研究出了各种技术:精确算法能在某些情况下找到最优解,近似算法能在多项式时间内给出一个有保证的解,元启发式算法能在大规模问题上给出不错的解。

接下来的章节,我们就来逐一了解这些技术。


线性规划入门:从几何直观到单纯形法

线性规划是什么?

在说组合优化之前,必须先聊一下它的”老祖宗”——线性规划(Linear Programming,LP)。

线性规划是优化领域最基础的问题类型。它的特点是:目标函数是线性的,所有约束条件也是线性的,变量可以取任意实数。听起来挺简单的对吧?但这个”简单”的问题却是现代运筹学的基石。

什么叫”线性的”?简单说就是没有平方、没有log、没有sin函数,就是一堆数字乘以变量,然后加起来。比下面这个就是典型的线性规划:

最大化:3x₁ + 5x₂

满足:
2x₁ + 3x₂ ≤ 12  (资源A的限制)
x₁ + 2x₂ ≤ 8    (资源B的限制)
x₁, x₂ ≥ 0      (不能生产负数的产品)

这里x₁和x₂可以理解成两种产品的生产数量,约束条件表示不能超过资源限制。目标是想办法让利润最大化。

可行域的几何图像

如果把x₁和x₂画成坐标系,横轴是x₁,纵轴是x₂,那么约束条件会划出一块区域。这块区域就是可行域——所有满足约束的解都在这个区域内。

线性规划的可行域有个很好的性质:它一定是凸多边形(更准确地说,是凸多面体)。凸的意思是,随便在区域里挑两个点,连一条线,线上的所有点都还在区域内。这个性质太重要了,它保证了一个关键结论:

线性规划的最优解必定出现在可行域的”角落”上。

这个”角落”在数学上叫做极点顶点。你可以想象一个凸多边形,最优点一定在某个顶点上,不可能在多边形内部。

单纯形法:沿着边界爬山

1947年,George Dantzig提出了单纯形法,这是求解线性规划的革命性算法。

单纯形法的思路很直观:

  1. 先找到可行域的一个顶点
  2. 看看从当前顶点出发,沿着边走到相邻顶点会不会让目标函数值更好
  3. 如果会,就走过去;否则,当前顶点就是最优解

就像在一个凸多边形里,你要找到让某个方向函数值最大的点。从一个角落开始,看看相邻的角落哪个更好,有就跳槽,没有就待着——因为再往里走只会更差。

def simplex(A, b, c):
    # 单纯形法的核心循环
    while True:
        # 计算检验数(reduced costs)
        # 检验数 < 0 说明从这个方向走能找到更好的解
        for each non-basic variable x_j:
            if reduced_cost[j] < 0:
                entering = j
                break
        
        if all reduced_cost >= 0:
            return current_solution  # 没有更好的方向了,这就是最优解
        
        # 找离基变量(最小比例测试)
        for each basic variable x_i:
            if A[:, entering] > 0:
                ratio = b[i] / A[i, entering]
                if ratio < min_ratio:
                    leaving = i
        
        # 交换进基离基变量
        swap(entering, leaving)

关于单纯形法,有几个有意思的事实:

  • 理论上:单纯形法是指数时间的。Klee和Minty构造了一个”坏例子”,让单纯形法走遍所有顶点才找到最优。
  • 实践中:单纯形法非常快,平均是多项式时间。
  • 平滑分析:后来有人用随机扰动的平滑分析证明,单纯形法的期望时间是多项式的。

所以不要被理论吓到,单纯形法在实际中表现极好,是求解线性规划的主流方法。

对偶理论:每个问题都有一个”镜像”

对偶理论是线性规划里一个深刻的概念。每个线性规划问题都有一个”对偶问题”,原问题和对偶问题的最优目标值是相等的。

还是上面那个最大化利润的例子:

原问题:
最大化:3x₁ + 5x₂
满足:2x₁ + 3x₂ ≤ 12, x₁ + 2x₂ ≤ 8, x₁,x₂ ≥ 0

对偶问题:
最小化:12y₁ + 8y₂  
满足:2y₁ + y₂ ≥ 3, 3y₁ + 2y₂ ≥ 5, y₁,y₂ ≥ 0

这里的y₁和y₂可以理解成资源的”单价”——如果你不想生产产品,也可以把资源卖掉。约束条件说的是,每种产品的利润不能超过生产它需要的资源价值。

对偶定理告诉我们:如果原问题有最优解,对偶问题也有最优解,而且两个目标值相等。这个性质在经济学上很有意义——“资源分配能赚多少钱”等于”资源本身值多少钱”。

对偶理论还有一个实际用途:提供下界。在后面要讲的分支定界法中,对偶问题(或者LP松弛)的解可以帮助我们剪枝,加快搜索速度。


整数规划:为什么”整数”让问题变难?

整数约束——看似简单的限制

整数规划(Integer Programming,IP)跟线性规划的区别只有一个:要求变量取整数值。

整数规划:
最大化:3x₁ + 5x₂
满足:2x₁ + 3x₂ ≤ 12, x₁ + 2x₂ ≤ 8
x₁, x₂ ∈ ℤ (必须整数!)

就这么一个小改动,问题就从一个”容易”的问题变成了”困难”的问题。

你说这不是很合理吗?现实世界本来就是离散的——买三件还是四件商品,派一个人还是两个人。但数学上,这个限制造成了本质性的困难。

线性规划 vs 整数规划

线性规划的最优解可以出现在可行域的任意位置。整数规划的可行解只能取离散的点,这些离散点散布在凸多面体的顶点、子顶点、甚至凸多面体内部的某些位置。

这就带来一个问题:最优解可能不在LP松弛的顶点上。

看一个具体例子:

最大化:x₁ + x₂
满足:x₁ + 2x₂ ≤ 3.5
x₁, x₂ ≥ 0,整数

LP松弛的最优解是x₁=3.5, x₂=0,目标值3.5。但这是小数解!整数可行解只能是(0,0)、(1,0)、(2,0)、(3,0)、(0,1)、(1,1)……其中最优的整数解是(3,0)或(1,1),目标值都是3。

注意,LP松弛的最优值(3.5)是整数规划最优值(3)的上界。这个gap虽然不大,但有时gap可以非常大。

NP难——问题的固有困难

整数规划是NP难的。这意味着:

  1. 除非P=NP,否则不存在多项式时间的精确算法
  2. 很多重要问题(旅行商、调度、指派)都可以规约到整数规划
  3. 实际问题通常规模一大就很难求最优解

但”NP难”不是说完全没有办法。我们有各种技术来求解整数规划:

  • 精确算法:分支定界、分支切割、割平面法,能在某些情况下找到最优解
  • 近似算法:有理论保证,但解不一定是最优
  • 启发式算法:没有理论保证,但通常能找到不错的解
  • 元启发式:遗传算法、模拟退火等,在大规模问题上表现不错

下面我们重点讲精确求解技术。


分支定界法:系统地搜索整数解空间

核心思想

分支定界(Branch and Bound)是求解整数规划的经典框架。它的思想可以用”分而治之”来概括:

  1. 分支(Branch):把原问题分成若干个子问题,每个子问题约束更多一点
  2. 定界(Bound):计算每个子问题能得到的最好结果(上限或下限)
  3. 剪枝(Prune):如果某个子问题不可能比已知的最优解更好,就扔掉它

这就好比在一棵巨大的决策树里搜索最优叶节点。树的根是原问题,每个分支代表对一个变量的约束(比如x₁≤2或x₁≥3),叶子代表一个具体的整数解。分支定界的任务是在不遍历所有叶子的情况下,找到最好的那个。

LP松弛:重要的下界来源

分支定界之所以有效,关键在于LP松弛

把整数约束”松弛”掉,变量可以从整数变成任意实数。LP松弛的可行域比原问题(整数约束的)大,所以LP最优值是整数规划最优值的下界(对于最小化问题)。

原问题:min c^Tx, Ax ≤ b, x ∈ ℤⁿ
LP松弛:min c^Tx, Ax ≤ b, x ≥ 0

为什么下界有用?因为如果某个子问题的LP松弛值已经比当前最好的整数解还要差,那这个子问题就可以直接剪掉——再怎么分支,它的整数解也不可能是最优的。

分支定界算法框架

def branch_and_bound(ilp):
    # 初始化全局最优解
    best_value = float('inf')
    best_solution = None
    
    # 待处理的问题队列
    node_queue = [ilp]
    
    while node_queue:
        # 取一个子问题来处理
        current_problem = node_queue.pop()
        
        # 求LP松弛
        lp_solution = solve_simplex(current_problem)
        
        # 定界:如果比当前最优还差,剪掉
        if lp_solution.value >= best_value:
            continue
        
        # 检查是否是整数解
        if is_integer(lp_solution):
            # 找到了更好的整数解!
            best_value = lp_solution.value
            best_solution = lp_solution
            continue
        
        # 分支:选一个非整数变量
        fractional_var = select_fractional_var(lp_solution)
        
        # 创建两个子问题
        left = current_problem.add_constraint(
            fractional_var <= floor(lp_solution[fractional_var])
        )
        right = current_problem.add_constraint(
            fractional_var >= ceil(lp_solution[fractional_var])
        )
        
        node_queue.push(left)
        node_queue.push(right)
    
    return best_solution

加速技巧

分支定界的效率取决于你怎么”分”和怎么”界”。好的策略能让搜索树小很多:

节点选择:先处理哪个子问题?

  • 深度优先:先往一边走到底,快速找到整数解,更新上界
  • 最佳优先:总是处理下界最好的节点,更接近最优但可能搜索深度很大
  • 混合策略:结合两者,开始用深度优先找到可行解,之后用最佳优先改进

变量选择:在哪个变量上分支?

  • 最分数变量:选小数部分最接近0.5的变量——这种变量”最不确定”
  • 最大影响变量:选目标函数系数绝对值最大的变量——对解影响最大
  • 强分支:计算几个候选变量的分支效果,选择能最大降低下界的

割平面:在每个节点添加割平面收紧LP松弛,让下界更紧。分支定界加上割平面就是分支切割算法,是现代求解整数规划的主流方法。


割平面法:添加有效不等式缩小搜索空间

割平面的直觉

割平面法(Cutting Plane Method)的思想很直观:LP松弛的可行域太大了,我们需要”切掉”一些不包含整数解的区域,让松弛域收紧。

想象你要在一个凸多边形里找一个整数点。LP最优解落在某个位置,但这个小数点不是整数点。割平面就是在LP最优解附近画一条”切割线”,把包含小数点的区域切掉,但保证不切掉任何整数可行点。

这样的切割多了,最优LP解就会慢慢”挤向”整数点,最终落在整数上。

Gomory切割

Ralph Gomory在1958年提出了第一个系统性的割平面方法。思路很巧妙:

从LP最优解表中取一行(对应一个基本变量)。这一行代表一个等式约束:

x_B + Σ a_ij * x_j = b_i

把每个系数拆成整数部分和小数部分。整数部分移到左边,小数部分移到右边:

x_B + Σ ⌊a_ij⌋ * x_j - ⌊b_i⌋ = {b_i} - Σ {a_ij} * x_j

左边是整数(因为所有项都是整数),右边有一个正的常数{b_i}(0到1之间),减去一堆非负项。所以右边必须是正的。

这就给出了一个Gomory切割

Σ {a_ij} * x_j ≥ {b_i}

这条不等式把包含当前LP最优解的区域切掉了——因为在那个点上,左边的值为0,而{b_i}是正数,不等式不成立。

割平面的收敛

纯割平面法的收敛性不太好。理论上可能需要指数多条割平面才能收敛到整数解。

实践中,割平面通常跟分支定界结合使用:

  1. 求LP松弛
  2. 如果解是整数,结束
  3. 生成割平面,加到问题里
  4. 回到步骤1

这就是分支切割框架(Branch and Cut),是现代MILP求解器的核心技术。CPLEX、Gurobi等商业求解器都实现了大量高效的割平面生成策略。


匈牙利算法:二分图最大匹配

问题场景

指派问题在现实中非常常见:

  • n个工人,n项工作,每个工人做每项工作有不同的效率/成本
  • n辆车,n个送货点,每辆车去每个点的成本不同
  • 运动员和比赛项目,每个运动员参加不同项目的能力不同

数学上,这就是二分图完美匹配问题:左边n个顶点代表工人,右边n个顶点代表工作,我们需要找一个完美匹配,让总权重最小(或最大)。

直接枚举有n!种匹配方式,当n大起来根本无法枚举。

匈牙利算法的巧妙之处

1955年,Harold Kuhn提出了匈牙利算法。算法核心只有三步,但每一步都有深刻的数学原理支撑:

  1. 行归约:每行减去该行的最小值
  2. 列归约:每列减去该列的最小值
  3. 覆盖所有零:用最少的线覆盖矩阵中的所有零元素

如果能覆盖的线数等于n,就找到了最优解;如果不够,就调整矩阵继续迭代。

def hungarian(cost_matrix):
    n = len(cost_matrix)
    
    # Step 1: 行归约
    for i in range(n):
        min_val = min(cost_matrix[i])
        cost_matrix[i] = [x - min_val for x in cost_matrix[i]]
    
    # Step 2: 列归约
    for j in range(n):
        min_val = min(cost_matrix[i][j] for i in range(n))
        for i in range(n):
            cost_matrix[i][j] -= min_val
    
    # Step 3-4: 迭代覆盖和调整
    while True:
        row_cover, col_cover = find_min_vertex_cover(cost_matrix)
        
        if len(row_cover) + len(col_cover) == n:
            # 找到最优匹配
            return extract_matching(cost_matrix)
        
        # 找到未覆盖元素中的最小值
        min_uncovered = find_min_uncovered(cost_matrix)
        
        # 调整矩阵
        for i in range(n):
            for j in range(n):
                if i in row_cover and j not in col_cover:
                    cost_matrix[i][j] -= min_uncovered
                elif i not in row_cover and j in col_cover:
                    cost_matrix[i][j] += min_uncovered

为什么有效?

匈牙利算法的正确性依赖于一个关键定理:如果覆盖零元素的最少直线数等于n,那么这些零元素对应一个最优匹配。

这背后的直觉是:行归约和列归约减少了矩阵的总成本,同时不改变任何完美匹配的相对成本差。当覆盖线数达到n时,说明我们可以通过只选择零元素(成本已经归零的边)来构造一个完美匹配,这个匹配的总成本等于归约减少的总量,因此一定是最优的。

时间复杂度是O(n³),这对于二分图匹配来说是非常高效的。


最大流最小割:网络流问题的经典定理

网络流问题

假设你有一个网络,从源点s发送物品到汇点t,每条边有容量限制,问最多能发送多少物品?

这就是最大流问题(Maximum Flow Problem)。形式化一下:

  • 给定一个有向图G=(V,E)
  • 每条边(e)有容量c(e) ≥ 0
  • 源点s,汇点t
  • 求从s到t的最大可行流量

可行流满足:

  1. 容量约束:每条边的流量 ≤ 容量
  2. 流量守恒:除了s和t,每个点的流入量 = 流出量

Ford-Fulkerson方法

求最大流的基本方法是Ford-Fulkerson算法:

  1. 找到一条从s到t的增广路径(还有剩余容量的路径)
  2. 沿这条路推送尽可能多的流量
  3. 重复直到没有增广路径

增广路径的寻找可以用BFS/DFS,这就是Edmonds-Karp算法,时间复杂度O(VE²)。

更好的实现是Dinic算法,利用分层图的概念,时间复杂度O(V²E)。

最大流最小割定理

这里有一个非常漂亮的定理:

最大流的流量 = 最小割的容量

是把图的顶点分成两部分S和T,满足s∈S,t∈T。割的容量是所有从S到T的边的容量之和。

这个定理的直观解释:任何流都必须穿过任何割,就像水管网络里的水流必须经过某个截面。最大流受到最小”瓶颈”的限制,而那个瓶颈就是最小割。

最大流最小割定理有很多应用:图像分割(Graph Cut)、网络可靠性、匹配问题的转化……


旅行商问题TSP:组合优化中的皇冠问题

TSP是什么?

旅行商问题(Traveling Salesman Problem)是组合优化领域最著名的问题之一:

给定n个城市和两两之间的距离,求访问所有城市恰好一次并返回起点的最短路径。

这个问题在实际中有大量应用:物流配送路线规划、电路板钻孔顺序优化、DNA测序……

TSP的复杂度

TSP是NP难的。即使是度量TSP(距离满足三角不等式),也是APX难的——不太可能有常数近似比的多项式时间算法。

但这不意味着我们放弃了。对于小规模TSP(几百个城市),可以用精确算法求解;对于大规模TSP,启发式算法能找到接近最优的解。

精确求解:分支定界

TSP的分支定界需要一些特殊的加速技术:

  • 下界计算:可以用1-树松弛(相当于把问题变成度量TSP的最小生成树),或者LKH算法中的启发式下界
  • 分支策略:固定或不固定某条边
  • 割平面:subtour elimination约束、DFJ formulation

对于几千个城市的TSP,现代求解器已经能在可接受时间内找到最优解。

近似算法:Christofides算法

对于度量TSP,有一个漂亮的多项式时间近似算法——Christofides算法,近似比是1.5:

  1. 求最小生成树
  2. 找到度为奇数的顶点
  3. 求这些顶点的最小完美匹配
  4. 合并得到欧拉回路
  5. 走捷径得到哈密顿回路

这是TSP近似算法的重大突破,至今仍是最好的已知近似比(除非P=NP)。

启发式求解

对于实际问题中的大规模TSP,启发式算法更实用:

  • 贪心最近邻:从某个城市开始,每次访问最近的未访问城市
  • 2-opt/3-opt:通过边交换改进解
  • Lin-Kernighan-Helsgaun(LKH):迭代地用k-opt改进,是最强的TSP启发式之一
  • 遗传算法、蚁群算法:群体智能方法

调度问题:加工顺序的优化

调度的现实需求

调度问题在制造业、项目管理、计算机系统中有广泛应用。核心是:在时间和资源约束下,安排作业的执行顺序以优化某个目标。

典型的调度问题参数:

  • 作业:需要完成的 tasks
  • 机器:执行作业的资源
  • 处理时间:每个作业在每台机器上需要的时间
  • 截止时间:作业必须完成的最晚时间
  • 释放时间:作业何时可以开始处理

调度问题的分类

根据机器环境和约束,调度问题有很多分类:

分类标准类型
机器环境单机、并行机、作业车间、开放车间
作业特征单工序、多工序;释放时间、截止时间
优化目标最小化最大完工时间、最小化总延迟、加权完工时间和
目标符号正则(只依赖完工时间)、非正则

P||C_max:并行机调度

最经典的调度问题之一:n个作业,m台相同的平行机器,每个作业处理时间p_j,目标是最小化最大完工时间(makespan)C_max。

下界容易计算:

LB = max{ max_j p_j, (Σ p_j) / m }
  • 第一个下界是最长作业时间,没有机器能在更短时间内完成它
  • 第二个下界是总工作量均分给所有机器的最快可能

LPT(Longest Processing Time first)算法

  1. 按处理时间降序排列作业
  2. 依次分配到当前负载最小的机器

LPT的近似比是4/3 - 1/(3m),对于大 m 接近 1.333。这意味着LPT解最多比最优解差33%。

单机调度:EDD规则

单机调度有个经典规则——最早截止时间优先(Earliest Due Date, EDD)

按作业的截止时间从早到晚排序。

EDD能最小化最大延迟(maximum lateness),这个结论很容易用邻位交换证明:

假设有两个相邻作业A、B,截止时间d_A ≤ d_B,但B在A前面。如果交换它们的顺序,C_A表示交换前A的完工时间,那么交换后B的完工时间会变小,A的完工时间会变大。由于d_A ≤ d_B,B的延迟一定减小,A的延迟可能增大但不会超过原来的B。整体上最大延迟不会变差。


组合优化在ML中的应用:结构化预测与特征选择

结构化预测

机器学习中的很多问题需要输出结构化对象:序列、树、图。这些输出的空间是组合的,无法用简单的回归或分类解决。

条件随机场(CRF) 把结构化预测建模为在图结构上的推断:

P(y|x) = (1/Z(x)) * exp(Σ_k θ_k * f_k(x, y))

推断等价于在标记图上求解最大后验(MAP)推断——本质上是一个组合优化问题。

序列标注是最常见的结构化预测任务:词性标注、命名实体识别。解码通常用Viterbi算法(动态规划),时间复杂度O(T * S²),T是序列长度,S是状态数。

特征选择

特征选择是从大量候选特征中挑出有用的子集,是降维的重要手段。

包装法(Wrapper)把特征选择当成子集搜索问题:

  • 前向搜索:每次加一个最好的特征
  • 后向搜索:每次删一个最差的特征
  • 递归特征消除:迭代地删除最不重要的特征

这本质上是组合优化问题,可以用分支定界加速搜索。

嵌入式方法(Lasso、树模型的特征重要性)在模型训练过程中隐式地进行特征选择,是更高效的做法。

图神经网络中的组合优化

图神经网络(GNN)的很多任务涉及组合结构:社区检测、子图匹配、图划分。

组合优化 + 深度学习是当前的研究热点:

  • Pointer Network:用注意力机制学习组合问题的求解
  • Graph Neural Network + 强化学习:学习更好的搜索策略
  • 神经组合优化:端到端学习问题求解器

整数规划实战:用PuLP/OR-Tools求解实际问题

PuLP:Python中的建模工具

PuLP是Python中轻量级的优化建模库,可以跟CBC、CPLEX、Gurobi等求解器配合使用。

安装

pip install pulp

基本用法

from pulp import *
 
# 创建问题
prob = LpProblem("My_Problem", LpMinimize)
 
# 决策变量
x = LpVariable.dicts("x", range(n), cat='Binary')
y = LpVariable("y", lowBound=0)
 
# 目标函数
prob += lpSum(c[i] * x[i] for i in range(n))
 
# 约束
prob += lpSum(a[i] * x[i] for i in range(n)) <= b
prob += y >= 0
 
# 求解
status = prob.solve(PULP_CBC_CMD(msg=0))
 
# 输出
print(f"Status: {LpStatus[status]}")
print(f"Optimal value: {value(prob.objective)}")

员工排班问题实战

排班问题在现实中很常见:医院、工厂、客服中心都需要合理安排员工工作时间。

问题描述

  • 每天需要一定数量的员工
  • 每个员工有工作时间限制和偏好
  • 目标是满足需求的前提下最小化人力成本
from pulp import *
 
# 参数
days = range(7)  # 一周7天
workers = range(10)  # 10个员工
need = [3, 2, 4, 3, 5, 5, 2]  # 每天需要的员工数
cost = {i: 100 for i in workers}  # 每个人工作一天的成本
max_days = {i: 5 for i in workers}  # 每人最多工作天数
 
prob = LpProblem("Staff_Scheduling", LpMinimize)
 
# 决策变量
x = LpVariable.dicts("work", 
                     [(i, d) for i in workers for d in days], 
                     cat='Binary')
 
# 目标函数:最小化总成本
prob += lpSum(cost[i] * x[(i, d)] 
              for i in workers for d in days)
 
# 约束1:满足每天需求
for d in days:
    prob += lpSum(x[(i, d)] for i in workers) >= need[d]
 
# 约束2:每个员工最多工作max_days天
for i in workers:
    prob += lpSum(x[(i, d)] for d in days) <= max_days[i]
 
# 求解
prob.solve(PULP_CBC_CMD(msg=1))
 
# 输出结果
print(f"状态: {LpStatus[prob.status]}")
print(f"最优成本: {value(prob.objective)}")
 
# 打印排班表
for d in days:
    assigned = [i for i in workers if value(x[(i, d)]) > 0.5]
    print(f"第{d+1}天: 员工 {assigned}")

车辆路径规划问题(VRP)

VRP是TSP的扩展:多辆车从仓库出发,服务一组客户点,求最小成本的路线。

# 简化的CVRP(Capacitated VRP)
from pulp import *
 
# 参数
n_customers = 20
n_vehicles = 4
vehicle_capacity = 50
customer_demand = {i: random.randint(5, 15) for i in range(1, n_customers+1)}
distance = {(i, j): random.randint(1, 100) 
            for i in range(n_customers+1) 
            for j in range(n_customers+1) if i != j}
 
prob = LpProblem("CVRP", LpMinimize)
 
# 决策变量
x = LpVariable.dicts("route", 
                     [(i, j, k) for i in range(n_customers+1) 
                      for j in range(n_customers+1) 
                      for k in range(n_vehicles)
                      if i != j],
                     cat='Binary')
 
u = LpVariable.dicts("load", 
                     [(i, k) for i in range(n_customers+1) 
                      for k in range(n_vehicles)],
                     lowBound=0)
 
# 目标函数
prob += lpSum(distance[(i, j)] * x[(i, j, k)] 
              for i in range(n_customers+1) 
              for j in range(n_customers+1) 
              for k in range(n_vehicles) 
              if i != j)
 
# MTZ formulation约束(篇幅所限略)
# ...
 
prob.solve()

OR-Tools:Google的优化工具箱

OR-Tools是Google开源的优化求解器,特别擅长约束规划和调度问题。

from ortools.sat.python import cp_model
 
model = cp_model.CpModel()
 
# 决策变量
x = [model.NewIntVar(0, 10, f'x{i}') for i in range(3)]
 
# 约束
model.Add(x[0] + x[1] + x[2] <= 15)
model.Add(x[0] * 2 + x[1] == 10)
 
# 目标
model.Minimize(x[0] + x[1] + x[2])
 
solver = cp_model.CpSolver()
status = solver.Solve(model)
 
if status == cp_model.OPTIMAL:
    print([solver.Value(x[i]) for i in range(3)])

OR-Tools的优势在于约束编程(CP-SAT),对复杂的组合约束支持很好,调度、排班、资源分配问题用它很方便。


启发式方法:贪心、局部搜索、模拟退火

启发式为什么必要?

精确算法(分支定界、割平面)的问题是:对于大规模问题,即使有很强的剪枝,搜索树还是可能太大。

启发式方法放弃最优性保证,转而追求”快速找到一个足够好的解”。这不是偷懒,而是现实主义——在工业界,一个能在1分钟内给出80分解的算法,往往比一个跑10小时给出100分解的算法更有用。

贪心算法

贪心是最简单的启发式:每步选当前最优,不考虑全局。

背包问题的贪心

def greedy_knapsack(items, capacity):
    # 按价值重量比排序
    items.sort(key=lambda x: x.value / x.weight, reverse=True)
    
    total_value = 0
    total_weight = 0
    selected = []
    
    for item in items:
        if total_weight + item.weight <= capacity:
            selected.append(item)
            total_weight += item.weight
            total_value += item.value
    
    return selected, total_value

贪心对某些问题效果不错(如最小生成树),但对背包问题,贪心解可能远非最优(特别是物品价值差异大时)。

局部搜索

局部搜索从一个可行解出发,不断做”小改进”:

def local_search(current_solution):
    while True:
        neighbor = get_best_neighbor(current_solution)
        if neighbor.value <= current_solution.value:
            break  # 没有改进,停止
        current_solution = neighbor
    return current_solution
 
def get_best_neighbor(solution):
    # 尝试所有"邻居"解,返回最好的
    best = solution
    for neighbor in generate_neighbors(solution):
        if neighbor.value > best.value:
            best = neighbor
    return best

邻居的定义很重要

  • TSP:交换两条边(2-opt)、移动一个城市的位置( relocate)
  • 调度:交换两个作业的顺序
  • 背包:加入或移除一个物品

局部搜索的问题是容易陷入局部最优——可能有很多更好的解,但局部搜索只看到邻居,无法跳出去。

模拟退火

模拟退火(Simulated Annealing)通过”概率接受差解”来跳出局部最优:

import random
import math
 
def simulated_annealing(initial_solution, initial_temp=1000, cooling_rate=0.995):
    current = initial_solution
    best = current
    temperature = initial_temp
    
    while temperature > 1:
        neighbor = get_random_neighbor(current)
        
        delta = neighbor.value - current.value
        
        # 如果是更好的解,或者按概率接受差解
        if delta > 0 or random.random() < math.exp(delta / temperature):
            current = neighbor
            
            # 更新最优
            if current.value > best.value:
                best = current
        
        temperature *= cooling_rate
    
    return best

物理类比:金属冷却时,原子可能跳到能量更高的位置(差解),随着温度降低,跳的概率越来越小,最终稳定在能量最低的状态(最优解)。

关键参数

  • 初始温度:太高会让算法几乎随机搜索,太低会让算法陷入局部最优
  • 冷却率:控制温度下降速度,常用0.95-0.999
  • 终止温度:一般设为1或0.1

模拟退火的关键洞察是:短期看是退步的移动,长期看可能导向更好的解


元启发式算法:遗传算法、粒子群在组合优化中的应用

遗传算法:模拟自然选择

遗传算法(Genetic Algorithm)受达尔文进化论启发:种群中适应性强的个体更可能存活并产生后代,后代继承父母的优点,可能产生更好的个体。

基本框架

def genetic_algorithm(population_size=100, generations=1000):
    # 初始化种群
    population = [create_random_solution() for _ in range(population_size)]
    
    for gen in range(generations):
        # 评估适应度
        fitness = [evaluate(ind) for ind in population]
        
        # 选择父母
        parents = selection(population, fitness, k=2)
        
        # 交叉产生后代
        children = []
        for _ in range(population_size // 2):
            child1, child2 = crossover(parents[0], parents[1])
            children.extend([child1, child2])
        
        # 变异
        children = [mutate(child) for child in children]
        
        # 替换种群
        population = selection(population + children, 
                               [evaluate(ind) for ind in population + children],
                               k=population_size)
    
    return best_individual(population)

TSP的遗传算法设计

  • 编码:路径表示,如[3,1,4,2]表示经过3→1→4→2
  • 交叉:顺序交叉(OX)、部分映射交叉(PMX)
  • 变异:交换两个城市、reverse一段子路径

粒子群优化

粒子群优化(Particle Swarm Optimization)模拟鸟群觅食:每只鸟(粒子)记住自己到过的最好位置和整个种群到过的最好位置,根据这两个信息调整飞行方向。

def particle_swarm_optimization(n_particles=30, n_iterations=1000):
    # 初始化
    particles = [{'position': random_position(),
                   'velocity': random_velocity(),
                   'best_position': None,
                   'best_value': -inf} 
                  for _ in range(n_particles)]
    global_best_position = None
    global_best_value = -inf
    
    for _ in range(n_iterations):
        for p in particles:
            # 评估
            value = evaluate(p['position'])
            
            # 更新个人最优
            if value > p['best_value']:
                p['best_position'] = p['position']
                p['best_value'] = value
            
            # 更新全局最优
            if value > global_best_value:
                global_best_position = p['position']
                global_best_value = value
            
            # 更新速度和位置
            r1, r2 = random(), random()
            p['velocity'] = (w * p['velocity'] 
                           + c1 * r1 * (p['best_position'] - p['position'])
                           + c2 * r2 * (global_best_position - p['position']))
            p['position'] = p['position'] + p['velocity']
    
    return global_best_position

参数说明

  • w:惯性权重,控制速度保持之前方向的程度
  • c1:认知系数,向个人最优学习的程度
  • c2:社会系数,向全局最优学习的程度

元启发式的实际应用建议

  1. 问题适配:没有”万能”的元启发式。TSP适合用GA/ACO(蚁群),连续优化用PSO/DE(差分进化)

  2. 参数调优:元启发式对参数敏感。可以用网格搜索、贝叶斯优化、或自适应参数控制

  3. 混合方法:实践中,往往把元启发式和局部搜索结合起来(memetic algorithm),效果最好

  4. 早停策略:如果解质量在多代内没有提升,可以提前停止


计算复杂度总结

理解组合优化问题的复杂度,对于选择求解方法至关重要:

问题复杂度最优算法近似比
线性规划(LP)P单纯形/内点法1
整数线性规划(ILP)NP难分支定界-
二分图匹配P匈牙利算法1
一般图匹配PEdmonds算法1
最大流PDinic/推送-重标记1
最小生成树PPrim/Kruskal1
旅行商(一般)NP难分支定界-
旅行商(度量)APX难Christofides1.5
调度 P‖C_maxAPX难LPT4/3
背包伪多项式动态规划完全多项式近似方案(FPTAS)
集合覆盖NP难贪心ln(n)
顶点覆盖NP难贪心/精确2

关键理解

  • P问题:有多项式时间精确算法
  • NP难问题:精确算法是指数时间的,只能用近似或启发式
  • APX难:不太可能有常数近似比的多项式算法
  • 近似比:近似解和最优解的比值上限,越接近1越好

关键词速览

关键词解释
线性规划目标函数和约束均为线性的优化问题,变量连续
整数线性规划变量受限为整数的LP,是NP难问题
单纯形法求解LP的多项式时间算法,实际运行很快
割平面法通过切割平面收紧松弛域,通常与其他方法结合
分支定界系统化枚举搜索树,是现代MILP求解器的基础
匈牙利算法O(n³)求解二分图最大权完美匹配
调度问题作业与机器的资源分配,有丰富的理论和应用
约束满足CSP的求解框架,包括回溯搜索和弧一致性
对偶理论LP的解的边界理论,提供下界的重要工具
近似算法有理论保证的启发式,性能有最优性证明
模拟退火受金属冷却启发的全局优化算法
遗传算法模拟自然选择,能有效搜索大解空间

相关文档