flowshop.m

离散算法

            [initialPop(i, funInd: funInd + 1)] = get2FunValues(initialPop(i, 1:jobCount), tProcess); % tComplete(i, 1:machCount, 1: jobCount),  , duedate
		end
		% 计算种群适应值
		[initialPop(:, fitInd: fitInd), initialPop(:, dCInd: dCInd), initialPop(:, nDInd: nDInd)] = ...
            getFitnessValues(initialPop(:, funInd: funInd + 1), o_shares, a_nonlinear, b_nonlinear);
		%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%%
		%       按适应值从高到低排序，得到新的种群stPop，并编号
		[stFit, idx] = sort(initialPop(:, fitInd)); % 按适应度这一列升序排序
		stPop = initialPop(idx, :);                 % 使用适应度列的游标顺序重排列initialPop
		for (i = 1: popSize)
            stPop(i, pBInd) = i; %记录粒子编号，代表了初始种群中按优到劣的序号，用于后面pBest排序
		end
		%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%%
		%       更新局部最优值和全局最优(多个)和历代Pareto最优解
		pBest = stPop;
		[gBest gBCount]=getStBests(stPop, fitInd, jobCount, dCInd, nDInd, fitInd, pBInd); % 返回的 gBest中，除了粒子编号统一改为 -1 以示区别，其他的数据不变

        [stFit, idx] = sort(initialPop(:, dCInd)); % 按 被优越数 这一列升序排序
        stDCPop = initialPop(idx, :);              % 用 被优越数 列的游标顺序重排列 allPop
		[ParetoOptimal pOCount]=getStBests(stDCPop, dCInd, jobCount, dCInd, nDInd, fitInd, pBInd); % 返回的 ParetoOptimal，除了粒子编号统一改为 -1 以示区别，其他的数据不变
        
% % % 		%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%%
% % % 		%       输出初始种群，并存为 bmp 文件       
% % %         output2PopFigure(10, proDes, stPop(:, funInd: funInd+1), otherParetoOptimalFun, 1,popSize,0, mid_figBounds); % otherParetoOptimal
% %         filename = strcat('initial_',proShotDes,'_p',num2str(popSize),'_',num2str(iRunTime),'_e',num2str(evaluations)); % 
% %         print('-depsc2', filename);
		
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		%%                  开始 每一轮的 PSO进化                   %%
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
		allPareto = 0;
        findAllD0 = 0;
		eval = 1;
        beforeEvalTimeStr = DATESTR(NOW);
		while (eval <= evaluations) && (allPareto == 0)
		%      eval		
		%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%%          
		%       初始新种群的newPop %%%为，保留作业序列、对应的目标值和粒子编号（一直不变）
            newPop = stPop;
		%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%%    
		%       更新当前的惯性权值 w
            if (eval <= w_varyfor + 1) & (eval > 1)
                w_now = w_now - inertdec; %Change inertia weight ＃＃ inertdec = (w_start-w_end)/w_varyfor;
            end
% %             w_now = 1;
            % －－1.1 计算该粒子的两个调整因子 （！！！待定，看是否要变化） 
            if (eval <= c1_varyfor + 1) & (eval > 1)
                c1_now = c1_now - c1_dec; % 对应个体最优
            end
            c2_now = 1 - c1_now; % 对应全局最优
% % 	        c1_now = 0.5;
% % 	        c2_now = 0.5;
            
            %－－2.1 计算当前的更新概率 upd_now
            if (eval <= upd_varyfor + 1) & (eval > 1)
                upd_now = upd_now - upd_dec; %Change inertia weight ＃＃ upd_dec = (upd_start-upd_end)/upd_varyfor;
            end
% % % 		    upd_now = 0.5;
		    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%                 更新种群中的每个粒子的位置               %%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
            for (particle = 1: popSize)
			    %%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%% 
                %       1.计算新位置
			    % －－1.1 根据惯性权值 w_now，变异粒子的序列
                [newPop(particle, 1: jobCount)] = getMutation(stPop(particle, 1: jobCount), w_now);                               
                % －－1.2 个体最优值与当前新位置交叉－－－－－－－－－－－－
                curPBInd = newPop(particle, pBInd);
                [pBested] = getCrossover(newPop(particle, 1: jobCount), pBest(curPBInd, 1: jobCount), c1_now);
                % －－1.3（方案1）随机选择一个 gBest，与当前新位置交叉－－－－－－－－－－－－
		        if (gBCount == 1)
		            gBi = 1;
		        else  
		            gBi = ceil(rand(1) * (gBCount));
		        end
		        gBested = -ones(1, jobCount + funCount);
		        [gBested(1, 1: jobCount)] = getCrossover(pBested, gBest(gBi, 1: jobCount), c2_now); 
                % －－1.3（方案2）随机选择一个 ParetoOptimal ，与当前新位置交叉－－－－－－－－－－－－
