遗传算法求解混合流水车间调度问题（HFSP）二：算法实现一
遗传算法求解混合流水车间调度问题（HFSP）二：算法实现一
遗传算法的设计
- 编码：对工件进行优先级编码，编码越小，优先级越高。
- 解码：按照工件优先级进行生产，求出整体完工时间。
- 目标函数值：整体完工时间。
- 适应度值：目标函数越小，适应度值越大。
- 选择：按照轮盘赌方法进行选择。
- 交叉：按照交叉概率，选择两个父代个体（父1和父2），交叉生成两个子代个体（子1和子2）。生成一个随机的交换空间，交换空间内，子1基因来自于父2，子2基因来自于父1；交换空间外，子1基因来自于父1，子2基因来自于父2。注意任意两个零件优先级不能相等。
- 变异：按照变异概率，随机交换两个基因。
- 终止条件：迭代固定次数。
计算结果
在本例中，有6个工件，有3道工序，每道工序有2台机器，下面是执行结果：
最优序列:
3 4 6 2 1 5
上图中，第1、2行是第1工序的2台设备，第3、4行是第2工序的2台设备，第5、6行是第3工序的两台设备，纵轴代表时间。按照最优序列[ 3 4 6 2 1 5]赋予每个零件优先级，一共用时25.
主函数
首先定义问题的参数：
piecetime = [2 4 6; ... % 设备加工时间
4 9 2; 4 2 8; 9 5 6; 5 2 7; 94 3];
equsize = [2 2 2]; % 每个工序设备数目
piecesize = size(piecetime, 1); % 工件数目
prosize = size(piecetime, 2); % 工序数目
在本例中，一共有6个设备，3道工序，设备必须按照1-2-3的工序进行加工，每个工序有2台机器。一共有6个工件。piecetime中行代表工件，列代表工序，内容是加工时间，比如第1行第1列，表示工件1在工序1加工时间为2.
定义遗传算法的参数：
popsize = 20; % 种群规模
cr = 0.6; % 交叉概率
mr = 0.05; % 变异概率
maxgen = 100; % 迭代次数
进行迭代：
pop = initpop(popsize, piecesize);
for gen = 1:maxgen
[objvalue, ptr,per] = calobjvalue(pop, piecetime, equsize);
fitness= calfitness(objvalue); % 计算适应度值
pop =selection(pop, fitness); % 选择
pop =crossover(pop, cr); % 交叉
pop =mutation(pop, mr); % 变异
end
主函数全部代码如下：
function main()
piecetime = [2 4 6; ... % 设备加工时间
4 9 2; 4 2 8; 9 56; 5 2 7; 9 4 3];
equsize = [2 2 2]; % 每个工序设备数目
piecesize = size(piecetime, 1); % 工件数目
prosize = size(piecetime, 2); % 工序数目
popsize = 20; % 种群规模
cr = 0.6; % 交叉概率
mr = 0.05; % 变异概率
maxgen = 100; % 迭代次数
bestobjvalue = zeros(1, maxgen);
bestpop = zeros(maxgen, piecesize);
avgobjvalue = zeros(1, maxgen);
bestptr = cell(1, maxgen);
bestper = cell(1, maxgen);
pop = initpop(popsize, piecesize);
for gen = 1:maxgen
[objvalue, ptr,per] = calobjvalue(pop, piecetime, equsize);
[bobjvalue,bindex] = min(objvalue);
bestobjvalue(1,gen) = bobjvalue;
bestpop(gen, :) =pop(bindex, :);
avgobjvalue(1,gen) = sum(objvalue) / popsize;
bestptr{1, gen} =ptr{1, bindex};
bestper{1, gen} =per{1, bindex};
fitness= calfitness(objvalue); % 计算适应度值
pop =selection(pop, fitness); % 选择
pop =crossover(pop, cr); % 交叉
pop =mutation(pop, mr); % 变异
end
[~, finalindex] = min(bestobjvalue);
finalptr = bestptr{1, finalindex};
finalper = bestper{1, finalindex};
fprintf("最优序列:\n");
disp(bestpop(finalindex, :));
gantt = makegantt(finalptr, finalper, equsize);
figure(1);
imagesc(gantt);
colorbar;
title("加工流程图");
figure(2);
plot(1:maxgen, bestobjvalue);
title("最优时间变化图");
xlabel("代数"); ylabel("最优时间");
figure(3);
plot(1:maxgen, avgobjvalue);
title("平均时间变化图");
xlabel("代数"); ylabel("平均时间");
end
计算目标函数值函数
在第一道工序中，所有的工件同时等待被加工，则按照优先级进行加工；在第二道和之后的工序中，由于上一道工序中工件完工时间不同，上一道工序先加工完的工件先进行本工序加工。
function [objvalue, ptr, per] = calobjvalue(pop, piecetime,equsize)
% 计算目标函数值
% pop input 种群
% piecetime input 工件加工时间
% equsize input 每个工序设备数量
% objvalue output 目标函数值（加工时间）
% ptr output 工件加工时间记录，cell
% per output 工件加工设备记录，cell
[popsize, piecesize] = size(pop);
prosize = size(equsize, 2);
objvalue = zeros(popsize, 1);
ptr = cell(1, popsize);
per = cell(1, popsize);
for i = 1:popsize
pieceweight =pop(i, :);
% 设备状态序列
% [工序1设备1工序1设备2 工序2设备1 工序2设备2 ……]
% 记录当前设备使用结束时间，默认为0表示未开始
equstu = zeros(1,sum(equsize));
% 对设备状态序列的工序分隔符
% 大于等于当前设备最小值的索引是当前设备所处的工序
% [2 35] 工序1有2台设备工序2有1台设备 工序3有2台设备
prosep =cumsum(equsize);
% 工件时间记录，记录每个工件每个工序的开始时间和结束时间
% 行表示工件，相邻两列表示开始加工时间和停止加工时间
% [1 2 2 3; 4 5 67]
% 表示工件1第1工序加工时间为1-2，第2工序加工时间为2-3
% 工件2第1工序加工时间为4-5，第2工序加工时间为6-7
piecetimerecord =zeros(piecesize, prosize*2);
% 工件设备记录，记录每个工件在工序中的加工设备
% 行数表示工件，列表示该零件在每个工序加工设备
% [1 2; 2 1]
% 表示工件1在第1工序加工设备为1，第2工序加工设备为2
% 工件2在第1工序加工设备为2，第2工序加工设备为1
pieceequrecord =zeros(piecesize, prosize);
% 对每一道工序
% 如果是第1道工序，对工件按优先级排序
% 其余工序上上一道工序完工时间对工件排序
% 对排序后的每一件工件
% 对该工序中可用机器按使用结束时间排序
% 使用使用结束时间最小的机器
% 加工开始时间为max{设备使用结束时间,零件上一工序完工时间}
% 加工结束时间=加工开始时间+加工时间
% 更新各个状态和记录矩阵
for pro =1:prosize
if(pro == 1)
[~,piecelist] = sort(pieceweight);
else
tempendtime = piecetimerecord(:, (pro-1)*2);
tempendtime = tempendtime';
[~,piecelist] = sort(tempendtime);
end
for pieceindex= 1:length(piecelist)
piece =piecelist(pieceindex);
equtimelist = equstu(prosep(pro)-equsize(pro)+1:prosep(pro));
[equtime,equlist] = sort(equtimelist);
equ =equlist(1);
if pro ==1
piecestarttime = 0;
else
piecestarttime = piecetimerecord(piece, pro*2-2);
end
starttime= max(equtime(1), piecestarttime) + 1;
endtime =starttime + piecetime(piece, pro) - 1;
equstuindex = prosep(pro)-equsize(pro)+equ;
equstu(equstuindex) = endtime;
piecetimerecord(piece, pro*2-1) = starttime;
piecetimerecord(piece, pro*2) = endtime;
pieceequrecord(piece, pro) = equ;
end
end
objvalue(i, 1) =max(max(piecetimerecord));
ptr{1, i} =piecetimerecord;
per{1, i} =pieceequrecord;
end
end
选择函数
使用轮盘赌方法进行选择：
function spop = selection(pop, fitness)
% 轮盘赌选择
% pop input 种群
% fitness input 适应度值
% spop output 选择后生成的种群
[popsize, piecesize] = size(pop);
spop = zeros(popsize, piecesize);
sumfit = sum(fitness);
fitness = fitness ./ sumfit;
fitness = cumsum(fitness);
r = rand(1, popsize);
r = sort(r);
j = 1;
for i = 1:popsize
while fitness(j)< r(i)
j = j + 1;
end
spop(i, :) =pop(j, :);
end
% 由于上面轮盘赌方法特殊性，一个个体在相邻位置多次重复，故随机排序
rr = randperm(popsize);
spop(:, :) = spop(rr, :);
end
交叉函数
按照交叉概率，选择两个父代个体（父1和父2），交叉生成两个子代个体（子1和子2）。生成一个随机的交换空间，交换空间内，子1基因来自于父2，子2基因来自于父1；交换空间外，子1基因来自于父1，子2基因来自于父2。注意任意两个零件优先级不能相等。
function cpop = crossover(pop, cr)
% 交叉
% pop input 种群
% cr input 交叉概率
% cpop output 交叉后种群
[popsize, piecesize] = size(pop);
cpop = pop;
if mod(popsize,2) ~= 0
nn = popsize - 1;
else
nn = popsize;
end
% 父代father mother, 子代son daughter
% 在rl:ru中，son继承mother，daughter继承father
% 其余位置son继承father，daughter继承mother
for i = 1:2:nn
if rand > cr
continue;
end
[rl, ru] =makerlru(piecesize);
father = pop(i, :);
mother = pop(i+1, :);
if father == mother
continue;
end
son = zeros(1, piecesize);
daughter = zeros(1,piecesize);
son(rl:ru) = mother(rl:ru);
daughter(rl:ru) =father(rl:ru);
j = 1;
for k = 1:piecesize
if k >= rl &&k <= ru
continue;
end
while ~isempty(find(son== father(j), 1))
j = j + 1;
end
son(k) = father(j);
end
j = 1;
for k = 1:piecesize
if k >= rl &&k <= ru
continue;
end
while~isempty(find(daughter == mother(j), 1))
j = j + 1;
end
daughter(k) = mother(j);
end
cpop(i, :) = son;
cpop(i+1, :) = daughter;
end
end
变异函数
随机交换染色体中两个位置的基因：
function mpop = mutation(pop, mr)
% 变异，交换两个随即位置的基因
% pop input 种群
% mr input 变异概率
% mpop output 变异后种群
[popsize, piecesize] = size(pop);
mpop = pop;
for i = 1:popsize
if rand > mr
continue;
end
r1 = randi(piecesize);
r2 = randi(piecesize);
temp = mpop(i, r1);
mpop(i, r1) = mpop(i, r2);
mpop(i, r2) = temp;
end
end