匈牙利算法 指派问题 Lingo&Matlab实现

Lingo实现（线性规划数学模型直接求解）：
5个人完成4项工作，时间成本最小。

dat1.txt数据内容：
3 8 2 10
8 7 2 9
6 4 2 7
8 4 2 3
9 10 6 9

model:
sets:
ren/1..5/;!5 persons;
wu/1..4/;!4 tasks;
link(ren,wu):c,x;
endsets

data:
c=@FILE('dat1.txt');
enddata

min=@sum(link:c*x);
@for(ren(i):@sum(wu(j):x(i,j))<=1);!如果人数和任务数相等，则为等号;
@for(wu(j):@sum(ren(i):x(i,j))=1);
@for(link:@bin(x));
end

Matlab实现（匈牙利算法，矩阵划线求解）：
成本矩阵：

% 测试用例
dat=[1 2 3 4 5;
    2 4 6 3 4;
    3 5 4 2 3;
    4 3 5 4 1;
    1,6,5,4,3];
[x,y]=Hungarian(dat)

% Hungarian.m
%指派问题的匈牙利算法，输入矩阵，a(ij)为i指派给j,第i人干第j个工作
function [Matching,Cost] = Hungarian(Perf)
% 
% 用于在给定MxN边缘的情况下找到最小边缘权重匹配的函数
% 使用匈牙利算法的权重矩阵权重。
% 
% Inf的边权重表示由其给出的顶点对
% 位置没有相邻的边缘。
% 
% MATCHING返回一个MxN矩阵，1表示匹配的位置
% COST返回最小匹配的成本
% 
% 
% 
% 初始化变量

 Matching = zeros(size(Perf));

% ？通过删除任何未连接的顶点来压缩性能矩阵
% ？会增加算法的速度

  % 找到每个已连接的列中的数字
    num_y = sum(~isinf(Perf),1);%?
  % 找到连接的每一行中的数字
    num_x = sum(~isinf(Perf),2);%?
    
  % 找到隔离的列（顶点）和行（顶点）
    x_con = find(num_x~=0);% find()返回矩阵或向量中非零元素的索引
    y_con = find(num_y~=0);
    
  % ?汇总凝聚性能矩阵
    P_size = max(length(x_con),length(y_con));
    P_cond = zeros(P_size);
    P_cond(1:length(x_con),1:length(y_con)) = Perf(x_con,y_con);
    if isempty(P_cond)
      Cost = 0;
      return
    end

    % 确保存在完美匹配
   %计算边缘矩阵的形式
      Edge = P_cond;
      Edge(P_cond~=Inf) = 0;
      % 找到Edge Matrix中的缺陷（CNUM）
      cnum = min_line_cover(Edge);
    
      % 投影额外的顶点和边缘，以便完美匹配
      % exists
      Pmax = max(max(P_cond(P_cond~=Inf)));%不等于无穷值的其他值的最大值
      P_size = length(P_cond)+cnum;
      P_cond = ones(P_size)*Pmax;
      P_cond(1:length(x_con),1:length(y_con)) = Perf(x_con,y_con);
   
%*************************************************
% 主程序：执行步骤的控制
%*************************************************
  exit_flag = 1;
  stepnum = 1;
  while exit_flag
    switch stepnum
      case 1 % 找到每行的最小数字并且每一行减去它
        [P_cond,stepnum] = step1(P_cond);
      case 2 % 划线
        [r_cov,c_cov,M,stepnum] = step2(P_cond);
      case 3 % 判断是不是划线覆盖所有值
        [c_cov,stepnum] = step3(M,P_size);
      case 4 % 调整划线，使覆盖所有0
        [M,r_cov,c_cov,Z_r,Z_c,stepnum] = step4(P_cond,r_cov,c_cov,M);
      case 5 % 未划线的0所在行没有星标0
        [M,r_cov,c_cov,stepnum] = step5(M,Z_r,Z_c,r_cov,c_cov);
      case 6 % 未划线的0所在行有星标0,变换矩阵以增加0
        [P_cond,stepnum] = step6(P_cond,r_cov,c_cov);
      case 7
        exit_flag = 0;
    end
  end