% %                 if (pOCount == 1)
% %                     pOi = 1;
% %                 else  
% %                     pOi = ceil(rand(1) * (pOCount));
% %                 end
% %                 gBested = -ones(1, jobCount + funCount);
% %                 [gBested(1, 1: jobCount)] = getCrossover(pBested, ParetoOptimal(pOi, 1: jobCount), c2_now);
				%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%% 
				%       2. 更新位置        
                %－－ （方案1）比较新解 gBested 与 原解 stPop 的优劣，并更新
                [gBested(1, funInd: funInd + 1)] = get2FunValues(gBested(1, 1:jobCount), tProcess); % tComplete(curPBInd, 1:machCount, 1: jobCount),  , duedate
                if (gBested(1, funInd) <= stPop(particle, funInd)) && (gBested(1, funInd + 1) <= stPop(particle, funInd + 1)) && ...
                   ((gBested(1, funInd) + gBested(1, funInd + 1)) < (stPop(particle, funInd) + stPop(particle, funInd + 1)))
                    newPop(particle, 1: jobCount) = gBested(1, 1: jobCount);
                else 
                    if ((rand(1) < upd_now)) % 如果新解较差且随机数小于更新概率，则更新，否则 newPop 不变
                        newPop(particle, 1: jobCount) = gBested(1, 1: jobCount);
                    end       
                end
                 %－－（方案2）直接替换
% % % % 				newPop(particle, 1: jobCount) = gBested(1, 1: jobCount);
				%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%% 
				%       3.计算粒子的函数值－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%    
                [newPop(particle, funInd: funInd + 1)] = get2FunValues(newPop(particle, 1:jobCount), tProcess);    %  tComplete(curPBInd, 1:machCount, 1: jobCount),  duedate
                %newPop
				%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%% 
				%       4.更新局部最优
                 if (newPop(particle, funInd) <= pBest(curPBInd, funInd)) && (newPop(particle, funInd + 1) <= pBest(curPBInd, funInd + 1)) && ...
                   ((newPop(particle, funInd) + newPop(particle, funInd + 1))<(pBest(curPBInd, funInd) + pBest(curPBInd, funInd + 1)))
                    pBest(curPBInd, :) = newPop(particle,:);
                 end
				%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%% 
				%       5.求全局最优解： 合并新种群和前一代全局最优，将适应度最高的一组 非劣解（！！去掉劣解） 作为新的全局最优解
                allPop=[newPop; gBest]; % allPop为当前的新种群和前一代的全局最优个体的并集
                allPop(:, dCInd: fitInd) = -1;
                [allPop(:, fitInd: fitInd), allPop(:, dCInd: dCInd), allPop(:, nDInd: nDInd)] = ...
                    getFitnessValues(allPop(:, funInd: funInd + 1), o_shares, a_nonlinear, b_nonlinear); % 求合并种群的适应度 
                [stFit, idx] = sort(allPop(:, fitInd)); % 按适应度这一列升序排序
                stAllPop = allPop(idx, :);              % 用适应度列的游标顺序重排列 allPop		
                % 更新全局最优(多个)
			    % % %     gBCount_before= gBCount;  % 用于后面替换解时，知道替换的范围（！！！可以直接读allPop的大小）
                [gBest gBCount] = getStBests(stAllPop, fitInd, jobCount, dCInd, nDInd, fitInd, pBInd);
            end % end for (particle = 1: popSize)   
		    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%        结束     更新种群中的每个粒子的位置               %%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
			%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%% 
			%      7. 将种群中适应度最小的gBCount_before个粒子替换成新的全局最优，！！！！先不作替换
                %gBCount_before
			% % % %     for (repPop = popSize+1 : popSize + gBCount_before)
			% % % %         repRan = 1;%rand(1);
			% % % %         if (repRan >= replaceRate) % 只有大于替换概率replaceRate时，才替换（为了最后的多样性，可进一步考虑变replaceRate）
			% % % %             badInd = stAllPop(repPop, pBInd);
			% % % %             if (badInd >= 1) && (badInd <= popSize) % 在正常范围，参见getStGBest中对全局最优解的粒子号定义（为-1）
			% % % %                 if (gBCount==1)
			% % % %                     selGBInd=1;
			% % % %                 else  
			% % % %                     selGBInd=ceil(rand(1)*(gBCount-1));
			% % % %                 end
			% % % %                 %strcat('Replace: ', num2str(badInd), ' to ', num2str(selGBInd))
			% % % %                 newPop(badInd, 1: jobCount + funCount) = gBest(selGBInd, 1: jobCount + funCount);  
			% % % %             end
			% % % %         end % end of if (repRan >= replaceRate)
			% % % %     end
            
			%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%% 
			%      8.计算粒子在当前种群newPop中的适应度
            [newPop(:, fitInd: fitInd), newPop(:, dCInd: dCInd), newPop(:, nDInd: nDInd)] = ...
                getFitnessValues(newPop(:, funInd: funInd + 1), o_shares, a_nonlinear, b_nonlinear); 
			%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%% 
			%     更新历代Pareto最优解集： 合并新种群和前一代Pareto最优，将被优越数为0的一组 非劣解 作为新的 Pareto最优
			% % %     [ParetoOptimal pOCount] = getBests(stPop, ParetoOptimal, dCInd)
            allPop=[newPop; ParetoOptimal]; % allPop为当前的新种群和前一代的全局最优个体的并集
            allPop(:, dCInd: fitInd) = -1;
            [allPop(:, fitInd: fitInd), allPop(:, dCInd: dCInd), allPop(:, nDInd: nDInd)] = ...
                getFitnessValues(allPop(:, funInd: funInd + 1), o_shares, a_nonlinear, b_nonlinear); % 求合并种群的适应度 
            [stFit, idx] = sort(allPop(:, dCInd)); % 按 被优越数 这一列升序排序