利用禁忌算法求解图节点着色问题

% 利用禁忌算法解决“图节点着色问题” % 算法参见《现代优化计算方法》P70 clear global K; % 划分集合数 K = 100; global h; % 禁忌表长度 h = 10; fstar = 0; % 目标值下界 Nbmax = 100; % 两目标值未改进的最大迭代次数 global NbNode; NbNode = 1000; % 节点数 global map; %节点图既可通过手写也可通过随机函数获得。节点图是对称矩阵，1表示两节点有连接，0表示不连接，对角线元素为0。 %map = [0,1,1,1,0; 1,0,0,0,1; 1,0,0,1,0; 1,0,1,0,0; 0,1,0,0,0]; %map = NodeMap(NbNode); load Map.mat; %S0 = [ 1, 3, 2, 1, 2 ]; S0 = RandS0( NbNode, K ); % 随机产生节点划分，产生NbNode个随机数，随机数在1-K间。注意：初始值的划分集合数与K值要相同 %load S0.mat; Nbiter = 0; Nbbestiter = 0; Sbestsol = S0; S = S0; SValue = f( S ); theBestValue = SValue; A = Aspir( SValue ); H = 0; %初始化禁忌表 NumH = 0; while ( (SValue > fstar) && (Nbiter-Nbbestiter<Nbmax) ) Nbiter = Nbiter + 1; % 求S的邻域。输出：由S到任一邻域的变化特征；任一邻域的目标值；邻域数目。 % VstarChange（1:3）：变量1，变量变化的节点序号；变量2，改变前节点对应的集合；变量3，改变后节点对应的集合。 [VstarChange, VstarValue, Num] = N( S, SValue ); % search the best solution of Vstar SstarValue = 1e5; % As soon as biger index = 1; while( 1 ) [ C, I ] = min( VstarValue ); i = I; bTabu = false; for j=1:1:NumH % 若 Vstar 被禁 if ( VstarChange(i, 1) == H(j, 1) ) %if ( VstarChange(i, 1) == H(j, 1) && VstarChange(i, 2) == H(j, 2) && VstarChange(i, 3) == H(j, 3) ) %若Vstar被禁且不被特赦 if VstarValue(i) > A bTabu = true; break; end end end if bTabu == true VstarValue(i) = inf; continue; end SstarValue = C; index = i; break; end %Sstar由S得到。 Sstar = S; Sstar( VstarChange(index, 1) ) = VstarChange(index, 3); %记录每次迭代的划分集的变化 VAccuChange(Nbiter, :) = VstarChange(index, : ); % 计算特赦值 if SstarValue <= Aspir( SValue ) % if ( f(Vstar) ) <= Aspir( f(S) ) A = Aspir( SstarValue ); elseif SValue <= Aspir( SstarValue ) A = Aspir( SValue ); end %更新禁忌表 [H, NumH] = updateH( H, NumH, VstarChange(index,:) ); %得到最优结果 if SstarValue < theBestValue Sbestsol = Sstar; Nbbestiter = Nbiter; theBestValue = SstarValue; end S = Sstar; SValue = SstarValue; if(mod(Nbiter,10) == 0 ) fprintf('The Num of iter is %d, currently the best Value is %d, the whole best value is %d\n', Nbiter, SstarValue, theBestValue); end end %计算并打印集合划分结果 Nb(1:K) = 0; k = 0; for i=1:1:NbNode j = Sbestsol(i); Nb(j) = Nb(j) + 1; V( j, Nb(j) ) = i; end V %最终的划分情况 theBestValue %最终的评价值