運籌優化（十五）--應用模型之分配問題

Eason.wxd發表於2019-01-21
優化模型

一般分配模型

一、問題描述

問題描述：N個人分配N項任務，一個人只能分配一項任務，一項任務只能分配給一個人，將一項任務分配給一個人是需要支付報酬，如何分配任務，保證支付的報酬總數最小。

問題數學描述：

　　

注意，這個規劃問題是整數線性規劃（ILP）問題，也就是說，兩個約束方程，保證每個任務被分配一次。決策變數僅允許取離散值0/1 

二、例項分析---窮舉法

在講將匈牙利演算法解決任務問題之前，先分析幾個具體例項。

以3個工作人員和3項任務為例項，下圖為薪酬圖表和根據薪酬圖表所得的cost矩陣。

　　

利用最簡單的方法(窮舉法)進行求解，計算出所有分配情況的總薪酬開銷，然後求最小值。

total_cost1 = 250 + 600 + 250 = 1100;  x00 = 1,x11 = 1,x22 = 1;

total_cost2 = 250 + 350 + 400 = 1000;  x00 = 1,x12 = 1,x21 = 1;

total_cost3 = 400 + 400 + 250 = 1050;  x01 = 1,x10 = 1,x22 = 1;

total_cost4 = 400 + 350 + 200 = 950;   x01 = 1,x12 = 1,x20 = 1;  //最優分配

total_cost5 = 350 + 400 + 400 = 1150; x02 = 1,x10 = 1,x21 = 1;

total_cost6 = 350 + 600 + 250 = 1150; x02 = 1,x11 = 1,x22 = 1;

對於任務數和人員數較少時，可利用窮舉法計算結果。

若將N任務分配給N個人員，其包含的所有分配情況數目為N!，N增大時，窮舉法將難以完成任務。

三、匈牙利演算法

下面簡要介紹匈牙利演算法。

其基本的理論基礎是針對cost矩陣，將cost矩陣的一行或一列資料加上或減去一個數，其最優任務分配求解問題不變。

　　

演算法的基本步驟如下：

　　

四、例項分析---匈牙利演算法

下面結合具體例項，分析匈牙利演算法如何解決任務分配問題。

以N = 4為例項，下圖為cost列表和cost矩陣。

　　

Step1.從第1行減去75，第2行減去35，第3行減去90，第4行減去45。

　　

Step2.從第1列減去0，第2列減去0，第3列減去0，第4列減去5。

　　

Step3.利用最少的水平線或垂直線覆蓋所有的0。

　　

Step4.由於水平線和垂直線的總數是3，少於4，進入Step5。

Step5.沒有被覆蓋的最小值是5，沒有被覆蓋的每行減去最小值5，被覆蓋的每列加上最小值5，然後跳轉到步驟3.

　　

Step3.利用最少的水平線或垂直線覆蓋所有的0。

　　

Step4.由於水平線和垂直線的總數是3，少於4，進入Step5。

Step5.沒有被覆蓋的最小值是20，沒有被覆蓋的每行減去最小值20，被覆蓋的每列加上最小值20，然後跳轉到步驟3.

　　

Step3.利用最少的水平線或垂直線覆蓋所有的0。

　　

Step4.由於水平線和垂直線的總數是4，演算法結束，分配結果如下圖所示。

　　

其中，黃色框表示分配結果，左邊矩陣的最優分配等價於左邊矩陣的最優分配。

注意：若分配問題中需要配對兩個集合大小不同，則小的一個可以利用偽成員進行擴張，這些偽成員可以分配給另外一個集合裡的所有成員，且相應的費用為零。

二次分配模型

我們在上面一般分配問題的約束下，最小化或者最大化如下形式的二次目標函式：

\sum_{i} \sum_{j} \sum_{k > i}\sum_{i \neq j} c_{i,j,k,l} x_{i,j}x_{k,l}

其中c_{i,j,k,l}是將i分配給j且k分配給l的收益或者成本

注意到，目標函式的每一項包含兩個分配決策：c_{i,j,k,l} x_{i,j}x_{k,l}，也就是隻有x_{i,j} = 1\, and \, x_{k,l} = 1才會有c_{i,j,k,l}的成本。

這種模型的應用，主要是在購物商場的店鋪裝修佈局。

注意：計算二次分配模型的全域性最優點是很困難的，這是因為目標函式變成了非線性，此時不再是整數規劃問題了，解這類問題通常採用啟發式演算法。

廣義分配模型

