稍微憋个招，聊聊为什么不能止步于会调求解器

文章背景

自从实践了整数规划建模+求解器计算的方式之后，一时间觉得找到了一把万能钥匙，啥组合优化问题都可以尝试以这种方式来解决。

但实际上，我也了解到，很多经典组合优化问题，都有比较经典的求解算法，比如针对背包问题的动态规划算法和针对指派问题的匈牙利算法等。

由此就引发了我的一个疑惑：这些经典算法的价值是什么？是否还有必要去学习和使用？

本文后续将通过详细探讨背包问题和指派问题的求解方案和方案效果，以期回答以上的疑惑。

背包问题

背包问题可以描述为：给定$n$个重量为$w_1,w_2,\cdot \cdot \cdot ,w_n$、价值为$v_1,v_2,\cdot \cdot \cdot ,v_n$的物品，有一个最大承重为$W$的背包，求这些物品中一个价值最高的子集，并且要能够装到背包中。

数学模型

整数规划

暂且不管经典算法，直接将其建模为整数规划问题。

定义$x_i$为第$i$个物品是否被放入背包，其值为0时，表示不放入，其值为1时，表示放入。

此时，可以建立如下的整数规划模型
m a x ∑ i = 1 n v i x i s.t ∑ i = 1 n w i x i ≤ W , i = 1 , 2 , . . . , n x i ∈ { 0 , 1 } , i = 1 , 2 , . . . , n max \quad \sum_{i=1}^nv_{i}x_{i} \\ \text{s.t} \quad \sum_{i=1}^nw_ix_{i}≤W, \quad i=1,2,...,n \\ \nonumber x_{i} \in \{0,1\} ,\quad i=1,2,...,n\\ maxi=1∑n​vi​xi​s.ti=1∑n​wi​xi​≤W,i=1,2,...,nxi​∈{0,1},i=1,2,...,n

动态规划

至于求解背包问题的经典算法，刷过LeetCode的人应该都知道是动态规划算法。鉴于动态规划的算法原理并不是本文的重点，所以这里只给个链接，有兴趣的人可以自己查看。

仿真实验

以下代码是基于Python实现的求解背包问题的整数规划算法和动态规划算法。通过调整$N$的值，可以改变背包问题的规模。因此我们可以很直观对比出两种算法在不同问题规模下的方案结果，包括最优解的质量和求解的速度，从而评估算法的能力。

from ortools.linear_solver import pywraplp
import numpy as np
import time


def calc_by_ortools(N, w, v, W):
    # 声明ortools求解器，使用SCIP算法
    solver = pywraplp.Solver.CreateSolver('SCIP')

    # 优化变量，0-1变量
    x = {}
    for j in range(N):
        x[j] = solver.IntVar(0, 1, 'x[%i]' % j)

    # 目标函数
    obj_expr = [v[j][0] * x[j] for j in range(N)]
    solver.Maximize(solver.Sum(obj_expr))

    # 约束条件
    cons_expr = [w[j][0] * x[j] for j in range(N)]
    solver.Add(solver.Sum(cons_expr) <= W)

    # 模型求解
    status = solver.Solve()

    # 打印模型结果
    if status == pywraplp.Solver.OPTIMAL:
        # 求解成功，打印最优目标函数值
        print('ortools, best_f =', solver.Objective().Value())

    else:
        # 求解不成功，提示未收敛
        print('not converge.')


def calc_by_dp(weight, value, bag_weight):
    # 初始化: 全为0
    dp = [0] * (bag_weight + 1)

    # 先遍历物品, 再遍历背包容量
    for i in range(len(weight)):
        for j in range(bag_weight, weight[i][0] - 1, -1):
            # 递归公式
            dp[j] = max(dp[j], dp[j - weight[i][0]] + value[i][0])
    print('dp, best_f =', dp[-1])


if __name__ == '__main__':
    # 设置随机种子，确保每次运行生成的随机数相同
    np.random.seed(0)

    # 设定物品数量N，重量w，价值v，背包可承重W
    N = 1000
    w = np.random.randint(1, 10, (N, 1))
    v = np.random.randint(1, 100, (N, 1))
    W = int(N / 10)
    print('N = ', N)

    # 使用ortools求解，并统计计算耗时
    t0 = time.time()
    calc_by_ortools(N, w, v, W)
    print('ortools计算耗时：{}'.format(time.time() - t0))

    # 使用动态规划方法求解，并统计计算耗时
    t1 = time.time()
    calc_by_dp(w, v, W)
    print('dp计算耗时：{}'.format(time.time() - t1))

