匈牙利算法——二分图匹配
第一步，求出每个正负残差点间的位置距离，构造系数矩阵C。
第二步，针对系数矩阵建立等效矩阵B，使各行各列都出现0元素，方法如下：
（1）把系数矩阵的每行元素减去该行的最小元素。
（2）再把所得矩阵的每列元素减去该列的最小元素。
第三步，进一步处理矩阵——划线（即匈牙利算法中覆盖所有0的线数量等于矩阵的维度，说明处理好了）。
划线的处理思路是，（1）先看行，如果这个行中的0的数量大于1，则可以选择划线，line_count可以+1，并且删除掉这一行。（2）所有行中0数量大于1的处理完之后，删除掉这些行形成新的矩阵，用同样的思路处理列，在形成新的矩阵。（3）这样新的矩阵中就只剩下一种情况，每行中只有一个0（每列中只有一个0），用行算就行。
第四步，分配
分配的思路是，行或列只有1个0的必然首先分配，并且这个0所在的行和列要不能出现0，所以直接行和列变成inf，原本想删除掉，发现有问题。最后可能会出现[0 0 ；0 0 ]这种情况，所以判断了多个0的时候直接取第一个。缺点就是没有办法显示多组解。
主程序：
1 clc; 2 clear; 3 % 主程序：基于匈牙利算法 4 % 输入的矩阵是前面计算好的COST（calculating_price.m) 5 % 需要对矩阵进行处理，基站列需要扩充，基站最多要满足M个CR用户的请求;行也需要调整。总之就是要将矩阵处理成方阵 6 7 %% 矩阵处理 - 处理成符合意义的方阵 8 % load('COST.mat'); 9 COST = [1 4 3;2 inf 1;3 6 9]; 10 [x,y] = size(COST); 11 COST_P = COST; 12 if x ~= y 13 % 赋值BS列 14 BS_column = repmat(COST(:,y),1,x-1); 15 COST_P = [COST_P BS_column]; 16 % 矩阵补零变成方阵 17 zero = zeros(1,x+y-1); 18 Add_zero = repmat(zero,y-1,1); 19 COST_P = [COST_P;Add_zero]; 20 end% 处理之后的COST_P矩阵应该是一个方阵 21 COST_expand = COST_P; 22 [m,n] = size(COST_P); 23 24 %% 处理矩阵 - 匈牙利算法 25 COST_P = rc(COST_P); 26 %% 计算覆盖的行数和列数 27 [line_sum,row_index1,row_sum1,column_index,column_sum,row_index2,row_sum2,matrix_out] = line_count(COST_P); 28 %% 处理划线的数量小于矩阵的维数 29 while (line_sum < m) 30 rest_min = min(min(matrix_out)); % 找到没有被覆盖的最小值，将没有被覆盖的每行减去最小值，被覆盖的每列加上最小值。转到上一步。 31 row_num = 1:1:m; 32 row_index = sort([row_index1 row_index2]); 33 row_num(row_index) = [];% 找到没有被覆盖的行 34 COST_P(row_num,:) = COST_P(row_num,:) - rest_min;% 没有被覆盖的行减去最小值 35 COST_P(:,column_index) = COST_P(:,column_index) + rest_min;% 被覆盖的列加上最小值 36 [line_sum,row_index1,row_sum1,column_index,column_sum,row_index2,row_sum2,matrix_out] = line_count(COST_P); 37 end 38 %% 输出匹配结果 39 [asig] = fenpei(COST_P); 40 if x ~= y 41 del_index =x+1:1:x+y-1; 42 asig(del_index) = [] 43 end 44 %% 计算COST 45 cost_sum = 0; 46 n = length(asig); 47 for i = 1:n 48 cost_sum = cost_sum + COST_expand(i,asig(i)); 49 end 50 cost_sum
行处理和列处理
1 function [pro_rc] = rc(matrix) 2 % 匈牙利算法中处理行和处理列，rc表示row和column 3 % 算法的第一步行减去当前行的最小值以及列减去当前列的最小值 4 % 输入参数是填入待处理的矩阵；输出参数是处理之后的矩阵 5 m = length(matrix); 6 for i=1:m 7 rowi_min = min(matrix(i,:)); 8 matrix(i,:) = matrix(i,:)-rowi_min; 9 end 10 for i=1:m 11 columni_min = min(matrix(:,i)); 12 matrix(:,i) = matrix(:,i)-columni_min; 13 end 14 pro_rc = matrix; 15 end
进一步处理矩阵——划线
1 function [line_sum,row_index,row_sum,column_index,column_sum,row_index2,row_sum2,matrix_out] = line_count(matrix) 2 % 对处理过的矩阵进行划线统计 3 % 在匈牙利算法中对行和列进行划线的时候，是可以先从行开始计算，先检索有不少于两个0的行，有的话就划掉，这一个行就直接从 4 % 矩阵中删除掉形成新的矩阵，接着循环；列也是这样。 5 % 等到矩阵中都没有不少于两个0的行和列的时候，就直接以行为准，划零，统计行数； 6 % line_sum是划线总数,row_sum是统计行划了几条线，column_sum是统计列划了几条线，matrix_out是看最后画完线的矩阵 7 [m,n] = size(matrix); %统计矩阵的维度m*n 8 line_sum = 0;%存放划0的数量 9 10 %% 处理行 11 row_index = []; 12 j = 1; 13 for i=1:m 14 a = length(find(matrix(i,:)==0)); 15 if a>1 16 line_sum = line_sum + 1; 17 row_index(j) = i; 18 j = j + 1; 19 end 20 end 21 row_sum = length(row_index); 22 matrix(row_index,:) =[]; %将输入矩阵删除掉行中含不少于2个0的行 23 24 %% 处理列 25 column_index = []; 26 k = 1; 27 for i=1:n 28 b = length(find(matrix(:,i)==0)); 29 if b>1 30 line_sum = line_sum +1; 31 column_index(k) = i; 32 k = k + 1; 33 end 34 end 35 column_sum = length(column_index); 36 matrix(:,column_index) = []; 37 38 %% 处理只有一个0 39 [m,n] = size(matrix);%上一步处理完的矩阵有可能不是方阵 40 row_index2 = []; 41 j = 1; 42 for i=1:m 43 a = length(find(matrix(i,:)==0)); 44 if a == 1 45 line_sum = line_sum +1; 46 row_index2(j) = i; 47 j = j+1; 48 end 49 end 50 row_sum2 = length(row_index2); 51 matrix(row_index2,:) = []; 52 matrix_out = matrix; 53 end
分配
1 function [asig] = fenpei(matrix) 2 % 函数用于最后的分配 3 % 函数输入是最后的0划线数量等于矩阵的维度的最终被处理过的矩阵 4 % 函数的输出是最后的分配方案，[2 3 1 4]表示User1被分配到Task2，依次类推 5 [m,n] = size(matrix); 6 asig = [];% 存放决策的矩阵,[4 2 3 1]表示第一个用户被分配到任务4，以此类推 7 p = 1; 8 q = 1; 9 temp = length(find(matrix==0)); 10 while (temp~=0) 11 %% 处理1个0 12 for i = 1:m 13 zero = length(find(matrix(i,:)==0)); 14 if zero == 1 15 row_0_1(p) = i;% row_0_1表示只有1个0的行 16 column_row_0_1(p) = find(matrix(i,:)==0);% column_row_0_1表示只与一个0的行对应的列 17 asig(i) = find(matrix(i,:)==0);% 用户i被分配任务 18 p = p+1; 19 end 20 end 21 22 for j = 1:n 23 zero = length(find(matrix(:,j)==0)); 24 if zero == 1 25 column_0_1(q) = j; 26 row_column_0_1(q) = find(matrix(:,j)==0); 27 asig(row_column_0_1(q)) = j;% 任务可被分配 28 q = q+1; 29 end 30 end 31 row_delete = [row_0_1 row_column_0_1]; 32 column_delete = [column_row_0_1 column_0_1]; 33 matrix(row_delete,:) = matrix(row_delete,:)+inf; 34 matrix(:,column_delete) = matrix(:,column_delete)+inf; 35 36 %% 处理多余1个0 37 for i = 1:m 38 zero = find(matrix(i,:)==0); 39 if length(zero) ~= 0 40 asig(i) = zero(1); 41 matrix(i,:) = matrix(i,:) + inf; 42 matrix(:,zero(1)) = matrix(:,zero(1)) + inf; 43 end 44 end 45 46 temp = length(find(matrix==0)); 47 end 48 end
==========================下面是我写的代码===============================================
主程序
1 function res=hungary(N) 2 %输入的矩阵应N*N的 3 C = [12 7 9 7 9; 4 8 9 6 6 6; 5 7 17 12 14 9; 6 15 14 6 6 10; 7 4 10 7 10 9]; 8 9 N = C; %复制矩阵，便于后续处理 10 11 [a,b]=size(N); 12 %====================第1步,每一行减去当前行最小值 13 for ii = 1:a 14 N(ii,:)= N(ii,:)-min( N(ii,:)); 15 end 16 %====================第2步,每一列减去当前列最小值 17 for ii = 1:a 18 N(:,ii)= N(:,ii)-min( N(:,ii)); 19 end 20 21 line = 0; %划线数 22 %==================================== 第5步，更新矩阵 23 %对没有被线划到的数中，找到最小的数 24 while line ~=a 25 [line, N_min,del_hang, del_lie, r]=lineCount(N); 26 if line~=a 27 for i = 1:a 28 if del_hang(i)~=i 29 N(i,:)= N(i,:)-N_min; 30 end 31 if del_lie(i)==i 32 N(:,i)=N(:,i) + N_min; 33 end 34 end 35 else 36 res = N; 37 end 38 end 39 40
划线
1 function [line,N_min,del_hang,del_lie, r] = lineCount(C) 2 M = C; %复制矩阵,便于做后续处理 3 A = C; 4 [a,b]=size(M); %统计矩阵的维度 5 Q = zeros(size(M)); % 0:未被画线 1：画了一次 2：画了2次（交点） 6 r = zeros(size(M)); % 0:非0元素 1:非独立0元素 2:独立0元素 7 8 depzero=0; % 独立0元素个数 9 10 x = zeros(1,a); % 画线时行是否被打勾，1是 0不是 11 y = zeros(1,a); % 画线时列是否被打勾，1是 0不是 12 13 %获取每行和每列0元素的个数， 14 rz = zeros(1,a); 15 cz = zeros(1,a); 16 17 %行列0元素个数 18 del_hang=zeros(a,1); 19 del_lie=zeros(a,1); 20 21 for i = 1:a 22 for j = 1:a 23 %找出非独立0元素位置 24 if(M(i,j)==0) 25 r(i,j)=1; 26 end 27 rz(i)=length(find(M(i,:)==0)); 28 %r0[i] = inf表示该行无0元素 29 if rz(i) == 0 30 rz(i) = inf; 31 end 32 cz(i)=length(find(M(:,i)==0)); 33 if cz(i) == 0 34 cz(i) = inf; 35 end 36 end 37 end 38 39 %===================第3步，试指派 40 %获取独立0元素个数 41 while find(M==0) 42 if min(rz)<=min(cz) %行最少0元素 43 %找含0最少那一行的行坐标 44 [~,rmin_index] = find(rz==min(min(rz))); 45 %列坐标 46 cmin_index=find(M(rmin_index(1),:)==0); 47 r(rmin_index(1),cmin_index(1))=2; %标记独立0元素位置 48 depzero = depzero+1; % 独立0元素个数 49 % 把这个0所在的行和列变成inf 50 M(rmin_index(1),:)=1./zeros(1,a); 51 M(:,cmin_index(1))=1./zeros(a,1); 52 %置0，下一循环要用 53 rmin_index=0; 54 cmin_index=0; 55 56 % 更新每行和每列0元素的个数 57 for i = 1:a 58 %for j = 1:a 59 rz(i)=length(find(M(i,:)==0)); 60 %r0[i] = inf表示该行无0元素 61 if rz(i) == 0 62 rz(i) = inf; 63 end 64 cz(i)=length(find(M(:,i)==0)); 65 if cz(i) == 0 66 cz(i) = inf; 67 end 68 end 69 70 else %列最少0元素 71 %找含0最少那一行的列坐标 72 [~,cmin_index] = find(cz==min(min(cz))); 73 %行坐标 74 rmin_index = find(M(:,cmin_index(1))==0); 75 r(rmin_index(1),cmin_index(1))=2; %标记独立0元素位置 76 depzero = depzero+1; % 独立0元素个数 77 % 把这个0所在的行和列变成inf 78 M(rmin_index(1),:)=1./zeros(1,a); 79 M(:,cmin_index(1))=1./zeros(a,1); 80 %置0，下一循环要用 81 rmin_index=0; 82 cmin_index=0; 83 84 % 更新每行和每列0元素的个数 85 for i = 1:a 86 rz(i)=length(find(M(i,:)==0)); 87 %r0[i] = inf表示该行无0元素 88 if rz(i) == 0 89 rz(i) = inf; 90 end 91 cz(i)=length(find(M(:,i)==0)); 92 if cz(i) == 0 93 cz(i) = inf; 94 end 95 end 96 end 97 98 end 99 100 %======================第4步，画0盖线 101 %打勾初始化，画线时是否被打勾，1是0不是 102 for i = 1:a 103 x(i) =1; 104 y(i) =0; 105 end 106 % 1.row 对所有不含独立0元素的行打勾 107 for i = 1:a 108 for j = 1:a 109 if(r(i,j)==2) 110 x(i)=0; 111 end 112 end 113 end 114 aq = 1; 115 while aq == 1 116 aq = 0; 117 % 2.col 对打勾的行中含0元素的未打勾的列打勾 118 for i = 1:a 119 if (x(i)==1) 120 for j = 1:a 121 if(C(i,j)==0 && y(j)==0) 122 y(j)=1; 123 aq = 1; 124 end 125 end 126 end 127 end 128 129 % 3.row 对打勾的列中含独立0元素的未打勾的行打勾 130 for j =1:a 131 if(y(j)==1) 132 for i = 1:a 133 if(C(i,j)==0 && x(i)==0 && r(i,j)==2) 134 x(i)=1; 135 aq = 1; 136 end 137 end 138 end 139 end 140 end 141 142 %没有打勾的行和有打勾的列画线，这就是覆盖所有0元素的最少直线数 143 line = 0; %存放划线的数量 144 h = 0; %行划线后上移h行 145 l = 0; %列划线后前移l列 146 %另一种画线 147 for i = 1:a 148 delr = i-h; 149 if(x(i)==0) 150 A(delr,:)=[]; 151 h = h+1; 152 del_hang(i)= i; %得到被覆盖的行数 153 line = line+1; 154 end 155 delc=i-l; 156 if (y(i)==1) 157 A(:,delc)=[]; 158 l = l+1; 159 del_lie(i) = i; %得到被覆盖的列数 160 line = line+1; 161 end 162 end 163 N_min = min(min(A)); %找到最小的数
测试
C = [12 7 9 7 9;
8 9 6 6 6;
7 17 12 14 9;
15 14 6 6 10;
4 10 7 10 9];
A = hungary(C);
参考链接：http://www.cnblogs.com/chenyg32/
