flowshop.m

离散算法

clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                 参数部分                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
runTime = 1;
% －－－多目标工作流问题参数－－－－－－－－－－－－－－－－－－－－－－－－－
funCount=2;       % 目标函数个数
jobCount = 100;
machCount = 10;
proFileName = ['tai', num2str(jobCount), '_', num2str(machCount)];
load(proFileName) % machCount jobCount proCount all_tProcess all_otherParetoOptimalFun all_otherPOFunCounts all_fBound;
selProInds = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; % 从中选出的 测试问题 下标  1, 2, 3, 4, 5, 6, 7, 8, 9, 10
selProCount = size(selProInds, 2); % 从中选出的 测试问题 个数

% －－－粒子群算法参数－－－－－－－－－－－－－－－－－－－－－－－－－
popSize = 50;     % 种群规模250
evaluations =15000 ; % 最大迭代次数20000
endEval = floor(0.9*evaluations); 
w_start = 0.9;   %Initial inertia weight's value
w_end = 0.3;       %Final inertia weight
w_varyfor = floor(0.7*evaluations); % evaluations - 1;Weight change step. Defines total number of evaluations for which weight is changed.
w_now = w_start;
inertdec = (w_start-w_end)/w_varyfor; % Inertia weight's change per iteration

c1_start = 0.9; % 最初的变异概率
c1_end = 0.15;
c1_varyfor = floor(0.7 * evaluations);
c1_now = c1_start;
c1_dec = (c1_start - c1_end) / c1_varyfor;


% 产生新解比原来差时，仍然更新为新解的概率，递减
upd_start = 0.5;   %Initial inertia weight's value
upd_end = 0.1;       %Final inertia weight
upd_varyfor = floor(0.7*evaluations); % Weight change step. Defines total number of evaluations for which weight is changed.
upd_now = upd_start;
upd_dec = (upd_start-upd_end)/upd_varyfor; % Inertia weight's change per iteration

replaceRate = 0.5;

% －－－适应度函数的参数－－－－－－－－－－－－－－－－－－－－－－－－－
o_shares = [0.01, 1]; % 各个目标空间上的共享半径，（1 × funCount）维）， ！！！此处固定为 2个目标
% ！！！以下两个参数只作为参数传入适应度的计算，实际适应度函数的计算公式固定，没用到该参数，可进一步研究
a_nonlinear = 2;    % 非线性结合被优越数domiCount的指数参数  
b_nonlinear = 1;    % 非线性结合邻近密度neighDensity的指数参数

% －－－数据结构中用到的参数－－－－－－－－－－－－－－－－－－－－－－－－－
% －－粒子的存储结构
% － [X1,...,Xn,   f1,...fk,    domiCount, neighDensity, Fitness,  pBInd      ]
% － [决策变量值， 目标函数值， 被优越数， 邻近密度，   适应度值， 在历史最优位置数组中的横坐标（粒子编号）]
% － －其中 pBInd 的初始编号为初始种群按从优到劣的顺序排列后，给各个粒子的编号，便于跟踪最优粒子的变化情况
parInfSize = jobCount + funCount + 3 + 1 ;% 表示粒子信息的一维数组长度
funInd = jobCount + 1;
dCInd = funInd + funCount;  % 被优越数值在种群数组中的第几列 
nDInd = dCInd + 1;          % 邻近密度值在种群数组中的第几列 
fitInd = nDInd + 1;         % 适应度值在种群数组中的第几列 
pBInd = fitInd + 1;          % 记录粒子编号的变量在第几列，用于后面pBest排序

gBCount=0;  % 全局最优解的个数
pOCount = 0; % 历代Pareto最优解的个数
gBest = [];                 % 记录适应度最优的最优解
ParetoOptimal = [];         % 记录历代搜索到的Pareto最优解

% －－－输出图中用到的参数－－－－－－－－－－－－－－－－－－－－－－－－－
% －－
figDiv = 10;   % 坐标分为几段
figF1Min = 0;
figF1Max = 0;
figF1Len = 0;
figF1Step = figF1Len / figDiv;

figF2Min = 0;
figF2Max = 0;
figF2Len = 0;
figF2Step = figF2Len / figDiv; 