以下表格是两个算法在不同$N$上的详细表现数据，其中ortools指的是整数规划算法，dp指的是动态规划算法。

从解的质量来看，两个算法都能找到全局最优解，因此是无差别的。

但是在求解效率上，两个算法就有很大的差异了：N<1000时，ortools的计算时间大于dp，但绝对数值都很小；N=1000时，ortools和dp的计算时间相差已经较小；继续增大$N$后，ortools的计算时间便小于dp，而且dp计算时间的增长速度显然快于ortools。

N	算法	最优解	耗时,s
10	ortools	89	0.0085
	dp	89	0.0000
100	ortools	616	0.0117
	dp	616	0.0007
1000	ortools	6154	0.0424
	dp	6154	0.0629
10000	ortools	60509	0.4257
	dp	60509	7.6769
100000	ortools	617258	5.111
	dp	617258	730.8

乍一看，这么对比，好像没什么问题。但是总结的时候突然想起来，ortools是基于C++写的，而dp是用Python写的，dp会不会在程序语言上吃了亏，所以换成Java再试试（别问我为什么不用C++，问就是不会）。

以下是Java版本的算法实现，整体逻辑和Python是一致的，就不赘述了。

import java.util.Random;
import com.google.ortools.Loader;
import com.google.ortools.linearsolver.MPConstraint;
import com.google.ortools.linearsolver.MPObjective;
import com.google.ortools.linearsolver.MPSolver;
import com.google.ortools.linearsolver.MPVariable;

public class ZeroOnePack {

    // 预加载本地库
     static {
        Loader.loadNativeLibraries();
    }

    public static void DP(int W, int N, int[] weight, int[] value){
        //动态规划
        int[] dp = new int[W +1];
        for(int i=1;i<N+1;i++){
            //逆序实现
            for(int j = W; j>=weight[i-1]; j--){
                dp[j] = Math.max(dp[j-weight[i-1]]+value[i-1],dp[j]);
            }
        }

        // 打印最优解
        System.out.println("DP, best_f: " + dp[W]);

    }

    public static void orToolsMethod(int W, int N, int[] weight, int[] value){
         // 声明求解器
         MPSolver solver = MPSolver.createSolver("SCIP");
         if (solver == null) {
             System.out.println("Could not create solver SCIP");
             return;
         }

         // 优化变量
         MPVariable[] x = new MPVariable[N];
         for (int j = 0; j < N; ++j) {
             x[j] = solver.makeIntVar(0.0, 1, "");
         }

         // 目标函数
         MPObjective objective = solver.objective();
         for (int j = 0; j < N; ++j) {
             objective.setCoefficient(x[j], value[j]);
         }

         // 约束条件
         objective.setMaximization();
         MPConstraint constraint = solver.makeConstraint(0, W, "");
         for (int j = 0; j < N; ++j) {
             constraint.setCoefficient(x[j], weight[j]);
         }
         // 模型求解
         MPSolver.ResultStatus resultStatus = solver.solve();

         if (resultStatus == MPSolver.ResultStatus.OPTIMAL) {
             // 求解成功，打印最优目标函数值
             System.out.println("ortools, best_f = " + objective.value());
         } else {
             // 求解不成功，提示未收敛
             System.err.println("The problem does not have an optimal solution.");
         }
     }

    public static void main(String[] args) {
         //设置随机种子，确保每次运行生成的随机数相同
         Random rand =new Random(0);

        // 设定物品数量N，重量weight，价值value，背包可承重W
         int N = 1000000;

         int[] weight=new int[N];
         for(int i=0;i<weight.length;i++){
             weight[i]= rand.nextInt(10) + 1;
         }

         int[] value=new int[N];
         for(int i=0;i<value.length;i++){
             value[i]= rand.nextInt(100) + 1;
         }

         int W = (int) N / 10;
         System.out.println("N = " + N);

         // 使用ortools求解，并统计计算耗时
         long start = System.currentTimeMillis();
         orToolsMethod(W, N, weight, value);
         System.out.println("cost time: " + (System.currentTimeMillis() - start) + " ms");

         // 使用动态规划方法求解，并统计计算耗时
         start = System.currentTimeMillis();
         DP(W, N, weight, value);
         System.out.println("cost time: " + (System.currentTimeMillis() - start) + " ms");

     }
}

