【0基础运筹学】【附代码】指派问题——MILP、匈牙利算法（两种形式）详解

小猪快跑爱摄影

已于 2024-09-12 11:41:03 修改

阅读量1.3k

点赞数 32

分类专栏： MILP 文章标签：算法

于 2024-09-02 16:22:01 首次发布

本文链接：https://blog.csdn.net/ymzhu385/article/details/141822685

版权

MILP 专栏收录该内容

8 篇文章 5 订阅

订阅专栏

文章目录

相关文献
参考文献
问题概述
匈牙利算法和KM算法
整数规划（MILP）

作者：小猪快跑

基础数学&计算数学，从事优化领域6年+，主要研究方向：MIP求解器、整数规划、随机规划、智能优化算法

本文将介绍指派问题常见算法——匈牙利算法、KM算法（两种形式）、MILP建模。

本文从简单的例子，用通俗的推导过程带着读者一步一步攻克看似复杂的KM算法。

如有错误，欢迎指正。如有更好的算法，也欢迎交流！！！——@小猪快跑

参考文献

M:/Work/CLASSES/COMP572_Fall2006/Notes/Hungarian/Matching.dvi (ust.hk)
[Assignment Problem and Hungarian Algorithm (topcoder.com)](https://www.topcoder.com/thrive/articles/Assignment Problem and Hungarian Algorithm)
HungarianAlgorithm.com - Solve the Assignment Problem

问题概述

指派问题是一种经典的优化问题，主要涉及将一定数量的任务分配给相同数量的员工（或资源），以最小（大）化总成本（或时间、距离等）。这种问题通常可以用一个n×n的矩阵来表示，矩阵中的每个元素 $c_{ij}$ 表示第 $i$ 个员工完成第 $j$ 个任务的成本。

例如，假设我们有 3 个员工和 3 个任务，需要为每个员工分配一个任务，使得总时间最短。

下面的矩阵显示了员工和任务的每种组合所需的时间（以分钟为单位）。工作用 J1、J2 和 J3 表示，员工用 W1、W2 和 W3 表示。

	J1	J2	J3
W1	52	10	15
W2	10	53	10
W3	25	10	32

在这里插入图片描述

匈牙利算法和KM算法

概述

匈牙利算法和KM算法实际上指的是同一种算法，它们都是用来解决指派问题的高效算法。这里的“匈牙利算法”通常是指由Harold Kuhn在1955年提出的算法版本，而“KM算法”则是指James Munkres在1957年对该算法所做的改进版本。

匈牙利算法：最初由Harold Kuhn提出，他基于匈牙利数学家Dénes Kőnig和Jenő Egerváry的工作，提出了这一算法。它提供了一种解决指派问题的有效方法，特别是当成本矩阵是正方形的时候。
KM算法（Kuhn-Munkres算法）：是匈牙利算法的一个改进版本，由James Munkres在1957年提出的。Munkres的贡献在于提供了一个更加详细的算法实现，使其更容易理解和实施。

匈牙利算法

我们先来看一个更简单的问题——二分图的最大匹配问题。

工作用 J1、J2 和 J3 表示，员工用 W1、W2 和 W3 ，如果员工和工作之间有线相连，则代表员工能胜任这份工作。

如图所示，员工 W1 能胜任 J1和 J3，员工 W2 胜任 J1、J2、J3，员工 W3 只能胜任 J3。

我们先给员工 W1 分配任务，假如是 J1

我们再给员工 W2 分配任务，假如是 J1，这时候 J1 就有两个员工做了，冲突了

一个很自然的思路就是让员工 W1 换成任务 J3

我们再给员工 W3 分配任务，假如是 J3，这时候 J3 就有两个员工做了，冲突了。

我们让员工 W1 换成任务 J1。这时候 J1 就有两个员工做了，冲突了。

我们让员工 W2 换成任务 J2。

<img src="【0基础运筹学】指派问题——MILP、匈牙利算法.assets/匈牙利算法_3.png" alt="匈牙利算法_3"

我们容易发现一个规律，在解决冲突（collision）的时候，我们都是从左边未匹配顶点开始，然后一左一右交错经过若干个匹配顶点，最终到达对面集合的一个未匹配顶点的路径（一般称之为增广路）。之后我们把绿色的线删除，就得到新的匹配。

KM算法 - 预备知识

我们尝试着先去理解几个简单的事实。

先从一个更简单的例子看起：两个员工和两个任务。如果员工 W1 今天效率高了，做任何工作都能快1分钟，我们容易发现最优匹配结果不会变化。

员工对其所有任务快慢X分钟不影响匹配结果。反之，任务对其所有员工快慢X分钟不影响匹配结果。

证明思路：如果上述问题任务可以重复做的话（每个员工仅能匹配一个任务），目标函数可以简单的看成` 左图{(W1, J1), (W1, J2)} + {(W2, J1), (W2, J2)} = {52, 10} + {10, 53} = 1 + {51, 9} + {10, 53} = 1 + 右图{(W1, J1), (W1, J2)} + {(W2, J1), (W2, J2)}`，其中{}表示只选择其中一个值。

原问题是左图的最小化时间指派问题。右图是把边权值取反，目标改为最大化时间，他的最优匹配结果和原问题一致。

KM算法形式1 - 调整顶标

【敲黑板，划重点了！！！】如果我们给每个顶点标记数字，且任意边的顶标和 >= 边权值，而我们刚好能找到一个完美匹配的任意边的顶标和 = 边权值，那么这个完美匹配就是解。

我们来具体解释一下上面这句话。比如下图中:
1 = (W1, J1) = W1 + J1
7 = (W2, J2) = W2 + J2
2 = W1 + J2 >= (W1, J2) = 2
6 = W2 + J1 >= (W2, J1) = 3
于是，(W1, J1), (W2, J2) 就是解

证明其实也很简单。
【任意边的顶标和 >= 边权值】 => 【任意匹配的顶标和 >= 边权值和】
我们的目标就是寻找【边权值和最大值】，所以如果能找到一个【匹配的顶标和 = 边权值】，那么这个匹配的边权值和一定就是最大值，就是我们要找的匹配。

于是我们现在来分析本文开始的指派问题。

我们先进行取反，问题变成MAX问题。之后通过员工对其任务增加相等的时间，得到如下等价问题。

我们最终想找到问题的解，其实就是找一个顶点标记数字，满足边的顶标和 >= 边权值，完美匹配的任意边的顶标和 = 边权值。所以问题就在于如何找到这样的顶点标记数字。

我们可以先找到边的顶标和 >= 边权值的一种顶点标记数字（这个只要数字足够大就可以），之后通过调整顶标得到完美匹配的任意边的顶标和 = 边权值

首先，我们分别对每个员工的所有边权值取最大，然后标记在顶点左边。右边的点全部标记为0。容易看出，边的顶标和 >= 边权值。

如果说每个员工分配到的工作都是自己最擅长的（边权值最大的），刚好人手一份工作，自然是最优解。

那我们不妨先试试看，尽量让每个员工都分配到自己最擅长的。

比如 W1 最擅长 J2，W2 擅长 J1，W3 擅长 J2。这时候我们发现 W1 和 W3 都擅长 J2，冲突了。

【敲黑板，划重点了！！！】

常见的方法就是让其中一个员工换一份工作。

比如说 W1 擅长 J2，其次 J3，效率只差了 5。

而 W3 擅长 J2，其次 J1，效率差了 15。

显然让 W1 换成 J3，总效率降低较小。只需 W1、W3 顶标 减少 5，J2 + 5。这样总效率降低了 5，完美匹配的任意边的顶标和 = 边权值最终我们就得到了最优解如图。

探索思路其实很简单，就是两条边的有一个公共顶点，如果公共顶点（J2）增加值，两个顶点（W1、W3）降低相同的值，边权值不变，但顶标和降低。

上述的 KM 算法是通过不断修改顶标，最终找到最优解。

KM算法形式2 - 调整边权值

接下来我们来看另一种形式的 KM 算法。有点类似对偶的感觉，通过改变边权值来寻找最优解。

跟之前一样，通过员工对其任务减少相等的时间，得到如下图二的等价问题，右图是前文描述的等价问题。

我们把边权值转换成矩阵来看：

在这里插入图片描述

我们现在用最少的（行 + 列）覆盖所有的零，找到未覆盖中最小元素 5，然后未打钩的行减去这个元素，打钩的列加上这个元素。（这是因为前文提到过：员工对其所有任务快慢X分钟不影响匹配结果。反之，任务对其所有员工快慢X分钟不影响匹配结果）

在这里插入图片描述

重复上述操作直到矩阵必须用所有行列才能覆盖所有零元素（这保证了我们一定能找到一个完美匹配）。我们把矩阵画成二分图，然后用匈牙利算法找到完美匹配。这就是最终解。

学术名词

二分图
匹配、最大匹配、完美匹配、最优匹配
最小覆盖
最大独立集
交替路、增广路

代码

我们使用 scipy 中的 linear_sum_assignment来计算，也可参考munkres · PyPI，HungarianAlg.cpp

import networkx as nx
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import linear_sum_assignment


def draw_graph(cost, ax, row_ind, col_ind):
    # 获取默认颜色列表
    colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

    u_shape = cost.shape[0]
    v_shape = cost.shape[1]

    edges = dict(zip(row_ind, col_ind))

    graph = nx.Graph()
    w_nodes = [f'W{i + 1}' for i in range(u_shape)]
    j_nodes = [f'J{i + 1}' for i in range(u_shape)]
    graph.add_nodes_from(w_nodes + j_nodes)
    for u in range(u_shape):
        for v in range(v_shape):
            graph.add_edge(f'W{u + 1}', f'J{v + 1}', weight=cost[u, v])
    pos = {}
    for u in range(u_shape):
        for v in range(v_shape):
            pos[f'W{u + 1}'] = (-1, u_shape - u)
            pos[f'J{v + 1}'] = (1, v_shape - v)

    # 绘制边
    for u in range(u_shape):
        for v in range(v_shape):
            point = (f'W{u + 1}', f'J{v + 1}')
            edge_color = colors[0] if edges[u] == v else 'k'
            nx.draw_networkx_edges(graph, pos, ax=ax, edgelist=[point], alpha=0.3, edge_color=edge_color, width=15,
                                   style='solid')
            nx.draw_networkx_edge_labels(graph, pos, ax=ax, edge_labels={point: cost[u, v]}, font_size=10,
                                         label_pos=0.2)

    # 绘制节点
    nx.draw_networkx_nodes(graph, pos, ax=ax, nodelist=w_nodes, node_color=colors[0], node_size=3000)
    nx.draw_networkx_nodes(graph, pos, ax=ax, nodelist=j_nodes, node_color=colors[1], node_size=3000)

    # 绘制节点标签
    nx.draw_networkx_labels(graph, pos, ax=ax, font_color='w', font_size=20)

    ax.set_title('Minimization Problem', fontsize=20)


def draw(cost, row_ind, col_ind):
    # 创建一个画布，并指定包含两个子图
    fig, axs = plt.subplots(1, 1, figsize=(6, cost.shape[0] * 2), dpi=72)
    for i, ax in np.ndenumerate(axs):
        ax.set_aspect('equal')
        # 显示图形
        ax.axis('off')
    draw_graph(cost, axs, row_ind, col_ind)

    # 显示图形
    plt.axis('off')
    # 调整布局
    plt.tight_layout()
    plt.savefig('figure/km_example.png', bbox_inches='tight', dpi=300)
    plt.show()


if __name__ == '__main__':
    cost = np.array([
        [82, 83, 69, 92],
        [77, 37, 49, 92],
        [11, 69, 5, 86],
        [8, 9, 98, 23],
    ])
    row_ind, col_ind = linear_sum_assignment(cost)
    draw(cost, row_ind, col_ind)

整数规划（MILP）

引入 0-1 决策变量 $x_{ij}$
$x_{ij} = \left\{ \begin{array}{l l} 1, & 如果员工 i 被分配工作 j \\ 0, & 其他 \end{array} \right.$
那么，该指派问题可以被建模为下面的整数规划模型：
$\begin{array}{r l l l} \min & \displaystyle\sum_{i \in W} \displaystyle\sum_{j \in J} x_{ij}c_{ij} \\ {s.t.} & \displaystyle\sum_{j \in J} x_{ij} = 1, & \forall i \in W \\ & \displaystyle\sum_{i \in W} x_{ij} = 1, & \forall j \in J \\ & x_{ij} \in \{0, 1 \}, & \forall j \in J, \forall i \in W \\ \end{array}$