規(guī)劃到任務(wù)分配的最優(yōu)匹配實(shí)戰(zhàn))
1. 項(xiàng)目概述從“分派難題”到匈牙利算法的優(yōu)雅解法在數(shù)學(xué)建模尤其是涉及資源分配、任務(wù)調(diào)度、人員匹配等優(yōu)化問(wèn)題時(shí)我們經(jīng)常會(huì)遇到一類特殊的約束決策變量必須是整數(shù)。這類問(wèn)題就是整數(shù)規(guī)劃。而“匈牙利算法”正是解決其中一類經(jīng)典問(wèn)題——指派問(wèn)題Assignment Problem——的一把利器。它高效、優(yōu)雅且原理直觀是數(shù)學(xué)建模競(jìng)賽和實(shí)際工程中處理最優(yōu)匹配問(wèn)題的必備工具。簡(jiǎn)單來(lái)說(shuō)指派問(wèn)題就是有n項(xiàng)任務(wù)要分配給n個(gè)代理人或機(jī)器每個(gè)代理完成每項(xiàng)任務(wù)的成本或效益已知。如何分配使得總成本最小或總效益最大比如5個(gè)工人操作5臺(tái)機(jī)床每個(gè)工人操作每臺(tái)機(jī)床的效率不同如何安排使總效率最高這就是一個(gè)典型的指派問(wèn)題。匈牙利算法能在多項(xiàng)式時(shí)間內(nèi)為這類問(wèn)題找到一個(gè)最優(yōu)的完美匹配。本文將從一個(gè)建模者的視角深入剖析整數(shù)規(guī)劃中的指派問(wèn)題并詳細(xì)拆解匈牙利算法的核心思想、實(shí)現(xiàn)步驟、代碼實(shí)操以及在實(shí)際建模中可能遇到的變體與陷阱。無(wú)論你是初次接觸數(shù)學(xué)建模的新手還是希望深化對(duì)組合優(yōu)化理解的老手這篇文章都將為你提供從理論到實(shí)戰(zhàn)的完整路徑。2. 整數(shù)規(guī)劃與指派問(wèn)題模型構(gòu)建與核心特征2.1 整數(shù)規(guī)劃的基本框架整數(shù)規(guī)劃是線性規(guī)劃的一個(gè)分支其決策變量被限制為整數(shù)。根據(jù)變量類型可分為純整數(shù)規(guī)劃所有變量為整數(shù)、混合整數(shù)規(guī)劃部分變量為整數(shù)和0-1整數(shù)規(guī)劃變量取0或1。指派問(wèn)題本質(zhì)上就是一個(gè)0-1整數(shù)規(guī)劃問(wèn)題。一個(gè)標(biāo)準(zhǔn)的線性規(guī)劃模型如下目標(biāo)函數(shù)最小化或最大化c^T * x約束條件A * x bx 0當(dāng)對(duì)變量x增加整數(shù)約束x ∈ Z時(shí)就變成了整數(shù)規(guī)劃。整數(shù)約束的引入使得問(wèn)題從連續(xù)的凸優(yōu)化變成了離散的組合優(yōu)化求解難度急劇上升。指派問(wèn)題是其中結(jié)構(gòu)特殊、存在高效專用算法的一類。2.2 指派問(wèn)題的數(shù)學(xué)模型假設(shè)有n個(gè)工人和n項(xiàng)工作c_{ij}表示第i個(gè)工人完成第j項(xiàng)工作的成本。我們引入0-1決策變量x_{ij}x_{ij} 1表示指派工人i去完成工作j。x_{ij} 0表示不指派。那么標(biāo)準(zhǔn)的指派問(wèn)題模型可以表述為目標(biāo)函數(shù)最小化總成本Min Z Σ_{i1}^{n} Σ_{j1}^{n} c_{ij} * x_{ij}約束條件每個(gè)工人只能做一項(xiàng)工作Σ_{j1}^{n} x_{ij} 1, 對(duì)于所有 i 1, 2, ..., n。每項(xiàng)工作只能由一個(gè)工人完成Σ_{i1}^{n} x_{ij} 1, 對(duì)于所有 j 1, 2, ..., n。0-1約束x_{ij} ∈ {0, 1}, 對(duì)于所有 i, j。這個(gè)模型是一個(gè)典型的二分圖完美匹配問(wèn)題。約束矩陣非常特殊全是0和1并且每行每列的和都為1。正是這種特殊的結(jié)構(gòu)使得匈牙利算法能夠繞過(guò)通用的整數(shù)規(guī)劃求解器如分支定界法以更高效的方式找到最優(yōu)解。注意模型默認(rèn)假設(shè)工人數(shù)和工作數(shù)相等即“平衡指派問(wèn)題”。在實(shí)際建模中常會(huì)遇到不相等的情況非平衡問(wèn)題我們會(huì)在后續(xù)章節(jié)討論如何處理。2.3 為什么不用窮舉或通用求解器對(duì)于n5的問(wèn)題所有可能的指派方案有5! 120種窮舉尚可接受。但當(dāng)n10時(shí)10! 3,628,800n15時(shí)15! ≈ 1.3萬(wàn)億。窮舉法顯然不可行。通用的整數(shù)規(guī)劃求解器如Gurobi, CPLEX當(dāng)然可以求解但對(duì)于大規(guī)模的純指派問(wèn)題匈牙利算法的時(shí)間復(fù)雜度為O(n^3)遠(yuǎn)低于通用求解器處理整數(shù)規(guī)劃問(wèn)題的復(fù)雜度。在數(shù)學(xué)建模競(jìng)賽中使用匈牙利算法不僅能保證正確性還能體現(xiàn)你對(duì)問(wèn)題特性和專用算法的掌握是加分項(xiàng)。3. 匈牙利算法核心原理從K?nig定理到增廣路匈牙利算法得名于匈牙利數(shù)學(xué)家Dénes K?nig和Jen? Egerváry的工作。其核心思想是通過(guò)矩陣的變換在不改變最優(yōu)解的前提下逐步“顯露出”一個(gè)完整的、成本為0的完美匹配。3.1 算法的理論基礎(chǔ)K?nig定理與等價(jià)變換算法基于一個(gè)關(guān)鍵原理系數(shù)矩陣的任一行或任一列同時(shí)加上或減去一個(gè)常數(shù)不改變指派問(wèn)題的最優(yōu)解。為什么考慮目標(biāo)函數(shù)Z Σ c_{ij} x_{ij}。如果我們對(duì)第i行所有元素都減去一個(gè)常數(shù)u_i那么新的目標(biāo)函數(shù)為Z Σ (c_{ij} - u_i) x_{ij} Σ c_{ij} x_{ij} - Σ u_i (Σ x_{ij})。由于約束條件Σ x_{ij} 1所以Σ u_i (Σ x_{ij}) Σ u_i是一個(gè)常數(shù)。因此Z和Z只相差一個(gè)常數(shù)它們的最優(yōu)解即x_{ij}的取值完全相同。對(duì)列的操作同理。這個(gè)性質(zhì)允許我們對(duì)成本矩陣進(jìn)行“化簡(jiǎn)”目標(biāo)是讓矩陣中出現(xiàn)盡可能多的零元素并且希望這些零元素的位置能構(gòu)成一個(gè)“獨(dú)立零元素集合”即不同行不同列的零這個(gè)集合就對(duì)應(yīng)著一個(gè)零成本的完美匹配如果存在的話。3.2 算法步驟的直觀理解標(biāo)準(zhǔn)的匈牙利算法通常包含以下幾步我們可以用一個(gè)生活化的類比來(lái)理解假設(shè)成本矩陣是一個(gè)“任務(wù)板”每個(gè)格子c_{ij}是工人i做工作j的“抱怨值”。我們的目標(biāo)是讓總“抱怨”最小。行歸約讓每個(gè)工人對(duì)自己最不擅長(zhǎng)的工作本行最小值的抱怨降為零。即每行減去該行的最小值。這樣每行至少出現(xiàn)一個(gè)零。這相當(dāng)于給每個(gè)工人發(fā)一筆“補(bǔ)貼”消除他們對(duì)自己最討厭工作的基礎(chǔ)抱怨。列歸約行歸約后有些工作可能仍然很“搶手”列中無(wú)零有些則很“冷門”列中有多個(gè)零。我們對(duì)列進(jìn)行同樣操作每列減去該列的最小值。這樣每行每列都至少有一個(gè)零。現(xiàn)在“任務(wù)板”上出現(xiàn)了很多零它們代表“零抱怨”的配對(duì)可能性。試指派與畫線覆蓋我們用最少的水平或垂直線覆蓋住所有的零。為什么這基于組合優(yōu)化中的K?nig定理二分圖中最大匹配數(shù)等于最小點(diǎn)覆蓋數(shù)。在這里“覆蓋所有零的線”對(duì)應(yīng)于點(diǎn)覆蓋。如果最少的線數(shù)等于矩陣的階數(shù)n說(shuō)明我們已經(jīng)找到了n個(gè)位于不同行不同列的零即一個(gè)完美匹配算法結(jié)束這些零的位置就是最優(yōu)指派。如果線數(shù)k n說(shuō)明當(dāng)前的零還不夠“獨(dú)立”無(wú)法直接構(gòu)成完美匹配。我們需要調(diào)整矩陣創(chuàng)造出新的零。矩陣調(diào)整在所有未被線覆蓋的元素中找到最小值min_val。將所有未被線覆蓋的元素減去min_val。將所有被兩條線交叉覆蓋的元素加上min_val。被一條線覆蓋的元素保持不變。 這個(gè)操作的精妙之處在于它保證了原有零元素如果被一條線覆蓋不會(huì)被破壞同時(shí)又在未被覆蓋的區(qū)域創(chuàng)造了新的零。并且它嚴(yán)格遵循了“行/列加減常數(shù)不改變解”的原則因?yàn)閷?duì)未被覆蓋的行或列進(jìn)行了整體減法同時(shí)對(duì)交叉點(diǎn)所在的列或行進(jìn)行了整體加法。重復(fù)迭代回到步驟3用新的矩陣重新畫線覆蓋直到覆蓋線數(shù)等于n為止。3.3 一個(gè)手算示例假設(shè)成本矩陣為工人\工作 | J1 | J2 | J3 ---------|----|----|---- W1 | 2 | 4 | 3 W2 | 5 | 6 | 1 W3 | 3 | 2 | 4步驟1行歸約。每行減最小值W1行減2 W2行減1 W3行減2。 得到0 2 1 4 5 0 1 0 2步驟2列歸約。每列減最小值J1列減0 J2列減0 J3列減0因?yàn)槊苛幸延?。矩陣不變。步驟3畫線覆蓋。嘗試用最少的線覆蓋所有0。先標(biāo)記只有一個(gè)0的列/行。J2列只有一個(gè)0第3行畫線覆蓋第3行。覆蓋后J3列的0第2行未被覆蓋畫線覆蓋J3列。 現(xiàn)在所有0都被覆蓋了用了2條線第3行和J3列。線數(shù)k2 n3。步驟4矩陣調(diào)整。未被覆蓋的元素是(1,1)0,(1,2)2,(2,1)4,(2,2)5。最小值min_val 0實(shí)際上(1,1)的0已被行線覆蓋這里需要仔細(xì)檢查畫線邏輯。讓我們重新規(guī)范地畫線 更系統(tǒng)的方法是先找獨(dú)立0即不同行不同列的0作為初始匹配。假設(shè)我們找到(1,1)0和(3,2)0匹配之。然后發(fā)現(xiàn)工人2無(wú)法匹配到0因?yàn)镴1和J2已被占用。此時(shí)需要用增廣路算法或畫線法。 畫線法對(duì)已匹配的0所在行畫線第1行、第3行。看這些行上的0所在的列第1列、第2列對(duì)這些列畫線。再看這些列上的0所在的行... 最終發(fā)現(xiàn)用線覆蓋第1行、第3行和第1列可以覆蓋所有0。共3條線。等等這不對(duì)線數(shù)不應(yīng)超過(guò)n。這說(shuō)明我的初始匹配沒(méi)找好。 實(shí)際上對(duì)于小矩陣更簡(jiǎn)單的方法是直接觀察。我們發(fā)現(xiàn)可以用兩條線覆蓋所有0覆蓋第3行覆蓋了(3,2)的0和覆蓋第1列覆蓋了(1,1)和(2,3)? 不對(duì)(2,3)不在第1列。(2,3)的0需要被覆蓋。所以嘗試覆蓋第2行和J3列覆蓋第2行覆蓋(2,3)和J3列覆蓋(2,3)和(1,3)(1,3)是1不是0。看來(lái)覆蓋所有0的最小線集是第3行和J3列。是的(3,2)的0被第3行覆蓋(2,3)的0被J3列覆蓋。(1,1)的0呢它沒(méi)有被覆蓋所以我們需要三條線第1行、第3行、J3列。線數(shù)k3等于n3不n3線數(shù)3等于n這意味著我們已經(jīng)找到了完美匹配匹配是(1,1),(2,3),(3,2)。總成本 2 1 2 5。檢查原始矩陣2125。這似乎是一個(gè)可行解。但我們還沒(méi)驗(yàn)證是否最優(yōu)。讓我們用另一種方法(1,1),(2,3),(3,2)。總成本5。有沒(méi)有更優(yōu)的(1,3)3,(2,1)5,(3,2)2總和10。(1,2)4,(2,3)1,(3,1)3總和8。看起來(lái)5確實(shí)是最小的。所以在這個(gè)簡(jiǎn)單例子中步驟2后其實(shí)已經(jīng)得到了最優(yōu)解雖然畫線邏輯有點(diǎn)繞。這個(gè)例子說(shuō)明了算法有時(shí)收斂很快。為了展示調(diào)整步驟我們故意找一個(gè)需要調(diào)整的例子。考慮矩陣3 7 5 4 8 6 5 9 7行歸約后0 4 2 0 4 2 0 4 2列歸約每列減0后不變。現(xiàn)在所有零都在第一列我們無(wú)法找到3個(gè)不同行不同列的零。最少用1條線覆蓋第一列就能蓋住所有零。k1 n3。進(jìn)行調(diào)整未被覆蓋區(qū)域最小值為4。未被覆蓋元素減4交叉點(diǎn)加4。得到新矩陣再迭代。這個(gè)過(guò)程清晰地展示了“創(chuàng)造新零”的過(guò)程。實(shí)操心得手工執(zhí)行匈牙利算法時(shí)畫線找最小覆蓋是最容易出錯(cuò)的一步。對(duì)于競(jìng)賽或編程實(shí)現(xiàn)更推薦使用基于深度優(yōu)先搜索DFS尋找增廣路的算法流程邏輯更清晰更容易編碼。下文將重點(diǎn)介紹這種實(shí)現(xiàn)方式。4. 匈牙利算法的代碼實(shí)現(xiàn)與逐行解析雖然手算有助于理解原理但在數(shù)學(xué)建模中我們幾乎總是通過(guò)編程來(lái)求解。下面以Python為例實(shí)現(xiàn)一個(gè)基于DFS增廣路的匈牙利算法用于求解最小化成本的指派問(wèn)題。4.1 算法核心二分圖最大權(quán)匹配的KM算法 vs. 匈牙利算法這里需要澄清一個(gè)常見混淆。我們通常所說(shuō)的“匈牙利算法”是指求解無(wú)權(quán)二分圖最大匹配的算法。而對(duì)于指派問(wèn)題最小化總成本我們通常使用Kuhn-Munkres算法KM算法它是匈牙利算法在加權(quán)二分圖上的推廣用于求解最大權(quán)完美匹配或最小權(quán)。當(dāng)所有權(quán)重非負(fù)時(shí)通過(guò)將最小化問(wèn)題轉(zhuǎn)化為最大化問(wèn)題例如用一個(gè)大數(shù)減去成本矩陣KM算法可以直接求解。但KM算法復(fù)雜度為O(n^3)且實(shí)現(xiàn)稍復(fù)雜。實(shí)際上對(duì)于最小成本指派問(wèn)題有一個(gè)更直接的轉(zhuǎn)化將成本矩陣的每個(gè)元素取相反數(shù)然后求最大權(quán)匹配。或者使用經(jīng)典的最小成本最大流算法。但還有一種更簡(jiǎn)潔的方式就是直接在我們化簡(jiǎn)后的“零矩陣”上尋找最大匹配即最多的獨(dú)立零元素。當(dāng)找到的匹配數(shù)等于n時(shí)這些零元素的位置就對(duì)應(yīng)著總成本最小的指派因?yàn)榻?jīng)過(guò)變換這些位置的當(dāng)前成本為零而變換不改變最優(yōu)解的結(jié)構(gòu)。下面給出的代碼是求解最小成本指派問(wèn)題的經(jīng)典實(shí)現(xiàn)它融合了矩陣變換歸約和DFS增廣路搜索。import numpy as np class AssignmentProblemSolver: def __init__(self, cost_matrix): 初始化求解器。 :param cost_matrix: 成本矩陣二維numpy數(shù)組shape為(n, n)。 self.n cost_matrix.shape[0] self.original_cost cost_matrix.copy() self.cost cost_matrix.copy().astype(float) # 使用浮點(diǎn)數(shù)以便進(jìn)行減法 # 記錄行、列約減值用于最終還原實(shí)際成本 self.row_reduction np.zeros(self.n) self.col_reduction np.zeros(self.n) # 匹配記錄col_of_row[i] j 表示行i與列j匹配row_of_col[j] i 同理 self.col_of_row -np.ones(self.n, dtypeint) self.row_of_col -np.ones(self.n, dtypeint) # 用于DFS搜索的輔助變量 self.visited_row None self.visited_col None def solve(self): 執(zhí)行匈牙利算法返回最優(yōu)指派和最小總成本。 # 步驟1: 行歸約 for i in range(self.n): min_val np.min(self.cost[i, :]) if min_val 0: # 如果最小值大于0才進(jìn)行歸約 self.cost[i, :] - min_val self.row_reduction[i] min_val # 步驟2: 列歸約 for j in range(self.n): min_val np.min(self.cost[:, j]) if min_val 0: self.cost[:, j] - min_val self.col_reduction[j] min_val # 步驟3: 嘗試尋找初始匹配 (貪心策略) for i in range(self.n): if self.col_of_row[i] -1: # 行i尚未匹配 self._dfs(i) # 如果初始匹配未找到所有匹配則需要進(jìn)入調(diào)整迭代 # 在實(shí)際的完整實(shí)現(xiàn)中這里應(yīng)包含一個(gè)循環(huán)當(dāng)匹配數(shù)小于n時(shí)執(zhí)行矩陣調(diào)整畫線、找最小值、更新矩陣并重新嘗試匹配。 # 為了代碼簡(jiǎn)潔和聚焦核心以下省略了完整的迭代調(diào)整循環(huán)直接假設(shè)初始匹配已成功對(duì)于許多經(jīng)過(guò)歸約的矩陣是成立的。 # 一個(gè)完整的實(shí)現(xiàn)需要包含 _adjust_matrix() 方法和循環(huán)。 # 計(jì)算最小總成本并生成指派方案 total_cost 0.0 assignments [] for i in range(self.n): j self.col_of_row[i] if j ! -1: total_cost self.original_cost[i, j] assignments.append((i, j)) else: # 理論上經(jīng)過(guò)完整算法后不應(yīng)出現(xiàn)未匹配的行 raise RuntimeError(f行 {i} 未找到匹配算法可能未收斂或需要完整迭代。) return assignments, total_cost def _dfs(self, i): 深度優(yōu)先搜索嘗試為行i尋找增廣路。 self.visited_row[i] True for j in range(self.n): if not self.visited_col[j] and abs(self.cost[i, j]) 1e-10: # 判斷是否為0考慮浮點(diǎn)誤差 self.visited_col[j] True # 如果列j未被匹配或者可以為列j的當(dāng)前匹配行找到新的匹配 if self.row_of_col[j] -1 or self._dfs(self.row_of_col[j]): self.col_of_row[i] j self.row_of_col[j] i return True return False # 注意這里省略了完整的 _adjust_matrix() 方法和外層循環(huán)。 # 一個(gè)生產(chǎn)級(jí)的實(shí)現(xiàn)需要它們來(lái)處理所有情況。 # 使用示例 if __name__ __main__: # 示例成本矩陣 cost_matrix np.array([ [2, 4, 3], [5, 6, 1], [3, 2, 4] ]) solver AssignmentProblemSolver(cost_matrix) assignments, min_cost solver.solve() print(最優(yōu)指派方案) for i, j in assignments: print(f 工人{(lán)i1} - 工作{j1} (成本{cost_matrix[i, j]})) print(f最小總成本{min_cost})4.2 代碼關(guān)鍵點(diǎn)解析與避坑指南浮點(diǎn)數(shù)精度問(wèn)題在矩陣變換中反復(fù)的加減可能導(dǎo)致浮點(diǎn)數(shù)誤差。代碼中判斷零時(shí)使用了abs(self.cost[i, j]) 1e-10這是一個(gè)必要的容錯(cuò)處理。在數(shù)學(xué)建模競(jìng)賽中如果成本矩陣是整數(shù)可以全程使用整數(shù)運(yùn)算以避免此問(wèn)題。初始匹配策略上述代碼在行、列歸約后直接對(duì)每一行嘗試DFS匹配。這是一種簡(jiǎn)單的貪心策略。更魯棒的實(shí)現(xiàn)應(yīng)該在DFS失敗后進(jìn)入矩陣調(diào)整階段即前述的手算步驟3和4并循環(huán)直到找到完美匹配。完整的KM算法或最小成本流算法會(huì)系統(tǒng)地處理這個(gè)過(guò)程。算法復(fù)雜度_dfs函數(shù)在最壞情況下會(huì)遍歷所有列并且可能遞歸調(diào)用。如果外層還需要循環(huán)調(diào)整矩陣最壞時(shí)間復(fù)雜度為O(n^4)。通過(guò)優(yōu)化如使用BFS查找增廣路即Hopcroft-Karp算法思想可以將二分圖最大匹配部分優(yōu)化到O(n^2.5)但KM算法的標(biāo)準(zhǔn)實(shí)現(xiàn)是O(n^3)。對(duì)于建模競(jìng)賽中n500的問(wèn)題O(n^3)的實(shí)現(xiàn)完全夠用。使用現(xiàn)成庫(kù)在實(shí)際建模和工程中除非有特殊需求或?qū)W習(xí)目的否則更推薦使用成熟的優(yōu)化庫(kù)。Python:scipy.optimize庫(kù)中的linear_sum_assignment函數(shù)它實(shí)現(xiàn)了高效的匈牙利算法Jonker-Volgenant算法是求解指派問(wèn)題的首選。from scipy.optimize import linear_sum_assignment row_ind, col_ind linear_sum_assignment(cost_matrix) min_cost cost_matrix[row_ind, col_ind].sum()MATLAB:assign函數(shù)或matchpairs函數(shù)。Lingo/LINDO: 直接建立整數(shù)規(guī)劃模型求解。重要提示在數(shù)學(xué)建模論文中如果你使用了scipy.optimize.linear_sum_assignment你仍然需要清晰地闡述匈牙利算法的基本原理。你可以寫“針對(duì)該指派問(wèn)題我們采用經(jīng)典的匈牙利算法進(jìn)行求解。在具體實(shí)現(xiàn)上我們調(diào)用了SciPy庫(kù)中的linear_sum_assignment函數(shù)該函數(shù)基于高效的Jonker-Volgenant算法能夠在多項(xiàng)式時(shí)間內(nèi)保證找到全局最優(yōu)解。” 這既體現(xiàn)了你對(duì)算法的理解也展示了你會(huì)利用高效工具。5. 數(shù)學(xué)建模中的實(shí)戰(zhàn)應(yīng)用與變體處理匈牙利算法不僅僅是解教科書上的標(biāo)準(zhǔn)問(wèn)題。在數(shù)學(xué)建模競(jìng)賽中問(wèn)題往往披著各種“外衣”需要你識(shí)別并轉(zhuǎn)化為指派問(wèn)題模型。5.1 經(jīng)典應(yīng)用場(chǎng)景識(shí)別任務(wù)分配這是最直接的應(yīng)用。如論文評(píng)審分配每位評(píng)審審閱幾篇論文總匹配度最高、出租車派單車與乘客的距離最小、教室安排課程與教室的適配度。路徑規(guī)劃與排序某些旅行商問(wèn)題TSP的近似解法中會(huì)用到指派問(wèn)題來(lái)構(gòu)建匹配。例如將城市兩兩配對(duì)然后連接這些配對(duì)形成路徑。資源調(diào)度在固定時(shí)間段內(nèi)將機(jī)器分配給加工任務(wù)使得總加工時(shí)間最短或利潤(rùn)最大。圖像處理與數(shù)據(jù)關(guān)聯(lián)在多目標(biāo)跟蹤中將上一幀的檢測(cè)框與當(dāng)前幀的檢測(cè)框進(jìn)行關(guān)聯(lián)關(guān)聯(lián)成本可以是邊界框的重疊度IoU的負(fù)數(shù)。平衡實(shí)驗(yàn)設(shè)計(jì)將實(shí)驗(yàn)對(duì)象如患者分配到不同的實(shí)驗(yàn)組和控制組使得各組在某些特征上盡可能平衡這可以轉(zhuǎn)化為一個(gè)最小化組間差異的指派問(wèn)題。5.2 非標(biāo)準(zhǔn)情況的處理技巧1. 非平衡指派問(wèn)題工人數(shù) ≠ 工作數(shù)工人多工作少引入“虛擬工作”其成本設(shè)為0如果是最小化問(wèn)題。這意味著多余的工人沒(méi)有被指派任務(wù)成本為0。工人少工作多引入“虛擬工人”其完成所有工作的成本設(shè)為0。這意味著多余的工作沒(méi)有被完成成本為0。注意虛擬行/列的成本設(shè)置取決于問(wèn)題目標(biāo)。如果是最大化效益問(wèn)題虛擬行/列的效益通常設(shè)為0或一個(gè)非常大的負(fù)數(shù)在最大化問(wèn)題中表示不選擇。2. 最大化問(wèn)題標(biāo)準(zhǔn)匈牙利算法解決最小化問(wèn)題。對(duì)于最大化問(wèn)題如最大效益、最大匹配度常用方法有方法一將效益矩陣B轉(zhuǎn)化為成本矩陣C M - B其中M是矩陣B中元素的最大值或一個(gè)足夠大的數(shù)。然后對(duì)C求解最小化指派。因?yàn)镸in Σ(M - b_{ij})x_{ij} M*n - Max Σ b_{ij}x_{ij}所以解相同。方法二直接對(duì)效益矩陣B取負(fù)值C -B然后求解最小化指派。3. 禁止指派某些工人不能做某些工作。處理方法是將對(duì)應(yīng)成本設(shè)為一個(gè)極大的數(shù)INF。在最小化問(wèn)題中算法會(huì)主動(dòng)避免選擇成本為INF的配對(duì)。在代碼實(shí)現(xiàn)中可以用一個(gè)遠(yuǎn)大于其他正常成本的值如1e9來(lái)代替INF。4. 多對(duì)一或一對(duì)多指派標(biāo)準(zhǔn)指派是一對(duì)一。如果允許一個(gè)工人做多項(xiàng)工作或者一項(xiàng)工作需要多個(gè)工人問(wèn)題就變成了廣義分配問(wèn)題Generalized Assignment Problem, GAP這比標(biāo)準(zhǔn)指派問(wèn)題復(fù)雜得多通常需要用到更高級(jí)的整數(shù)規(guī)劃或啟發(fā)式算法如遺傳算法、模擬退火。此時(shí)匈牙利算法不再直接適用。5. 有額外約束的指派例如除了成本最小還要求某些工人必須被分配到一起或者某些工作必須在其他工作之后完成。這些約束破壞了二分圖匹配的結(jié)構(gòu)需要將其建模為更復(fù)雜的整數(shù)規(guī)劃問(wèn)題使用通用求解器如Gurobi, CPLEX或定制算法。5.3 建模實(shí)例數(shù)學(xué)建模競(jìng)賽題改編問(wèn)題某市有5個(gè)突發(fā)公共事件應(yīng)急點(diǎn)現(xiàn)有5支救援隊(duì)。已知各救援隊(duì)到達(dá)各應(yīng)急點(diǎn)的預(yù)計(jì)時(shí)間小時(shí)。由于專業(yè)設(shè)備限制第2支救援隊(duì)無(wú)法前往第3個(gè)應(yīng)急點(diǎn)。如何分配救援隊(duì)使得總響應(yīng)時(shí)間最短建模步驟定義決策變量x_{ij} 1表示派遣救援隊(duì)i到應(yīng)急點(diǎn)j否則為0。建立成本矩陣c_{ij}為行駛時(shí)間。對(duì)于禁止指派救援隊(duì)2 - 應(yīng)急點(diǎn)3令c_{23} INF一個(gè)大數(shù)如999。目標(biāo)函數(shù)Min Σ Σ c_{ij} x_{ij}。約束條件標(biāo)準(zhǔn)的指派問(wèn)題約束每行每列和為1。求解使用匈牙利算法或scipy.optimize.linear_sum_assignment求解。結(jié)果分析檢查最優(yōu)解中x_{23}是否為0驗(yàn)證禁止指派是否被遵守。計(jì)算總響應(yīng)時(shí)間。實(shí)操心得在論文寫作中將原始問(wèn)題抽象成矩陣形式是關(guān)鍵一步。建議在論文中清晰地畫出成本矩陣表格并對(duì)特殊值如INF加以說(shuō)明。這能讓評(píng)委一眼看出你正確理解了問(wèn)題并進(jìn)行了恰當(dāng)?shù)霓D(zhuǎn)化。6. 常見問(wèn)題、調(diào)試技巧與算法局限6.1 算法實(shí)現(xiàn)中的常見陷阱浮點(diǎn)誤差導(dǎo)致匹配失敗如前所述在判斷c_{ij} 0時(shí)使用絕對(duì)容差abs(cost[i][j]) 1e-10而不是cost[i][j] 0。非方陣處理不當(dāng)對(duì)于非平衡問(wèn)題務(wù)必先將其補(bǔ)全為方陣并正確設(shè)置虛擬行/列的成本。如果目標(biāo)是最大化補(bǔ)全時(shí)需要格外小心。無(wú)限循環(huán)在自編的完整迭代算法中如果矩陣調(diào)整步驟的邏輯有誤可能導(dǎo)致無(wú)法增加匹配數(shù)從而陷入無(wú)限循環(huán)。確保每次調(diào)整后至少有一個(gè)新的零元素在未被覆蓋的區(qū)域產(chǎn)生。誤用最大化算法直接將最大化問(wèn)題的矩陣輸入給最小化算法會(huì)得到錯(cuò)誤結(jié)果。務(wù)必先進(jìn)行轉(zhuǎn)化。6.2 調(diào)試與驗(yàn)證小規(guī)模驗(yàn)證用3x3或4x4的矩陣手動(dòng)計(jì)算與程序結(jié)果對(duì)比。這是最有效的調(diào)試方法。檢查解的可行性確保得到的指派方案滿足“每個(gè)代理恰好一個(gè)任務(wù)”的約束。計(jì)算row_ind和col_ind是否都是[0, 1, ..., n-1]的一個(gè)排列。與暴力枚舉對(duì)比對(duì)于n很小如n8的問(wèn)題可以編寫暴力枚舉所有排列的程序驗(yàn)證匈牙利算法給出的解是否確實(shí)是最優(yōu)的。使用庫(kù)函數(shù)交叉驗(yàn)證用scipy.optimize.linear_sum_assignment的結(jié)果來(lái)驗(yàn)證自己編寫的算法。6.3 匈牙利算法的局限與替代方案盡管匈牙利算法高效但它有其適用范圍僅適用于線性目標(biāo)函數(shù)總成本必須是各配對(duì)成本的和。如果是非線性如成本與配對(duì)順序有關(guān)則不適用。一對(duì)一嚴(yán)格約束這是核心假設(shè)。一對(duì)多、多對(duì)多需要其他模型。單目標(biāo)優(yōu)化只能處理最小化總成本或最大化總效益。多目標(biāo)指派問(wèn)題需要其他方法如目標(biāo)規(guī)劃、進(jìn)化算法。替代算法拍賣算法Auction Algorithm另一種求解指派問(wèn)題的經(jīng)典算法思想直觀模擬拍賣過(guò)程在某些情況下并行性好。最小成本最大流將指派問(wèn)題建模為網(wǎng)絡(luò)流問(wèn)題。源點(diǎn)連接所有工人容量1成本0工人連接所有工作容量1成本為c_ij工作連接匯點(diǎn)容量1成本0。求解從源到匯的最小成本最大流。這是一個(gè)更通用的框架可以處理更多變體。整數(shù)規(guī)劃求解器對(duì)于復(fù)雜約束的指派問(wèn)題直接使用Gurobi、CPLEX等求解器建模求解是最穩(wěn)妥的方式。雖然可能不如專用算法快但能保證在復(fù)雜約束下找到最優(yōu)解如果問(wèn)題可解。在我多年的建模和編程經(jīng)驗(yàn)中處理指派問(wèn)題的首選路徑是首先判斷是否是標(biāo)準(zhǔn)的一對(duì)一、線性成本問(wèn)題。如果是毫不猶豫地使用scipy.optimize.linear_sum_assignment。如果問(wèn)題帶有特殊約束如資源容量、先后順序則將其建立為整數(shù)規(guī)劃模型調(diào)用專業(yè)求解器。匈牙利算法的價(jià)值在于其優(yōu)美的理論和作為構(gòu)建更復(fù)雜算法基礎(chǔ)組件的作用但在實(shí)際應(yīng)用中我們更應(yīng)注重正確、高效地解決問(wèn)題而非重復(fù)造輪子。理解匈牙利算法就像是掌握了一把打開組合優(yōu)化大門的鑰匙。它讓你看到對(duì)于具有特殊結(jié)構(gòu)的整數(shù)規(guī)劃問(wèn)題存在比蠻力搜索和通用求解器更巧妙的道路。這種“發(fā)現(xiàn)結(jié)構(gòu)、利用結(jié)構(gòu)”的思維才是數(shù)學(xué)建模中最寶貴的財(cái)富。當(dāng)你下次遇到分配、匹配、調(diào)度類問(wèn)題時(shí)不妨先想一想這能不能抽象成一個(gè)二分圖能不能用匈牙利算法的思想來(lái)近似或求解這種思考習(xí)慣往往能讓你在競(jìng)賽中脫穎而出。 匈牙利算法：從整數(shù)規(guī)劃到任務(wù)分配的最優(yōu)匹配實(shí)戰(zhàn)
1. 項(xiàng)目概述從“分派難題”到匈牙利算法的優(yōu)雅解法在數(shù)學(xué)建模尤其是涉及資源分配、任務(wù)調(diào)度、人員匹配等優(yōu)化問(wèn)題時(shí)我們經(jīng)常會(huì)遇到一類特殊的約束決策變量必須是整數(shù)。這類問(wèn)題就是整數(shù)規(guī)劃。而“匈牙利算法”正是解決其中一類經(jīng)典問(wèn)題——指派問(wèn)題Assignment Problem——的一把利器。它高效、優(yōu)雅且原理直觀是數(shù)學(xué)建模競(jìng)賽和實(shí)際工程中處理最優(yōu)匹配問(wèn)題的必備工具。簡(jiǎn)單來(lái)說(shuō)指派問(wèn)題就是有n項(xiàng)任務(wù)要分配給n個(gè)代理人或機(jī)器每個(gè)代理完成每項(xiàng)任務(wù)的成本或效益已知。如何分配使得總成本最小或總效益最大比如5個(gè)工人操作5臺(tái)機(jī)床每個(gè)工人操作每臺(tái)機(jī)床的效率不同如何安排使總效率最高這就是一個(gè)典型的指派問(wèn)題。匈牙利算法能在多項(xiàng)式時(shí)間內(nèi)為這類問(wèn)題找到一個(gè)最優(yōu)的完美匹配。本文將從一個(gè)建模者的視角深入剖析整數(shù)規(guī)劃中的指派問(wèn)題并詳細(xì)拆解匈牙利算法的核心思想、實(shí)現(xiàn)步驟、代碼實(shí)操以及在實(shí)際建模中可能遇到的變體與陷阱。無(wú)論你是初次接觸數(shù)學(xué)建模的新手還是希望深化對(duì)組合優(yōu)化理解的老手這篇文章都將為你提供從理論到實(shí)戰(zhàn)的完整路徑。2. 整數(shù)規(guī)劃與指派問(wèn)題模型構(gòu)建與核心特征2.1 整數(shù)規(guī)劃的基本框架整數(shù)規(guī)劃是線性規(guī)劃的一個(gè)分支其決策變量被限制為整數(shù)。根據(jù)變量類型可分為純整數(shù)規(guī)劃所有變量為整數(shù)、混合整數(shù)規(guī)劃部分變量為整數(shù)和0-1整數(shù)規(guī)劃變量取0或1。指派問(wèn)題本質(zhì)上就是一個(gè)0-1整數(shù)規(guī)劃問(wèn)題。一個(gè)標(biāo)準(zhǔn)的線性規(guī)劃模型如下目標(biāo)函數(shù)最小化或最大化c^T * x約束條件A * x bx 0當(dāng)對(duì)變量x增加整數(shù)約束x ∈ Z時(shí)就變成了整數(shù)規(guī)劃。整數(shù)約束的引入使得問(wèn)題從連續(xù)的凸優(yōu)化變成了離散的組合優(yōu)化求解難度急劇上升。指派問(wèn)題是其中結(jié)構(gòu)特殊、存在高效專用算法的一類。2.2 指派問(wèn)題的數(shù)學(xué)模型假設(shè)有n個(gè)工人和n項(xiàng)工作c_{ij}表示第i個(gè)工人完成第j項(xiàng)工作的成本。我們引入0-1決策變量x_{ij}x_{ij} 1表示指派工人i去完成工作j。x_{ij} 0表示不指派。那么標(biāo)準(zhǔn)的指派問(wèn)題模型可以表述為目標(biāo)函數(shù)最小化總成本Min Z Σ_{i1}^{n} Σ_{j1}^{n} c_{ij} * x_{ij}約束條件每個(gè)工人只能做一項(xiàng)工作Σ_{j1}^{n} x_{ij} 1, 對(duì)于所有 i 1, 2, ..., n。每項(xiàng)工作只能由一個(gè)工人完成Σ_{i1}^{n} x_{ij} 1, 對(duì)于所有 j 1, 2, ..., n。0-1約束x_{ij} ∈ {0, 1}, 對(duì)于所有 i, j。這個(gè)模型是一個(gè)典型的二分圖完美匹配問(wèn)題。約束矩陣非常特殊全是0和1并且每行每列的和都為1。正是這種特殊的結(jié)構(gòu)使得匈牙利算法能夠繞過(guò)通用的整數(shù)規(guī)劃求解器如分支定界法以更高效的方式找到最優(yōu)解。注意模型默認(rèn)假設(shè)工人數(shù)和工作數(shù)相等即“平衡指派問(wèn)題”。在實(shí)際建模中常會(huì)遇到不相等的情況非平衡問(wèn)題我們會(huì)在后續(xù)章節(jié)討論如何處理。2.3 為什么不用窮舉或通用求解器對(duì)于n5的問(wèn)題所有可能的指派方案有5! 120種窮舉尚可接受。但當(dāng)n10時(shí)10! 3,628,800n15時(shí)15! ≈ 1.3萬(wàn)億。窮舉法顯然不可行。通用的整數(shù)規(guī)劃求解器如Gurobi, CPLEX當(dāng)然可以求解但對(duì)于大規(guī)模的純指派問(wèn)題匈牙利算法的時(shí)間復(fù)雜度為O(n^3)遠(yuǎn)低于通用求解器處理整數(shù)規(guī)劃問(wèn)題的復(fù)雜度。在數(shù)學(xué)建模競(jìng)賽中使用匈牙利算法不僅能保證正確性還能體現(xiàn)你對(duì)問(wèn)題特性和專用算法的掌握是加分項(xiàng)。3. 匈牙利算法核心原理從K?nig定理到增廣路匈牙利算法得名于匈牙利數(shù)學(xué)家Dénes K?nig和Jen? Egerváry的工作。其核心思想是通過(guò)矩陣的變換在不改變最優(yōu)解的前提下逐步“顯露出”一個(gè)完整的、成本為0的完美匹配。3.1 算法的理論基礎(chǔ)K?nig定理與等價(jià)變換算法基于一個(gè)關(guān)鍵原理系數(shù)矩陣的任一行或任一列同時(shí)加上或減去一個(gè)常數(shù)不改變指派問(wèn)題的最優(yōu)解。為什么考慮目標(biāo)函數(shù)Z Σ c_{ij} x_{ij}。如果我們對(duì)第i行所有元素都減去一個(gè)常數(shù)u_i那么新的目標(biāo)函數(shù)為Z Σ (c_{ij} - u_i) x_{ij} Σ c_{ij} x_{ij} - Σ u_i (Σ x_{ij})。由于約束條件Σ x_{ij} 1所以Σ u_i (Σ x_{ij}) Σ u_i是一個(gè)常數(shù)。因此Z和Z只相差一個(gè)常數(shù)它們的最優(yōu)解即x_{ij}的取值完全相同。對(duì)列的操作同理。這個(gè)性質(zhì)允許我們對(duì)成本矩陣進(jìn)行“化簡(jiǎn)”目標(biāo)是讓矩陣中出現(xiàn)盡可能多的零元素并且希望這些零元素的位置能構(gòu)成一個(gè)“獨(dú)立零元素集合”即不同行不同列的零這個(gè)集合就對(duì)應(yīng)著一個(gè)零成本的完美匹配如果存在的話。3.2 算法步驟的直觀理解標(biāo)準(zhǔn)的匈牙利算法通常包含以下幾步我們可以用一個(gè)生活化的類比來(lái)理解假設(shè)成本矩陣是一個(gè)“任務(wù)板”每個(gè)格子c_{ij}是工人i做工作j的“抱怨值”。我們的目標(biāo)是讓總“抱怨”最小。行歸約讓每個(gè)工人對(duì)自己最不擅長(zhǎng)的工作本行最小值的抱怨降為零。即每行減去該行的最小值。這樣每行至少出現(xiàn)一個(gè)零。這相當(dāng)于給每個(gè)工人發(fā)一筆“補(bǔ)貼”消除他們對(duì)自己最討厭工作的基礎(chǔ)抱怨。列歸約行歸約后有些工作可能仍然很“搶手”列中無(wú)零有些則很“冷門”列中有多個(gè)零。我們對(duì)列進(jìn)行同樣操作每列減去該列的最小值。這樣每行每列都至少有一個(gè)零。現(xiàn)在“任務(wù)板”上出現(xiàn)了很多零它們代表“零抱怨”的配對(duì)可能性。試指派與畫線覆蓋我們用最少的水平或垂直線覆蓋住所有的零。為什么這基于組合優(yōu)化中的K?nig定理二分圖中最大匹配數(shù)等于最小點(diǎn)覆蓋數(shù)。在這里“覆蓋所有零的線”對(duì)應(yīng)于點(diǎn)覆蓋。如果最少的線數(shù)等于矩陣的階數(shù)n說(shuō)明我們已經(jīng)找到了n個(gè)位于不同行不同列的零即一個(gè)完美匹配算法結(jié)束這些零的位置就是最優(yōu)指派。如果線數(shù)k n說(shuō)明當(dāng)前的零還不夠“獨(dú)立”無(wú)法直接構(gòu)成完美匹配。我們需要調(diào)整矩陣創(chuàng)造出新的零。矩陣調(diào)整在所有未被線覆蓋的元素中找到最小值min_val。將所有未被線覆蓋的元素減去min_val。將所有被兩條線交叉覆蓋的元素加上min_val。被一條線覆蓋的元素保持不變。 這個(gè)操作的精妙之處在于它保證了原有零元素如果被一條線覆蓋不會(huì)被破壞同時(shí)又在未被覆蓋的區(qū)域創(chuàng)造了新的零。并且它嚴(yán)格遵循了“行/列加減常數(shù)不改變解”的原則因?yàn)閷?duì)未被覆蓋的行或列進(jìn)行了整體減法同時(shí)對(duì)交叉點(diǎn)所在的列或行進(jìn)行了整體加法。重復(fù)迭代回到步驟3用新的矩陣重新畫線覆蓋直到覆蓋線數(shù)等于n為止。3.3 一個(gè)手算示例假設(shè)成本矩陣為工人\工作 | J1 | J2 | J3 ---------|----|----|---- W1 | 2 | 4 | 3 W2 | 5 | 6 | 1 W3 | 3 | 2 | 4步驟1行歸約。每行減最小值W1行減2 W2行減1 W3行減2。 得到0 2 1 4 5 0 1 0 2步驟2列歸約。每列減最小值J1列減0 J2列減0 J3列減0因?yàn)槊苛幸延?。矩陣不變。步驟3畫線覆蓋。嘗試用最少的線覆蓋所有0。先標(biāo)記只有一個(gè)0的列/行。J2列只有一個(gè)0第3行畫線覆蓋第3行。覆蓋后J3列的0第2行未被覆蓋畫線覆蓋J3列。 現(xiàn)在所有0都被覆蓋了用了2條線第3行和J3列。線數(shù)k2 n3。步驟4矩陣調(diào)整。未被覆蓋的元素是(1,1)0,(1,2)2,(2,1)4,(2,2)5。最小值min_val 0實(shí)際上(1,1)的0已被行線覆蓋這里需要仔細(xì)檢查畫線邏輯。讓我們重新規(guī)范地畫線 更系統(tǒng)的方法是先找獨(dú)立0即不同行不同列的0作為初始匹配。假設(shè)我們找到(1,1)0和(3,2)0匹配之。然后發(fā)現(xiàn)工人2無(wú)法匹配到0因?yàn)镴1和J2已被占用。此時(shí)需要用增廣路算法或畫線法。 畫線法對(duì)已匹配的0所在行畫線第1行、第3行。看這些行上的0所在的列第1列、第2列對(duì)這些列畫線。再看這些列上的0所在的行... 最終發(fā)現(xiàn)用線覆蓋第1行、第3行和第1列可以覆蓋所有0。共3條線。等等這不對(duì)線數(shù)不應(yīng)超過(guò)n。這說(shuō)明我的初始匹配沒(méi)找好。 實(shí)際上對(duì)于小矩陣更簡(jiǎn)單的方法是直接觀察。我們發(fā)現(xiàn)可以用兩條線覆蓋所有0覆蓋第3行覆蓋了(3,2)的0和覆蓋第1列覆蓋了(1,1)和(2,3)? 不對(duì)(2,3)不在第1列。(2,3)的0需要被覆蓋。所以嘗試覆蓋第2行和J3列覆蓋第2行覆蓋(2,3)和J3列覆蓋(2,3)和(1,3)(1,3)是1不是0。看來(lái)覆蓋所有0的最小線集是第3行和J3列。是的(3,2)的0被第3行覆蓋(2,3)的0被J3列覆蓋。(1,1)的0呢它沒(méi)有被覆蓋所以我們需要三條線第1行、第3行、J3列。線數(shù)k3等于n3不n3線數(shù)3等于n這意味著我們已經(jīng)找到了完美匹配匹配是(1,1),(2,3),(3,2)。總成本 2 1 2 5。檢查原始矩陣2125。這似乎是一個(gè)可行解。但我們還沒(méi)驗(yàn)證是否最優(yōu)。讓我們用另一種方法(1,1),(2,3),(3,2)。總成本5。有沒(méi)有更優(yōu)的(1,3)3,(2,1)5,(3,2)2總和10。(1,2)4,(2,3)1,(3,1)3總和8。看起來(lái)5確實(shí)是最小的。所以在這個(gè)簡(jiǎn)單例子中步驟2后其實(shí)已經(jīng)得到了最優(yōu)解雖然畫線邏輯有點(diǎn)繞。這個(gè)例子說(shuō)明了算法有時(shí)收斂很快。為了展示調(diào)整步驟我們故意找一個(gè)需要調(diào)整的例子。考慮矩陣3 7 5 4 8 6 5 9 7行歸約后0 4 2 0 4 2 0 4 2列歸約每列減0后不變。現(xiàn)在所有零都在第一列我們無(wú)法找到3個(gè)不同行不同列的零。最少用1條線覆蓋第一列就能蓋住所有零。k1 n3。進(jìn)行調(diào)整未被覆蓋區(qū)域最小值為4。未被覆蓋元素減4交叉點(diǎn)加4。得到新矩陣再迭代。這個(gè)過(guò)程清晰地展示了“創(chuàng)造新零”的過(guò)程。實(shí)操心得手工執(zhí)行匈牙利算法時(shí)畫線找最小覆蓋是最容易出錯(cuò)的一步。對(duì)于競(jìng)賽或編程實(shí)現(xiàn)更推薦使用基于深度優(yōu)先搜索DFS尋找增廣路的算法流程邏輯更清晰更容易編碼。下文將重點(diǎn)介紹這種實(shí)現(xiàn)方式。4. 匈牙利算法的代碼實(shí)現(xiàn)與逐行解析雖然手算有助于理解原理但在數(shù)學(xué)建模中我們幾乎總是通過(guò)編程來(lái)求解。下面以Python為例實(shí)現(xiàn)一個(gè)基于DFS增廣路的匈牙利算法用于求解最小化成本的指派問(wèn)題。4.1 算法核心二分圖最大權(quán)匹配的KM算法 vs. 匈牙利算法這里需要澄清一個(gè)常見混淆。我們通常所說(shuō)的“匈牙利算法”是指求解無(wú)權(quán)二分圖最大匹配的算法。而對(duì)于指派問(wèn)題最小化總成本我們通常使用Kuhn-Munkres算法KM算法它是匈牙利算法在加權(quán)二分圖上的推廣用于求解最大權(quán)完美匹配或最小權(quán)。當(dāng)所有權(quán)重非負(fù)時(shí)通過(guò)將最小化問(wèn)題轉(zhuǎn)化為最大化問(wèn)題例如用一個(gè)大數(shù)減去成本矩陣KM算法可以直接求解。但KM算法復(fù)雜度為O(n^3)且實(shí)現(xiàn)稍復(fù)雜。實(shí)際上對(duì)于最小成本指派問(wèn)題有一個(gè)更直接的轉(zhuǎn)化將成本矩陣的每個(gè)元素取相反數(shù)然后求最大權(quán)匹配。或者使用經(jīng)典的最小成本最大流算法。但還有一種更簡(jiǎn)潔的方式就是直接在我們化簡(jiǎn)后的“零矩陣”上尋找最大匹配即最多的獨(dú)立零元素。當(dāng)找到的匹配數(shù)等于n時(shí)這些零元素的位置就對(duì)應(yīng)著總成本最小的指派因?yàn)榻?jīng)過(guò)變換這些位置的當(dāng)前成本為零而變換不改變最優(yōu)解的結(jié)構(gòu)。下面給出的代碼是求解最小成本指派問(wèn)題的經(jīng)典實(shí)現(xiàn)它融合了矩陣變換歸約和DFS增廣路搜索。import numpy as np class AssignmentProblemSolver: def __init__(self, cost_matrix): 初始化求解器。 :param cost_matrix: 成本矩陣二維numpy數(shù)組shape為(n, n)。 self.n cost_matrix.shape[0] self.original_cost cost_matrix.copy() self.cost cost_matrix.copy().astype(float) # 使用浮點(diǎn)數(shù)以便進(jìn)行減法 # 記錄行、列約減值用于最終還原實(shí)際成本 self.row_reduction np.zeros(self.n) self.col_reduction np.zeros(self.n) # 匹配記錄col_of_row[i] j 表示行i與列j匹配row_of_col[j] i 同理 self.col_of_row -np.ones(self.n, dtypeint) self.row_of_col -np.ones(self.n, dtypeint) # 用于DFS搜索的輔助變量 self.visited_row None self.visited_col None def solve(self): 執(zhí)行匈牙利算法返回最優(yōu)指派和最小總成本。 # 步驟1: 行歸約 for i in range(self.n): min_val np.min(self.cost[i, :]) if min_val 0: # 如果最小值大于0才進(jìn)行歸約 self.cost[i, :] - min_val self.row_reduction[i] min_val # 步驟2: 列歸約 for j in range(self.n): min_val np.min(self.cost[:, j]) if min_val 0: self.cost[:, j] - min_val self.col_reduction[j] min_val # 步驟3: 嘗試尋找初始匹配 (貪心策略) for i in range(self.n): if self.col_of_row[i] -1: # 行i尚未匹配 self._dfs(i) # 如果初始匹配未找到所有匹配則需要進(jìn)入調(diào)整迭代 # 在實(shí)際的完整實(shí)現(xiàn)中這里應(yīng)包含一個(gè)循環(huán)當(dāng)匹配數(shù)小于n時(shí)執(zhí)行矩陣調(diào)整畫線、找最小值、更新矩陣并重新嘗試匹配。 # 為了代碼簡(jiǎn)潔和聚焦核心以下省略了完整的迭代調(diào)整循環(huán)直接假設(shè)初始匹配已成功對(duì)于許多經(jīng)過(guò)歸約的矩陣是成立的。 # 一個(gè)完整的實(shí)現(xiàn)需要包含 _adjust_matrix() 方法和循環(huán)。 # 計(jì)算最小總成本并生成指派方案 total_cost 0.0 assignments [] for i in range(self.n): j self.col_of_row[i] if j ! -1: total_cost self.original_cost[i, j] assignments.append((i, j)) else: # 理論上經(jīng)過(guò)完整算法后不應(yīng)出現(xiàn)未匹配的行 raise RuntimeError(f行 {i} 未找到匹配算法可能未收斂或需要完整迭代。) return assignments, total_cost def _dfs(self, i): 深度優(yōu)先搜索嘗試為行i尋找增廣路。 self.visited_row[i] True for j in range(self.n): if not self.visited_col[j] and abs(self.cost[i, j]) 1e-10: # 判斷是否為0考慮浮點(diǎn)誤差 self.visited_col[j] True # 如果列j未被匹配或者可以為列j的當(dāng)前匹配行找到新的匹配 if self.row_of_col[j] -1 or self._dfs(self.row_of_col[j]): self.col_of_row[i] j self.row_of_col[j] i return True return False # 注意這里省略了完整的 _adjust_matrix() 方法和外層循環(huán)。 # 一個(gè)生產(chǎn)級(jí)的實(shí)現(xiàn)需要它們來(lái)處理所有情況。 # 使用示例 if __name__ __main__: # 示例成本矩陣 cost_matrix np.array([ [2, 4, 3], [5, 6, 1], [3, 2, 4] ]) solver AssignmentProblemSolver(cost_matrix) assignments, min_cost solver.solve() print(最優(yōu)指派方案) for i, j in assignments: print(f 工人{(lán)i1} - 工作{j1} (成本{cost_matrix[i, j]})) print(f最小總成本{min_cost})4.2 代碼關(guān)鍵點(diǎn)解析與避坑指南浮點(diǎn)數(shù)精度問(wèn)題在矩陣變換中反復(fù)的加減可能導(dǎo)致浮點(diǎn)數(shù)誤差。代碼中判斷零時(shí)使用了abs(self.cost[i, j]) 1e-10這是一個(gè)必要的容錯(cuò)處理。在數(shù)學(xué)建模競(jìng)賽中如果成本矩陣是整數(shù)可以全程使用整數(shù)運(yùn)算以避免此問(wèn)題。初始匹配策略上述代碼在行、列歸約后直接對(duì)每一行嘗試DFS匹配。這是一種簡(jiǎn)單的貪心策略。更魯棒的實(shí)現(xiàn)應(yīng)該在DFS失敗后進(jìn)入矩陣調(diào)整階段即前述的手算步驟3和4并循環(huán)直到找到完美匹配。完整的KM算法或最小成本流算法會(huì)系統(tǒng)地處理這個(gè)過(guò)程。算法復(fù)雜度_dfs函數(shù)在最壞情況下會(huì)遍歷所有列并且可能遞歸調(diào)用。如果外層還需要循環(huán)調(diào)整矩陣最壞時(shí)間復(fù)雜度為O(n^4)。通過(guò)優(yōu)化如使用BFS查找增廣路即Hopcroft-Karp算法思想可以將二分圖最大匹配部分優(yōu)化到O(n^2.5)但KM算法的標(biāo)準(zhǔn)實(shí)現(xiàn)是O(n^3)。對(duì)于建模競(jìng)賽中n500的問(wèn)題O(n^3)的實(shí)現(xiàn)完全夠用。使用現(xiàn)成庫(kù)在實(shí)際建模和工程中除非有特殊需求或?qū)W習(xí)目的否則更推薦使用成熟的優(yōu)化庫(kù)。Python:scipy.optimize庫(kù)中的linear_sum_assignment函數(shù)它實(shí)現(xiàn)了高效的匈牙利算法Jonker-Volgenant算法是求解指派問(wèn)題的首選。from scipy.optimize import linear_sum_assignment row_ind, col_ind linear_sum_assignment(cost_matrix) min_cost cost_matrix[row_ind, col_ind].sum()MATLAB:assign函數(shù)或matchpairs函數(shù)。Lingo/LINDO: 直接建立整數(shù)規(guī)劃模型求解。重要提示在數(shù)學(xué)建模論文中如果你使用了scipy.optimize.linear_sum_assignment你仍然需要清晰地闡述匈牙利算法的基本原理。你可以寫“針對(duì)該指派問(wèn)題我們采用經(jīng)典的匈牙利算法進(jìn)行求解。在具體實(shí)現(xiàn)上我們調(diào)用了SciPy庫(kù)中的linear_sum_assignment函數(shù)該函數(shù)基于高效的Jonker-Volgenant算法能夠在多項(xiàng)式時(shí)間內(nèi)保證找到全局最優(yōu)解。” 這既體現(xiàn)了你對(duì)算法的理解也展示了你會(huì)利用高效工具。5. 數(shù)學(xué)建模中的實(shí)戰(zhàn)應(yīng)用與變體處理匈牙利算法不僅僅是解教科書上的標(biāo)準(zhǔn)問(wèn)題。在數(shù)學(xué)建模競(jìng)賽中問(wèn)題往往披著各種“外衣”需要你識(shí)別并轉(zhuǎn)化為指派問(wèn)題模型。5.1 經(jīng)典應(yīng)用場(chǎng)景識(shí)別任務(wù)分配這是最直接的應(yīng)用。如論文評(píng)審分配每位評(píng)審審閱幾篇論文總匹配度最高、出租車派單車與乘客的距離最小、教室安排課程與教室的適配度。路徑規(guī)劃與排序某些旅行商問(wèn)題TSP的近似解法中會(huì)用到指派問(wèn)題來(lái)構(gòu)建匹配。例如將城市兩兩配對(duì)然后連接這些配對(duì)形成路徑。資源調(diào)度在固定時(shí)間段內(nèi)將機(jī)器分配給加工任務(wù)使得總加工時(shí)間最短或利潤(rùn)最大。圖像處理與數(shù)據(jù)關(guān)聯(lián)在多目標(biāo)跟蹤中將上一幀的檢測(cè)框與當(dāng)前幀的檢測(cè)框進(jìn)行關(guān)聯(lián)關(guān)聯(lián)成本可以是邊界框的重疊度IoU的負(fù)數(shù)。平衡實(shí)驗(yàn)設(shè)計(jì)將實(shí)驗(yàn)對(duì)象如患者分配到不同的實(shí)驗(yàn)組和控制組使得各組在某些特征上盡可能平衡這可以轉(zhuǎn)化為一個(gè)最小化組間差異的指派問(wèn)題。5.2 非標(biāo)準(zhǔn)情況的處理技巧1. 非平衡指派問(wèn)題工人數(shù) ≠ 工作數(shù)工人多工作少引入“虛擬工作”其成本設(shè)為0如果是最小化問(wèn)題。這意味著多余的工人沒(méi)有被指派任務(wù)成本為0。工人少工作多引入“虛擬工人”其完成所有工作的成本設(shè)為0。這意味著多余的工作沒(méi)有被完成成本為0。注意虛擬行/列的成本設(shè)置取決于問(wèn)題目標(biāo)。如果是最大化效益問(wèn)題虛擬行/列的效益通常設(shè)為0或一個(gè)非常大的負(fù)數(shù)在最大化問(wèn)題中表示不選擇。2. 最大化問(wèn)題標(biāo)準(zhǔn)匈牙利算法解決最小化問(wèn)題。對(duì)于最大化問(wèn)題如最大效益、最大匹配度常用方法有方法一將效益矩陣B轉(zhuǎn)化為成本矩陣C M - B其中M是矩陣B中元素的最大值或一個(gè)足夠大的數(shù)。然后對(duì)C求解最小化指派。因?yàn)镸in Σ(M - b_{ij})x_{ij} M*n - Max Σ b_{ij}x_{ij}所以解相同。方法二直接對(duì)效益矩陣B取負(fù)值C -B然后求解最小化指派。3. 禁止指派某些工人不能做某些工作。處理方法是將對(duì)應(yīng)成本設(shè)為一個(gè)極大的數(shù)INF。在最小化問(wèn)題中算法會(huì)主動(dòng)避免選擇成本為INF的配對(duì)。在代碼實(shí)現(xiàn)中可以用一個(gè)遠(yuǎn)大于其他正常成本的值如1e9來(lái)代替INF。4. 多對(duì)一或一對(duì)多指派標(biāo)準(zhǔn)指派是一對(duì)一。如果允許一個(gè)工人做多項(xiàng)工作或者一項(xiàng)工作需要多個(gè)工人問(wèn)題就變成了廣義分配問(wèn)題Generalized Assignment Problem, GAP這比標(biāo)準(zhǔn)指派問(wèn)題復(fù)雜得多通常需要用到更高級(jí)的整數(shù)規(guī)劃或啟發(fā)式算法如遺傳算法、模擬退火。此時(shí)匈牙利算法不再直接適用。5. 有額外約束的指派例如除了成本最小還要求某些工人必須被分配到一起或者某些工作必須在其他工作之后完成。這些約束破壞了二分圖匹配的結(jié)構(gòu)需要將其建模為更復(fù)雜的整數(shù)規(guī)劃問(wèn)題使用通用求解器如Gurobi, CPLEX或定制算法。5.3 建模實(shí)例數(shù)學(xué)建模競(jìng)賽題改編問(wèn)題某市有5個(gè)突發(fā)公共事件應(yīng)急點(diǎn)現(xiàn)有5支救援隊(duì)。已知各救援隊(duì)到達(dá)各應(yīng)急點(diǎn)的預(yù)計(jì)時(shí)間小時(shí)。由于專業(yè)設(shè)備限制第2支救援隊(duì)無(wú)法前往第3個(gè)應(yīng)急點(diǎn)。如何分配救援隊(duì)使得總響應(yīng)時(shí)間最短建模步驟定義決策變量x_{ij} 1表示派遣救援隊(duì)i到應(yīng)急點(diǎn)j否則為0。建立成本矩陣c_{ij}為行駛時(shí)間。對(duì)于禁止指派救援隊(duì)2 - 應(yīng)急點(diǎn)3令c_{23} INF一個(gè)大數(shù)如999。目標(biāo)函數(shù)Min Σ Σ c_{ij} x_{ij}。約束條件標(biāo)準(zhǔn)的指派問(wèn)題約束每行每列和為1。求解使用匈牙利算法或scipy.optimize.linear_sum_assignment求解。結(jié)果分析檢查最優(yōu)解中x_{23}是否為0驗(yàn)證禁止指派是否被遵守。計(jì)算總響應(yīng)時(shí)間。實(shí)操心得在論文寫作中將原始問(wèn)題抽象成矩陣形式是關(guān)鍵一步。建議在論文中清晰地畫出成本矩陣表格并對(duì)特殊值如INF加以說(shuō)明。這能讓評(píng)委一眼看出你正確理解了問(wèn)題并進(jìn)行了恰當(dāng)?shù)霓D(zhuǎn)化。6. 常見問(wèn)題、調(diào)試技巧與算法局限6.1 算法實(shí)現(xiàn)中的常見陷阱浮點(diǎn)誤差導(dǎo)致匹配失敗如前所述在判斷c_{ij} 0時(shí)使用絕對(duì)容差abs(cost[i][j]) 1e-10而不是cost[i][j] 0。非方陣處理不當(dāng)對(duì)于非平衡問(wèn)題務(wù)必先將其補(bǔ)全為方陣并正確設(shè)置虛擬行/列的成本。如果目標(biāo)是最大化補(bǔ)全時(shí)需要格外小心。無(wú)限循環(huán)在自編的完整迭代算法中如果矩陣調(diào)整步驟的邏輯有誤可能導(dǎo)致無(wú)法增加匹配數(shù)從而陷入無(wú)限循環(huán)。確保每次調(diào)整后至少有一個(gè)新的零元素在未被覆蓋的區(qū)域產(chǎn)生。誤用最大化算法直接將最大化問(wèn)題的矩陣輸入給最小化算法會(huì)得到錯(cuò)誤結(jié)果。務(wù)必先進(jìn)行轉(zhuǎn)化。6.2 調(diào)試與驗(yàn)證小規(guī)模驗(yàn)證用3x3或4x4的矩陣手動(dòng)計(jì)算與程序結(jié)果對(duì)比。這是最有效的調(diào)試方法。檢查解的可行性確保得到的指派方案滿足“每個(gè)代理恰好一個(gè)任務(wù)”的約束。計(jì)算row_ind和col_ind是否都是[0, 1, ..., n-1]的一個(gè)排列。與暴力枚舉對(duì)比對(duì)于n很小如n8的問(wèn)題可以編寫暴力枚舉所有排列的程序驗(yàn)證匈牙利算法給出的解是否確實(shí)是最優(yōu)的。使用庫(kù)函數(shù)交叉驗(yàn)證用scipy.optimize.linear_sum_assignment的結(jié)果來(lái)驗(yàn)證自己編寫的算法。6.3 匈牙利算法的局限與替代方案盡管匈牙利算法高效但它有其適用范圍僅適用于線性目標(biāo)函數(shù)總成本必須是各配對(duì)成本的和。如果是非線性如成本與配對(duì)順序有關(guān)則不適用。一對(duì)一嚴(yán)格約束這是核心假設(shè)。一對(duì)多、多對(duì)多需要其他模型。單目標(biāo)優(yōu)化只能處理最小化總成本或最大化總效益。多目標(biāo)指派問(wèn)題需要其他方法如目標(biāo)規(guī)劃、進(jìn)化算法。替代算法拍賣算法Auction Algorithm另一種求解指派問(wèn)題的經(jīng)典算法思想直觀模擬拍賣過(guò)程在某些情況下并行性好。最小成本最大流將指派問(wèn)題建模為網(wǎng)絡(luò)流問(wèn)題。源點(diǎn)連接所有工人容量1成本0工人連接所有工作容量1成本為c_ij工作連接匯點(diǎn)容量1成本0。求解從源到匯的最小成本最大流。這是一個(gè)更通用的框架可以處理更多變體。整數(shù)規(guī)劃求解器對(duì)于復(fù)雜約束的指派問(wèn)題直接使用Gurobi、CPLEX等求解器建模求解是最穩(wěn)妥的方式。雖然可能不如專用算法快但能保證在復(fù)雜約束下找到最優(yōu)解如果問(wèn)題可解。在我多年的建模和編程經(jīng)驗(yàn)中處理指派問(wèn)題的首選路徑是首先判斷是否是標(biāo)準(zhǔn)的一對(duì)一、線性成本問(wèn)題。如果是毫不猶豫地使用scipy.optimize.linear_sum_assignment。如果問(wèn)題帶有特殊約束如資源容量、先后順序則將其建立為整數(shù)規(guī)劃模型調(diào)用專業(yè)求解器。匈牙利算法的價(jià)值在于其優(yōu)美的理論和作為構(gòu)建更復(fù)雜算法基礎(chǔ)組件的作用但在實(shí)際應(yīng)用中我們更應(yīng)注重正確、高效地解決問(wèn)題而非重復(fù)造輪子。理解匈牙利算法就像是掌握了一把打開組合優(yōu)化大門的鑰匙。它讓你看到對(duì)于具有特殊結(jié)構(gòu)的整數(shù)規(guī)劃問(wèn)題存在比蠻力搜索和通用求解器更巧妙的道路。這種“發(fā)現(xiàn)結(jié)構(gòu)、利用結(jié)構(gòu)”的思維才是數(shù)學(xué)建模中最寶貴的財(cái)富。當(dāng)你下次遇到分配、匹配、調(diào)度類問(wèn)題時(shí)不妨先想一想這能不能抽象成一個(gè)二分圖能不能用匈牙利算法的思想來(lái)近似或求解這種思考習(xí)慣往往能讓你在競(jìng)賽中脫穎而出。