以下表格是两个算法（Java版本）在不同$N$上的详细表现数据，这里需要注意一下第四列，耗时的单位是ms，Python版本里是s。

因为使用了Java，dp的计算效率显著增长，例如$N=100000$时，java耗时2s，但是python耗时高达730s。

不过即使如此，算法对比的结论依然不变：两个算法都可以找到最优解；但是在计算效率上，ortools先落后再领先。

N	算法	最优解	耗时,ms
10	ortools	78	9
	dp	78	0
100	ortools	481	10
	dp	481	0
1000	ortools	6224	17
	dp	6224	3
10000	ortools	60442	81
	dp	60442	22
100000	ortools	603439	2405
	dp	603439	2039
200000	ortools	1207108	3128
	dp	1207108	6994
1000000	ortools	6037100	15614
	dp	6037100	135898

结果分析

可以使用动态规划算法能够得到背包问题的全局最优解，是因为背包问题满足最优化原理和无后效性原则。在求解过程中，动态规划算法的时间复杂度是$O(nW)$，但由于$W$只是一个输入数据，它可以表示成输入规模$n$的指数形式，所以是个伪多项式算法，即不是多项式算法。因此，随着$N$的增加，dp的计算时间增加的非常快。

而整数规划能够得到全局最优解，是毋庸置疑的；但由于本身也不是多项式算法，所以随着$N$的增大，其计算耗时也增加了很多，只不过从对比数据来看，相比动态规划算法，整数规划的效率还是更高一些。

指派问题

分析完了背包问题，再来研究一下指派问题。

指派问题可以描述为：$n$个人分配$n$项任务，一个人只能分配一项任务，一项任务只能分配给一个人，将一项任务分配给一个人是需要支付报酬，求如何分配任务，保证支付的报酬总数最小。

数学模型

整数规划

设定支付的报酬为矩阵$C_{n\times n}$，其中$c_{i,j}$表示第$i$个人分配第$j$项任务时需要支付的报酬。

定义$x_{i,j}$为第$i$个人是否被分配第$j$项任务，其值为0时，表示不被分配，其值为1时，表示被分配。

此时，可以建立如下的数学规划模型
m i n ∑ i = 1 n ∑ j = 1 n c i , j x i , j s.t ∑ j = 1 n x i , j = 1 , i = 1 , 2 , . . . , n ∑ i = 1 n x i , j = 1 , j = 1 , 2 , . . . , n x i , j ∈ { 0 , 1 } , i , j = 1 , 2 , . . . , n min \quad \sum_{i=1}^n\sum_{j=1}^nc_{i,j}x_{i,j} \\ \text{s.t} \quad \sum_{j=1}^nx_{i,j}=1, \quad i=1,2,...,n \\ \nonumber \sum_{i=1}^nx_{i,j}=1, \quad j=1,2,...,n \\ \nonumber x_{i,j} \in \{0,1\} ,\quad i,j=1,2,...,n\\ mini=1∑n​j=1∑n​ci,j​xi,j​s.tj=1∑n​xi,j​=1,i=1,2,...,ni=1∑n​xi,j​=1,j=1,2,...,nxi,j​∈{0,1},i,j=1,2,...,n

匈牙利算法

指派问题的经典算法是匈牙利算法。不过该算法实现起来就不像动态规划算法那么容易，所以还是找现成的工具包更省事一些，本文用的是scipy.optimize包中的linear_sum_assignment模块。该算法的原理可以参考文献：On implementing 2D rectangular assignment algorithms，据称其本质还是匈牙利算法，但由于不影响文中结论，所以笔者自身未仔细考究哈~

仿真实验

以下代码是基于Python实现的求解指派问题的整数规划算法和匈牙利算法。通过调整$N$的值，可以改变指派问题的规模。因此我们可以很容易对比出这两种算法在不同问题规模下的方案结果，包括解的质量和求解的速度，从而评估算法的能力。

from ortools.linear_solver import pywraplp
from scipy.optimize import linear_sum_assignment
import numpy as np
import time