%删除所有虚拟卫星和目标并取消激活
%匹配原始性能矩阵的大小。
Matching(x_con,y_con) = M(1:length(x_con),1:length(y_con));
Cost = sum(sum(Perf(Matching==1)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   STEP 1: 找到每行的最小数字并且每一行减去它
% Find the smallest number of zeros in each row
%           and subtract that minimum from its row
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [P_cond,stepnum] = step1(P_cond)

  P_size = length(P_cond);% 返回行数和列数较大的那个值
  
  % Loop throught each row
  for ii = 1:P_size
    rmin = min(P_cond(ii,:));
    P_cond(ii,:) = P_cond(ii,:)-rmin;
  end

  stepnum = 2;
  
%**************************************************************************  
%   STEP 2: 划线：找到一个0，如果这个0所在行没有被星标的0，从这个0开始。重复每个0.
% Find a zero in P_cond. If there are no starred zeros in its
%           column or row start the zero. Repeat for each zero
%**************************************************************************

function [r_cov,c_cov,M,stepnum] = step2(P_cond)

% Define variables
  P_size = length(P_cond);
  r_cov = zeros(P_size,1);  % 显示是否覆盖行的向量
  c_cov = zeros(P_size,1);  % 显示是否覆盖列的向量
  M = zeros(P_size);        %一个掩码，显示一个位置是星号还是已准备好
  
  for ii = 1:P_size
    for jj = 1:P_size
      if P_cond(ii,jj) == 0 && r_cov(ii) == 0 && c_cov(jj) == 0
        M(ii,jj) = 1;
        r_cov(ii) = 1;
        c_cov(jj) = 1;
      end
    end
  end
  
% Re-initialize the cover vectors
  r_cov = zeros(P_size,1);  % A vector that shows if a row is covered
  c_cov = zeros(P_size,1);  % A vector that shows if a column is covered
  stepnum = 3;
  
%**************************************************************************
%   STEP 3: 从星标的0开始划线，如果所有列都被覆盖，那么这个匹配是最大匹配。
% Cover each column with a starred zero. If all the columns are
%           covered then the matching is maximum
%**************************************************************************

function [c_cov,stepnum] = step3(M,P_size)

  c_cov = sum(M,1);
  if sum(c_cov) == P_size
    stepnum = 7;
  else
    stepnum = 4;
  end
  
%**************************************************************************
%   STEP 4: 调整划线，使覆盖所有0：找到一个未被覆盖的0，置顶它。置顶0所在行没有
%           星标0，则进入第5步————。否则，有星标0，对这一行划线，并且取消包含这个星标0的列的线。
%           重复这一步骤直到没有剩余的未被覆盖的0。保留这个最小覆盖值并且进入步骤6.
%           Find a noncovered zero and prime it.  If there is no starred
%           zero in the row containing this primed zero, Go to Step 5.  
%           Otherwise, cover this row and uncover the column containing 
%           the starred zero. Continue in this manner until there are no 
%           uncovered zeros left. Save the smallest uncovered value and 
%           Go to Step 6.
%**************************************************************************
function [M,r_cov,c_cov,Z_r,Z_c,stepnum] = step4(P_cond,r_cov,c_cov,M)

P_size = length(P_cond);

zflag = 1;
while zflag  
    % Find the first uncovered zero
      row = 0; col = 0; exit_flag = 1;
      ii = 1; jj = 1;
      while exit_flag
          if P_cond(ii,jj) == 0 && r_cov(ii) == 0 && c_cov(jj) == 0
            row = ii;
            col = jj;
            exit_flag = 0;
          end      
          jj = jj + 1;      
          if jj > P_size; jj = 1; ii = ii+1; end      
          if ii > P_size; exit_flag = 0; end      
      end

    % If there are no uncovered zeros go to step 6
      if row == 0
        stepnum = 6;
        zflag = 0;
        Z_r = 0;
        Z_c = 0;
      else
        % Prime the uncovered zero
        M(row,col) = 2;
        % If there is a starred zero in that row
        % Cover the row and uncover the column containing the zero
          if sum(find(M(row,:)==1)) ~= 0
            r_cov(row) = 1;
            zcol = find(M(row,:)==1);
            c_cov(zcol) = 0;
          else
            stepnum = 5;
            zflag = 0;
            Z_r = row;
            Z_c = col;
          end            
      end
end
  
%**************************************************************************
% STEP 5: 调整划线
%         Z0表示第4步中发现的未被覆盖的置顶0；Z1表示Z0所在列的星标0（如果有）；
%         Z2表示Z1所在行的置顶0（总是会有一个）.
%         直到置顶0所在列没有星标0时终止。
%         取消星标这个系列的每个星标0，星标这个系列的每个置顶0
%         擦除所有置顶，取消所有划线。
%         返回步骤3.
%         Construct a series of alternating primed and starred zeros as
%         follows.  Let Z0 represent the uncovered primed zero found in Step 4.
%         Let Z1 denote the starred zero in the column of Z0 (if any). 
%         Let Z2 denote the primed zero in the row of Z1 (there will always
%         be one).  Continue until the series terminates at a primed zero
%         that has no starred zero in its column.  Unstar each starred 
%         zero of the series, star each primed zero of the series, erase 
%         all primes and uncover every line in the matrix.  Return to Step 3.
%**************************************************************************

function [M,r_cov,c_cov,stepnum] = step5(M,Z_r,Z_c,r_cov,c_cov)

  zflag = 1;
  ii = 1;
  while zflag 
    % Find the index number of the starred zero in the column
    rindex = find(M(:,Z_c(ii))==1); %置顶0所在列的星标0
    if rindex > 0
      % Save the starred zero
      ii = ii+1;
      % Save the row of the starred zero
      Z_r(ii,1) = rindex;
      % The column of the starred zero is the same as the column of the 
      % primed zero
      Z_c(ii,1) = Z_c(ii-1);
    else
      zflag = 0;
    end
    
    % Continue if there is a starred zero in the column of the primed zero
    if zflag == 1;
      % Find the column of the primed zero in the last starred zeros row
      cindex = find(M(Z_r(ii),:)==2);% 上一个星标0所在行有无置顶0
      ii = ii+1;
      Z_r(ii,1) = Z_r(ii-1);
      Z_c(ii,1) = cindex;    
    end    
  end
  
  % UNSTAR all the starred zeros in the path and STAR all primed zeros
  for ii = 1:length(Z_r)
    if M(Z_r(ii),Z_c(ii)) == 1
      M(Z_r(ii),Z_c(ii)) = 0; %星标0取消星标
    else
      M(Z_r(ii),Z_c(ii)) = 1; %置顶0全部星标
    end
  end
  
  % Clear the covers
  r_cov = r_cov.*0;
  c_cov = c_cov.*0;
  
  % Remove all the primes
  M(M==2) = 0;

stepnum = 3;

% *************************************************************************
% STEP 6: 已划线行每个元素加上未覆盖的最小值，然后从未划线的列中减去。不改变星标0、
%         置顶0和划线，返回步骤4.
%         Add the minimum uncovered value to every element of each covered
%         row, and subtract it from every element of each uncovered column.  
%         Return to Step 4 without altering any stars, primes, or covered lines.
%**************************************************************************

function [P_cond,stepnum] = step6(P_cond,r_cov,c_cov)
a = find(r_cov == 0);
b = find(c_cov == 0);
minval = min(min(P_cond(a,b)));

P_cond(find(r_cov == 1),:) = P_cond(find(r_cov == 1),:) + minval;
P_cond(:,find(c_cov == 0)) = P_cond(:,find(c_cov == 0)) - minval;

stepnum = 4;

function cnum = min_line_cover(Edge)

  % Step 2
    [r_cov,c_cov,M,stepnum] = step2(Edge);
  % Step 3
    [c_cov,stepnum] = step3(M,length(Edge));
  % Step 4
    [M,r_cov,c_cov,Z_r,Z_c,stepnum] = step4(Edge,r_cov,c_cov,M);
  % Calculate the deficiency
    cnum = length(Edge)-sum(r_cov)-sum(c_cov);

参考资料

1、匈牙利算法——MATLAB代码https://blog.csdn.net/qq_36980029/article/details/105866390
2、分派问题（匈牙利算法）与MATLAB实现https://blog.csdn.net/qq_36980029/article/details/96626332
3、https://blog.csdn.net/Richard_pl/article/details/79291463
4、https://blog.csdn.net/tommy0095/article/details/104466364/
5、https://blog.csdn.net/Wonz5130/article/details/80678410