我們之前的分配問題約束是每個i必須被分配到一個j中，反之亦然。現在假設每一個i必須被分配給一些事物j，每一個事物j可以與多個i配對。具體來說，定義如下變數：

bj 表示j的容量

si,j表示如果i被分配j，那麼j需要消耗的容量、空間、或者類似的量

ci,j表示i分配給j的成本或者收益

此時，用來找到在不超過容量限制的情況下最佳分配所有事物i的方案的模型，就是廣義分配模型，形式如下：

min \, or \, max \sum_{i}\sum_{j} c_{i,j} x_{i,j}

s.t. \sum_{j} x_{i,j} = 1 \, \forall i

\sum_{j} s_{i,j}x_{i,j} \leq b_{j} \, \forall j

x_{i,j} =0 \, or \, 1 \, \forall i,j

廣義分配模型，還是屬於整數規劃問題，在規模適中時可以用整數規劃求解方式求解。

匹配模型

最後一種分配模型是消除集合之間的區別，即尋找同一類事物之間的最佳配對，此時，我們把一般分配模型中的j換成i',同時，約束條件變為：

\sum_{i' < i} x_{i',i} + \sum_{i' > i} x_{i,i'} = 1 \, \forall i

x_{i,i'} =0 \, or \, 1 \, \forall i,i' > i

一般來說，令i‘ > i的目的在於消去重複計數。

其應用如揚聲器匹配模型，同一批揚聲器由於生產工藝的差別，搭配成一個立體聲系統時，兩個揚聲器之間存在干擾失真，建模的意義在於，如何搭配兩個揚聲器，使得整體的干擾失真最小。

下面補充一點二分圖匹配問題，仔細體會，跟上面分配問題還是不太一樣，上面定義為整數線性規劃問題，下面採用圖的方式，差別在於，上面分配問題已經考慮了分配成本或者收益，下面圖分兩種，一種帶權（就是上面分配問題），一種不帶權。

以下參考：https://www.renfei.org/blog/bipartite-matching.html

無權二分圖

二分圖：簡單來說，如果圖中點可以被分為兩組，並且使得所有邊都跨越組的邊界，則這就是一個二分圖。準確地說：把一個圖的頂點劃分為兩個不相交集U和V，使得每一條邊都分別連線U,V的頂點。如果存在這樣的劃分，則此圖為一個二分圖。二分圖的一個等價定義是：不含有「含奇數條邊的環」的圖。圖 1 是一個二分圖。為了清晰，我們以後都把它畫成圖 2 的形式。

匹配：在圖論中，一個「匹配」（matching）是一個邊的集合，其中任意兩條邊都沒有公共頂點。例如，圖 3、圖 4 中紅色的邊就是圖 2 的匹配。

我們定義匹配點、匹配邊、未匹配點、非匹配邊，它們的含義非常顯然。例如圖 3 中 1、4、5、7 為匹配點，其他頂點為未匹配點；1-5、4-7為匹配邊，其他邊為非匹配邊。

最大匹配：一個圖所有匹配中，所含匹配邊數最多的匹配，稱為這個圖的最大匹配。圖 4 是一個最大匹配，它包含 4 條匹配邊。

完美匹配：如果一個圖的某個匹配中，所有的頂點都是匹配點，那麼它就是一個完備匹配。如果兩個集合的頂點數相同，且是完備匹配則成為一個完美匹配，圖 4 是一個完美匹配。顯然，完美匹配一定是最大匹配（完美匹配的任何一個點都已經匹配，新增一條新的匹配邊一定會與已有的匹配邊衝突）。但並非每個圖都存在完美匹配。

舉例來說：如下圖所示，如果在某一對男孩和女孩之間存在相連的邊，就意味著他們彼此喜歡。是否可能讓所有男孩和女孩兩兩配對，使得每對兒都互相喜歡呢？圖論中，這就是完美匹配問題。如果換一個說法：最多有多少互相喜歡的男孩/女孩可以配對兒？這就是最大匹配問題。

基本概念講完了。求解最大匹配問題的一個演算法是匈牙利演算法，下面講的概念都為這個演算法服務。

交替路：從一個未匹配點出發，依次經過非匹配邊、匹配邊、非匹配邊…形成的路徑叫交替路。

增廣路：從一個未匹配點出發，走交替路，如果途徑另一個未匹配點（出發的點不算），則這條交替路稱為增廣路（agumenting path）。例如，圖 5 中的一條增廣路如圖 6 所示（圖中的匹配點均用紅色標出）：

增廣路有一個重要特點：非匹配邊比匹配邊多一條。因此，研究增廣路的意義是改進匹配。只要把增廣路中的匹配邊和非匹配邊的身份交換即可。由於中間的匹配節點不存在其他相連的匹配邊，所以這樣做不會破壞匹配的性質。交換後，圖中的匹配邊數目比原來多了 1 條。

我們可以通過不停地找增廣路來增加匹配中的匹配邊和匹配點。找不到增廣路時，達到最大匹配（這是增廣路定理）。匈牙利演算法正是這麼做的。在給出匈牙利演算法 DFS 和 BFS 版本的程式碼之前，先講一下匈牙利樹。

匈牙利樹一般由 BFS 構造（類似於 BFS 樹）。從一個未匹配點出發執行 BFS（唯一的限制是，必須走交替路），直到不能再擴充套件為止。例如，由圖 7，可以得到如圖 8 的一棵 BFS 樹：

       

這棵樹存在一個葉子節點為非匹配點（7 號），但是匈牙利樹要求所有葉子節點均為匹配點，因此這不是一棵匈牙利樹。如果原圖中根本不含 7 號節點，那麼從 2 號節點出發就會得到一棵匈牙利樹。這種情況如圖 9 所示（順便說一句，圖 8 中根節點 2 到非匹配葉子節點 7 顯然是一條增廣路，沿這條增廣路擴充後將得到一個完美匹配）。

下面給出匈牙利演算法的 DFS 和 BFS 版本的程式碼：

// 頂點、邊的編號均從 0 開始
// 鄰接表儲存

struct Edge
{
    int from;
    int to;
    int weight;

    Edge(int f, int t, int w):from(f), to(t), weight(w) {}
};

vector<int> G[__maxNodes]; /* G[i] 儲存頂點 i 出發的邊的編號 */
vector<Edge> edges;
typedef vector<int>::iterator iterator_t;
int num_nodes;
int num_left;
int num_right;
int num_edges;
int matching[__maxNodes]; /* 儲存求解結果 */
int check[__maxNodes];

bool dfs(int u)
{
    for (iterator_t i = G[u].begin(); i != G[u].end(); ++i) { // 對 u 的每個鄰接點
        int v = edges[*i].to;
        if (!check[v]) {     // 要求不在交替路中
            check[v] = true; // 放入交替路
            if (matching[v] == -1 || dfs(matching[v])) {
                // 如果是未蓋點，說明交替路為增廣路，則交換路徑，並返回成功
                matching[v] = u;
                matching[u] = v;
                return true;
            }
        }
    }
    return false; // 不存在增廣路，返回失敗
}

int hungarian()
{
    int ans = 0;
    memset(matching, -1, sizeof(matching));
    for (int u=0; u < num_left; ++u) {
        if (matching[u] == -1) {
            memset(check, 0, sizeof(check));
            if (dfs(u))
                ++ans;
        }
    }
    return ans;
}
queue<int> Q;
int prev[__maxNodes];
int Hungarian()
{
    int ans = 0;
    memset(matching, -1, sizeof(matching));
    memset(check, -1, sizeof(check));
    for (int i=0; i<num_left; ++i) {
        if (matching[i] == -1) {
            while (!Q.empty()) Q.pop();
            Q.push(i);
            prev[i] = -1; // 設 i 為路徑起點
            bool flag = false; // 尚未找到增廣路
            while (!Q.empty() && !flag) {
                int u = Q.front();
                for (iterator_t ix = G[u].begin(); ix != G[u].end() && !flag; ++ix) {
                    int v = edges[*ix].to;
                    if (check[v] != i) {
                        check[v] = i;
                        Q.push(matching[v]);
                        if (matching[v] >= 0) { // 此點為匹配點
                            prev[matching[v]] = u;
                        } else { // 找到未匹配點，交替路變為增廣路
                            flag = true;
                            int d=u, e=v;
                            while (d != -1) {
                                int t = matching[d];
                                matching[d] = e;
                                matching[e] = d;
                                d = prev[d];
                                e = t;
                            }
                        }
                    }
                }
                Q.pop();
            }
            if (matching[i] != -1) ++ans;
        }
    }
    return ans;
}

匈牙利演算法的要點如下

  1. 從左邊第 1 個頂點開始，挑選未匹配點進行搜尋，尋找增廣路。

    1. 如果經過一個未匹配點，說明尋找成功。更新路徑資訊，匹配邊數 +1，停止搜尋。
    2. 如果一直沒有找到增廣路，則不再從這個點開始搜尋。事實上，此時搜尋後會形成一棵匈牙利樹。我們可以永久性地把它從圖中刪去，而不影響結果。

  2. 由於找到增廣路之後需要沿著路徑更新匹配，所以我們需要一個結構來記錄路徑上的點。DFS 版本通過函式呼叫隱式地使用一個棧，而 BFS 版本使用 prev 陣列。

效能比較

兩個版本的時間複雜度均為O(V⋅E)O(V⋅E)。DFS 的優點是思路清晰、程式碼量少，但是效能不如 BFS。測試了兩種演算法的效能。對於稀疏圖，BFS 版本明顯快於 DFS 版本；而對於稠密圖兩者則不相上下。在完全隨機資料 9000 個頂點 4,0000 條邊時前者領先後者大約 97.6%，9000 個頂點 100,0000 條邊時前者領先後者 8.6%, 而達到 500,0000 條邊時 BFS 僅領先 0.85%。

補充定義和定理：

最大匹配數：最大匹配的匹配邊的數目

最小點覆蓋數：選取最少的點，使任意一條邊至少有一個端點被選擇

最大獨立數：選取最多的點，使任意所選兩點均不相連

最小路徑覆蓋數：對於一個 DAG（有向無環圖），選取最少條路徑，使得每個頂點屬於且僅屬於一條路徑。路徑長可以為 0（即單個點）。

定理1：最大匹配數 = 最小點覆蓋數（這是 Konig 定理）

定理2：最大匹配數 = 最大獨立數

定理3：最小路徑覆蓋數 = 頂點數 - 最大匹配數

帶權二分圖

在一個二分圖內，左頂點為X，右頂點為Y，現對於每組左右連線XiYj有權wij，求一種匹配使得所有wij的和最大。也就是最大權匹配一定是完備匹配。如果兩邊的點數相等則是完美匹配。如果點數不相等，其實可以虛擬一些點，使得點數相等，也成為了完美匹配。最大權匹配還可以用最大流去解決。從這裡可以看出，本文開頭提到的一般分配問題，就是帶權二分圖的完美匹配問題。

Kuhn－Munkras演算法流程：

　　(1)初始化可行頂標的值

　　(2)用匈牙利演算法尋找完備匹配

　　(3)若未找到完備匹配則修改可行頂標的值

　　(4)重複(2)(3)直到找到相等子圖的完備匹配為止　

詳細的KM演算法這裡不提了，本文的主要目的在於區分清楚這幾個模型概念，且找到每個問題的最佳求解方式。

KM演算法是求最大權完備匹配，如果要求最小權完備匹配怎麼辦？方法很簡單，只需將所有的邊權值取其相反數，求最大權完備匹配，匹配的值再取相反數即可。

KM演算法的執行要求是必須存在一個完備匹配，如果求一個最大權匹配(不一定完備)該如何辦？依然很簡單，把不存在的邊權值賦為0。

KM演算法求得的最大權匹配是邊權值和最大，如果我想要邊權之積最大，又怎樣轉化？還是不難辦到，每條邊權取自然對數，然後求最大和權匹配，求得的結果a再算出e^a就是最大積匹配。至於精度問題則沒有更好的辦法了。

求最小(大)權匹配的費用流建模方法

求最小(大)權匹配，可以用最小(大)費用最大流的方法。和二分圖最大匹配的構圖方法類似，新增附加源S和附加匯T，從S向二分圖X集合中每個頂點連線一條權值為0，容量為1的有向邊，從Y集合中每個頂點向T也連線一條權值為0，容量為1的有向邊。然後把原有的邊變成容量為1，權值不變的有向邊。求從S到T的最小(大)費用最大流，就能求得最小(大)權匹配。

上述建模求最大權匹配的方法求得的一定是最佳匹配(如果存在完備匹配)，因為S到X集合每條邊全部滿流。如下圖所示，最小費用最大流為2。

要求最大權匹配(不一定完備匹配)。如下圖，只需再引入一個頂點A，從X集合的每個頂點向A連線一條容量為1，權值為0的邊，然後再由A向T連線一條權值為0，容量不小於|X|的邊，求最大費用最大流，這時是100。

最小權匹配也類似，不過新加的邊權要為一個極大值，大於所有已有邊權值。

KM演算法與費用流的比較

從理論上分析，KM演算法的時間複雜度比費用流要好，但是實際上和較好的費用流演算法比起來執行效率是差不多的，KM演算法優勢僅僅在於程式設計容易。對於十分稀疏的圖，許多優秀的費用流演算法效率是很高的。這並不說明KM演算法不如費用流，畢竟在資訊學競賽中，程式設計的複雜度也是一個相當重要的需要考慮的因素。

相關文章