midShowTime = 2; % 中间图输出的个数
midShowStep = floor(evaluations / midShowTime);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                 参数部分      结束                    %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
proParetoOptimalFun = -ones(40, funCount + 1,selProCount); % 初始为每个问题的 PTrue 个数 上限为 40 个
proPOFunCount =  -ones(selProCount, 1);
proMyPSOPTrueFunCounts = -ones(selProCount, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%          开始    各个 测试问题计算                       %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beforeProTimeStr = '（无）';
for (iSelProInd = 1: selProCount)		
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%                读入 测试问题 部分                        %%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    selProNum = selProInds(1, iSelProInd);  % 测问题编号
	tProcess = all_tProcess(:, :, selProNum); % 测试问题的 处理时间
	otherPOFunCounts= all_otherPOFunCounts(selProNum) ; % 测试问题的 其他算法得到的最优解 个数
    otherParetoOptimalFun = [];
	otherParetoOptimalFun = all_otherParetoOptimalFun(1:otherPOFunCounts,:,selProNum); % 测试问题的 其他算法得到的最优解
    f1OPOMin = min(otherParetoOptimalFun(:, 1)); % 测试问题 其他人全部最优解的边界，用于最终图输出
    f2OPOMin = min(otherParetoOptimalFun(:, 2));
    f1OPOMax = max(otherParetoOptimalFun(:, 1)); 
    f2OPOMax = max(otherParetoOptimalFun(:, 2));
    
	mid_f1Min = all_fBound(selProNum, 1); % 测试问题的边界值（包含一个随机序列和其他人全部最优解的边界），用于初始图和中间图输出
	mid_f2Min = all_fBound(selProNum, 2);
	mid_f1Max = all_fBound(selProNum, 3);
	mid_f2Max = all_fBound(selProNum, 4);   	            
    proDes = strcat(' Taillard -  ',num2str(jobCount),'/',num2str(machCount),', Problem ',num2str(selProNum));
    proShotDes = strcat('tai',num2str(jobCount),'_',num2str(machCount),' - ',num2str(selProNum));
    strcat(proDes, '开始：',DATESTR(NOW), '上个函数：', beforeProTimeStr)   
    beforeProTimeStr = DATESTR(NOW);
    % －－10.1 得到初始图和中间图输出的边界  
	mid_figF1Step = round((mid_f1Max - mid_f1Min) / (figDiv - 2 ));
	mid_figF1Min = mid_f1Min - mid_figF1Step;
	mid_figF1Max = mid_f1Max + mid_figF1Step/4 ; % f1Min + figF1Step * figDiv;
	mid_figF1Len = mid_figF1Max - mid_figF1Min;
	mid_figF2Step = round((mid_f2Max - mid_f2Min) / (figDiv - 2 ));
	mid_figF2Min = mid_f2Min - mid_figF2Step;
	mid_figF2Max = mid_f2Max + mid_figF2Step/4;
	mid_figF2Len = mid_figF2Max - mid_figF2Min;
	mid_figBounds = [mid_figF1Min, mid_figF1Max, mid_figF1Step;
                 mid_figF2Min, mid_figF2Max, mid_figF2Step];
% % 	%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%%
% % 	%       得到图输出的边界
% % 	figF1Step = round((f1Max - f1Min) / (figDiv - 2 ));
% % 	figF1Min = f1Min - figF1Step;
% % 	figF1Max = f1Max + figF1Step/4 ; % f1Min + figF1Step * figDiv;
% % 	figF1Len = figF1Max - figF1Min;
% % 	figF2Step = round((f2Max - f2Min) / (figDiv - 2 ));
% % 	figF2Min = f2Min - figF2Step;
% % 	figF2Max = f2Max + figF2Step/4;
% % 	figF2Len = figF2Max - figF2Min;
% % 	figBounds = [figF1Min, figF1Max, figF1Step;
% %                  figF2Min, figF2Max, figF2Step];
             
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%                开始多轮 进化                             %%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    runTimeParetoOptimals = []; % 记录每一轮计算得到的 最优解
	beforeTimeStr = '（无）';
	for (iRunTime = 1: runTime)	
% %         strcat('新一轮开始( ',proDes, '： 第',num2str(runTime), ' - ',num2str(iRunTime),'轮)', ...
% %             '－－现在时间：',DATESTR(NOW), '－上一轮开始时间：', beforeTimeStr)
        beforeTimeStr = DATESTR(NOW);
        w_now = w_start;
        c1_now = c1_start;
        upd_now = upd_start;
	%     iRunTime
		%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－%%
		%       初始化种群，计算函数值和适应度值
		initialPop = ones(popSize, pBInd);
		for i=1:popSize
            initialPop(i,1:jobCount) = randperm(jobCount); %[2,3,1,4,5];
            curJobSeq = initialPop(i, 1:jobCount);		
            % 计算完成时间和函数值  ！！！暂时去掉 完成时间 的返回