分配問題與Hungarian算法

分配問題

指派問題

匈牙利算法

匈牙利方法是一種能夠在多項式時間內解決分配問題(assignment problem)的組合優化算法。它由Harold Kuhn 與1955年發展並提出，由於該算法很大程度上依賴於先前兩位匈牙利數學家:Denes Konig 和 Jeno Egervary,所以被命名為“匈牙利方法”。
1957年James Munkres重新審視了這個方法，證明發現該方法是嚴格polynomial的，所以之后該方法也被稱為Kuhn-Munkres 算法或者Munkres分配算法。原始的匈牙利算法的時間復雜度是,然而之后Edmonds和Karp，以及Tomizawa獨立發現經過一定的修改，該算法能改達到的時間復雜度。 Ford和Fulkerson將該方法擴展到一般運輸問題的求解上。 2006年，研究發現Carl Custav Jacobi在19實際就解決了assignment問題，並且在其逝世后的1890年求解過程被以拉丁語形式發表。。。

指派問題

匈牙利法解決的指派問題應該具有兩個約束條件

workes 和tasks的數目應該相同，即o2o問題。
求解的是最小化問題，如工作時間的最小化、費用的最小化等等

指派問題示例：
有三個workers: Jim, Steve和Alan，現在有3個工作：clean the bathroom, sweep the floors和wash the windows需要交給這三個人，每個人只能完成一個任務，對應的cost matrix如下

---	Clean bathroom	Sweep floors	Wash windows
Jim	$2	$3	$3
Steve	$3	$2	$3
Alan	$3	$3	$2

那么如何分配任務是開銷最小就是一個指派問題

匈牙利法步驟

問題： 假定某單位有甲、乙、丙、丁、戊五個員工，現需要完成A、B、C、D、E五項任務，每個員工完成某項任務的時間如下圖所示，應該如何分配任務，才能保證完成任務所需要的時間開銷最小？

1476015762594.jpg

解：

寫出系數矩陣
更新系數矩陣，使系數矩陣的每一行每一列都減去該行該列的最小值，保證每一行每一列都有0元素出現,參見定理2.
選擇只有一個0元素的行或列將該0元素標注為獨立0元素，並將該0元素所在的列或行中0元素划掉，直至找不到滿足條件的行或列，需要注意的是在循環時，划掉的0元素不再視為0元素。比如我們找下面這個系數矩陣的標注0元素和划掉的0元素

那么我們用1表示0元素，0表示非0元素，用2標注獨立0元素，-2表示划掉0元素，依次得到的中間結果為

可以發現我們在標注第一行第2個0元素時，並沒有把已經划掉的第一個0當作0元素看待。

划蓋0線。
a. 首先找到不含有獨立0元素的行，將行標注 (4)
b. 找到標注行中划掉的0元素所在的列，將列標注 (2,3)
c. 將標注列中獨立0元素所在的行標注 (4, 1, 2)
d. 重復b,c步驟知道所有的標注行中不再存在划掉的0元素
(行：4 ,1, 2, 3 ; 列: 2, 3, 1)
e. 在所有標注的列上做蓋0線，在所有未被標注的行上做蓋0線，我們在矩陣中用3表示被蓋0線覆蓋，則

可以發現所有的0元素都被覆蓋掉，所以稱為蓋0線。
根據定理1，如果此時蓋0線的個數等於矩陣的維數，相當於找到n個獨立0元素，則跳到步驟7，否則步驟5更新系數矩陣。
更新系數矩陣。
a. 找到未被蓋0線覆蓋的元素中的最小值
b. 所有未被蓋0線覆蓋的元素減去最小值
c. 所有蓋0線交叉處的元素值加上最小值
重復步驟4，5
計算最優解。
類似步驟3中找獨立0元素的方法在C中找到n個不同行不同列的0元素所對應的位置。

定理1. 系數矩陣C中獨立零元素的最多個數等於能覆蓋所有零元素的最少線數。 —— D. Konig
定理2. 若將分配問題的系數矩陣每一行及每一列分別減去各行及各列的最小元素，則新的分配問題與原分配問題有相同的最優解，只是最優值差一個常數。

當然，注意找蓋〇線的方法並不是唯一的，比如下述方法

同上一方法
同上一方法
在C中尋找未被蓋〇線覆蓋的存在0元素且數組最少的行(列)，標注一個0元素作為獨立0元素，該0元素所在的列(行)做蓋0線。重復此步驟，直至找不到符合條件的行或列，已經被找過的行(列)就不再找了！
若是蓋0線數組等於矩陣維度，則跳到7
同上一方法
重復3，4，5
同上一方法

當然還有別的方法，比如從包含最多0元素的行或列開始做蓋0線直到將所有的0元素覆蓋掉等

這里代碼實現了上面兩個方法，注釋掉的是第二個方法的核心

 
            function [cost,CMatrix]=Assignment(C,ismin)  
              
% Assignment problem solved by hungarian method.  
              
%  
              
% input:  
              
% C - 系數矩陣，可以適應workers和tasks數目不同的情形  
              
% ismin - 1表示最小化問題，0表示最大化問題  
              
% ouput:  
              
% cost - 最終花費代價  
              
% CMatrix - 對應的匹配矩陣，元素1所在位置c_{ij}表示j task分配給 i worker。  
              
%  
              
[m,n]=size(C);  
              
if ismin==0  
              
 C=max(C(:))-C;  
              
end  
              
%workes 和tasks數目不相同  
              
if m<n  
              
 C=[C;zeros(n-m,n)];  
              
elseif m>n  
              
 C=[C zeros(m,m-n)];  
              
end  
              