def calc_by_ortools(C):
    # 声明ortools求解器，使用SCIP算法
    solver = pywraplp.Solver.CreateSolver('SCIP')
    m = C.shape[0]
    n = C.shape[1]

    # 优化变量，0-1变量
    x = {}
    for i in range(m):
        for j in range(n):
            x[i, j] = solver.IntVar(0, 1, 'x[%i,%i]' % (i, j))

    # 目标函数
    obj_expr = [C[i][j] * x[i, j] for i in range(m) for j in range(n)]
    solver.Minimize(solver.Sum(obj_expr))

    # 约束条件
    for i in range(m):
        cons_expr = [x[i, j] for j in range(n)]
        solver.Add(solver.Sum(cons_expr) == 1)

    for j in range(n):
        cons_expr = [x[i, j] for i in range(m)]
        solver.Add(solver.Sum(cons_expr) == 1)

    # 模型求解
    status = solver.Solve()

    # 打印模型结果
    if status == pywraplp.Solver.OPTIMAL:

        # 求解成功，打印最优目标函数值
        print('ortools, best_f =', solver.Objective().Value())

    else:
        # 求解不成功，提示未收敛
        print('not converge.')


def calc_by_scipy(C):
    # 调用工具包：linear_sum_assignment
    row_ind, col_ind = linear_sum_assignment(C)
    # 打印最优目标函数值
    print('scipy, best_f =', cost[row_ind, col_ind].sum())


if __name__ == '__main__':
    # 设置随机种子，确保每次运行生成的随机数相同
    np.random.seed(0)

    # 设定报酬矩阵的维度
    N = 1000
    # 报酬范围是10~100间的随机值
    cost = np.random.randint(10, 100, (N, N))
    print('N = ', N)

    # 使用ortools求解，并统计计算耗时
    t0 = time.time()
    calc_by_ortools(cost)
    print('ortools计算耗时：{}'.format(time.time() - t0))

    # 使用求解scipy中的 modified Jonker-Volgenant algorithm求解，并统计计算耗时
    t1 = time.time()
    calc_by_scipy(cost)
    print('scipy计算耗时：{}'.format(time.time() - t1))

以下表格是两个算法在不同$N$上的详细表现数据，其中ortools指的是整数规划算法，scipy指的是匈牙利算法。

从解的质量来看，两个算法也总是能找到全局最优解，因此是无差别的。

在求解时间上，ortools始终高于scipy；而且随着$N$的增大，scipy的求解时间增长较慢，而ortools却增长的非常快。

N	算法	最优解	耗时,s
10	ortools	222	0.0136
	scipy	222	0
50	ortools	621	0.1599
	scipy	621	0.0001
100	ortools	1087	0.9516
	scipy	1087	0.0003
300	ortools	3034	9.9593
	scipy	3034	0.0047
500	ortools	5004	24.89
	scipy	5004	0.0118
1000	ortools	10000	177.5
	scipy	10000	0.0396

结果分析

整数规划算法的结果就不分析了，和背包问题基本一致，简单看一下匈牙利算法。

从算法原理上来说，它是针对指派问题的特点，找到的一个多项式算法，所以耗时非常短。

总结

其实分析分析后，结论已经呼之欲出了：将组合优化问题建模为整数规划问题来求解，本质上使用的是一种通用方案，只是由于很多公司都致力于迭代优化求解器的效率，所以目前来看，这个通用方案的整体表现还不错；但那些针对特定问题的特定算法，可以理解为一种个性化解决方案，旨在通过利用问题自身的特征，探查更高效的解决方案。

基于这个理解，在实际问题的求解中，应该优先将问题建模为有个性化求解算法同时问题复杂度是多项式的经典问题；其次是建模为有个性化求解算法但问题复杂度不是多项式的经典问题，此时需要对比经典算法和整数规划的效率和精度；最后才是直接建模为整数规划问题。

参考文献

背包问题，动态规划Python代码：https://blog.csdn.net/m0_51370744/article/details/127120649

背包问题，动态规划java代码：https://blog.csdn.net/baidu_41602099/article/details/110383230

指派问题和匈牙利算法:https://zhuanlan.zhihu.com/p/103125599

指派问题scipy算法原理：https://sci-hub.se/10.1109/TAES.2016.140952