copyC=C;  
              
d=max(m,n);% 最終系數矩陣的維度  
              
C=C-repmat(min(C,[],2),1,d);  
              
C=C-repmat(min(C,[],1),d,1);  
              
%% 方法一  
              
while 1  
              
 A=int8((C==0));  
              
 nIZeros=0; % 獨立0元素的個數  
              
 while 1  
              
 r=sum(A==1,2); % 每一行0元素的個數  
              
 [~,idr]=find(r'==1);%找到只有一個0元素的行  
              
 if ~isempty(idr) % 如果找到這樣的行  
              
 tr=A(idr(1),:);  
              
 [~,idc]=find(tr==1);%找到0元素所在列  
              
 A(idr(1),idc)=2;%標注獨立元素  
              
 tc=A(:,idc);  
              
 tc(idr(1))=2;  
              
 [~,idr]=find(tc'==1);%找到獨立0元素所在列的其他0元素  
              
 A(idr,idc)=-2;%划掉獨立0元素所在列的其余0元素  
              
 nIZeros=nIZeros+1;  
              
 else  
              
 c=sum(A==1,1); % 每一列0元素的個數  
              
 [~,idc]=find(c==1);%找到只含有一個0元素的列  
              
 if ~isempty(idc)% 找到這樣的列  
              
 tc=A(:,idc(1));  
              
 [~,idr]=find(tc'==1);%0元素所在的行  
              
 A(idr,idc(1))=2;%標注獨立0元素  
              
 tr=A(idr,:);  
              
 tr(idc(1))=2;  
              
 [~,idc]=find(tr==1);%獨立0元素所在行的其他0元素  
              
 A(idr,idc)=-2;%划掉獨立0元素所在行的其余0元素  
              
 nIZeros=nIZeros+1;  
              
 else  
              
 break;  
              
 end  
              
 end  
              
 end  
              
 if nIZeros==d  
              
 %計算最優解  
              
 CMatrix=(A==2);  
              
 if ismin==1  
              
 cost=sum(copyC(:).*CMatrix(:));  
              
 else  
              
 cost = sum((max(copyC(:))-copyC(:)).*CMatrix(:));  
              
 end  
              
 CMatrix=CMatrix(1:m,1:n);  
              
 break;%找到d個獨立0元素，則跳出循環  
              
 else% 獨立0元素個數不足，就要找蓋0線了  
              
 r=sum(A==2,2);  
              
 [~,idr]=find(r'==0);%不含有獨立0元素的行  
              
 idrr=idr;  
              
 idcc=[];  
              
 while 1  
              
 tr=A(idrr,:);  
              
 [~,idc]=find(tr==-2);%不含獨立0元素的行中划掉的0元素所在列  
              
 if isempty(idc),break;end  
              
 tc=A(:,unique(idc));  
              
 [idrr,~]=find(tc==2);%這些列中標注的0元素所在行  
              
 idr=[idr,idrr'];  
              
 idcc=[idcc,idc];  
              
 end  
              
 idry=1:d;  
              
 idry(idr)=[];%蓋0線所在的行索引  
              
 TempC=C;%存儲非覆蓋元素  
              
 TempC(idry,:)=[];  
              
 TempC(:,idcc)=[];  
              
 minUnOverlap=min(TempC(:));  
              
 %更新系數矩陣  
              
 C=C-minUnOverlap;  
              
 C(idry,:)=C(idry,:)+minUnOverlap;  
              
 C(:,idcc)=C(:,idcc)+minUnOverlap;  
              
 end  
              
end  
              
%% 方法二  
              
% while 1  
              
% CMatrix=zeros(d);  
              
% nLines=0;  
              
% A=(C==0);  
              
% idx=[];  
              
% idy=[];  
              
% sr=[];  
              
% sc=[];  
              
% while 1  
              
% r=sum(A,2);  
              
% c=sum(A,1);  
              
% r(sr)=0;  
              
% c(sc)=0;  
              
% trc=[r(:);c(:)];  
              
% [trc,idtrc]=sort(trc,1,'ascend');  
              
% [~,idn0]=find(trc'>0);  
              
% if ~isempty(idn0)  
              
% id=idtrc(idn0(1));  
              
% if id>d  
              
% tc=A(:,id-d);  
              
% [~,idr]=find(tc'==1);  
              
% A(idr(1),:)=0;   
              
% nLines=nLines+1;  
              
% idy=[idy,idr(1)];  
              
% CMatrix(idr(1),id-d)=1;  
              
% sc=[sc,id-d];  
              
% else  
              
% tr=A(id,:);  
              
% [~,idc]=find(tr==1);  
              
% A(:,idc(1))=0;  
              
% nLines=nLines+1;  
              
% idx=[idx,idc(1)];  
              
% CMatrix(id,idc(1))=1;  
              
% sr=[sr,id];  
              
% end   
              
% else  
              
% break;  
              
% end  
              
% end  
              
% if nLines==d  
              
% if ismin  
              
% cost=sum(copyC(:).*CMatrix(:));  
              
% else  
              
% cost=sum((max(copyC(:))-copyC(:)).*CMatrix(:));  
              
% end  
              
% CMatrix=CMatrix(1:m,1:n);  
              
% break;  
              
% else  
              
% tempC=C;  
              
% tempC(idy,:)=[];  
              
% tempC(:,idx)=[];  
              
% minUnOverlap=min(tempC(:));  
              
% C=C-minUnOverlap;  
              
% C(idy,:)=C(idy,:)+minUnOverlap;  
              
% C(:,idx)=C(:,idx)+minUnOverlap;  
              
% end  
              
% end  
              
end