当前位置:主页 > 软件编程 > C代码 >

matlab鸟群算法求解车间调度问题详解及实现源码

时间:2022-06-23 11:39:25 | 栏目:C代码 | 点击:

一、车间调度简介

1 车间调度定义

车间调度是指根据产品制造的合理需求分配加工车间顺序,从而达到合理利用产品制造资源、提高企业经济效益的目的。车间调度问题从数学上可以描述为有n个待加工的零件要在m台机器上加工。问题需要满足的条件包括每个零件的各道工序使用每台机器不多于1次,每个零件都按照一定的顺序进行加工。

2 传统作业车间调度

传统作业车间带调度实例

在这里插入图片描述

有若干工件,每个工件有若干工序,有多个加工机器,但是每道工序只能在一台机器上加工。对应到上面表格中的实例就是,两个工件,工件J1有三道工序,工序Q11只能在M3上加工,加工时间是5小时。
约束是对于一个工件来说,工序的相对顺序不能变。O11->O12->O13。每时刻,每个工件只能在一台机器上加工;每个机器上只能有一个工件。
调度的任务则是安排出工序的加工顺序,加工顺序确定了,因为每道工序只有一台机器可用,加工的机器也就确定了。
调度的目的是总的完工时间最短(也可以是其他目标)。举个例子,比如确定了O21->O22->O11->O23->O12->O13的加工顺序之后,我们就可以根据加工机器的约束,计算出总的加工时间。
M2加工O21消耗6小时,工件J2当前加工时间6小时。
M1加工O22消耗9小时,工件J2当前加工时间6+9=15小时。
M3加工O11消耗5小时,工件J1当前加工时间5小时。
M4加工O23消耗7小时,工件J2加工时间15+7=22小时。
M1加工O12消耗11小时,但是要等M1加工完O22之后才开始加工O12,所以工件J1的当前加工时间为max(5,9)+11=20小时。
M5加工O13消耗8小时,工件J2加工时间20+8=28小时。
总的完工时间就是max(22,28)=28小时。

3 柔性作业车间调度

柔性作业车间带调度实例(参考自高亮老师论文
《改进遗传算法求解柔性作业车间调度问题》——机械工程学报)

在这里插入图片描述

相比于传统作业车间调度,柔性作业车间调度放宽了对加工机器的约束,更符合现实生产情况,每个工序可选加工机器变成了多个,可以由多个加工机器中的一个加工。比如上表中的实例,J1的O12工序可以选择M2和M4加工,加工时间分别是8小时和4小时,但是并不一定选择M4加工,最后得出来的总的完工时间就更短,所以,需要调度算法求解优化。

相比于传统作业车间,柔性车间作业调度的调度任务不仅要确定工序的加工顺序,而且需要确定每道工序的机器分配。比如,确定了O21->O22->O11->O23->O12->O13的加工顺序,我们并不能相应工序的加工机器,所以还应该确定对应的[M1、M3、M5]->[M1、M2、M3]->[M1、M2、M3、M4、M5]->[M2、M3、M4、M5]->[M2、M4]->[M1、M3、M4、M5]的机器组合。调度的目的还是总的完工时间最短(也可以是其他目标,比如机器最大负荷最短、总的机器负荷最短)

二、蝴蝶优化算法(MBO)简介

1 介绍

蝴蝶优化算法(butterfly optimization algorithm, BOA)是Arora 等人于2019年提出的一种元启发式智能算法。该算法受到了蝴蝶觅食和交配行为的启发,蝴蝶接收/感知并分析空气中的气味,以确定食物来源/交配伙伴的潜在方向。
蝴蝶利用它们的嗅觉、视觉、味觉、触觉和听觉来寻找食物和伴侣,这些感觉也有助于它们从一个地方迁徙到另一个地方,逃离捕食者并在合适的地方产卵。在所有感觉中,嗅觉是最重要的,它帮助蝴蝶寻找食物(通常是花蜜)。蝴蝶的嗅觉感受器分散在蝴蝶的身体部位,如触角、腿、触须等。这些感受器实际上是蝴蝶体表的神经细胞,被称为化学感受器。它引导蝴蝶寻找最佳的交配对象,以延续强大的遗传基因。雄性蝴蝶能够通过信息素识别雌性蝴蝶,信息素是雌性蝴蝶发出的气味分泌物,会引起特定的反应。
通过观察,发现蝴蝶对这些来源的位置有非常准确的判断。此外,它们可以辨识出不同的香味,并感知它们的强度。蝴蝶会产生与其适应度相关的某种强度的香味,即当蝴蝶从一个位置移动到另一个位置时,它的适应度会相应地变化。当蝴蝶感觉到另一只蝴蝶在这个区域散发出更多的香味时,就会去靠近,这个阶段被称为全局搜索。另外一种情况,当蝴蝶不能感知大于它自己的香味时,它会随机移动,这个阶段称为局部搜索。

2 香味

为了理解BOA中的香味是如何计算的,首先需要理解,像气味、声音、光、温度等这样的模态是如何计算的。感知、处理这些模态需要知道三个重要的术语:感觉模态C、刺激强度I和幂指数a。在感觉模态中,感觉意味着测量能量的形式并以类似方式对其进行处理,而模态是指传感器使用的原始输入。不同的形态可以是气味,声音,光线,温度,在BOA中,模态是香味。I是物理刺激的大小。在BOA中,I与蝴蝶/解决方案的适应度相关。这意味着,当一只蝴蝶散发出更多的香味时,周围的其他蝴蝶可以感知到并被吸引。幂是强度增加的指数。参数a允许正则表达式、线性响应和响应压缩。响应扩展是当I增加时,香味(f)比I增长更快。响应压缩是当I增加时,f比I增长慢。线性响应是当I增加时,f成比例地增加。经实验证明,有时随着刺激的增强,昆虫对刺激变化的敏感性变得越来越低。因此在BOA中,为了估计I的大小,使用了响应压缩。
蝴蝶的自然现象基于两个重要问题:I的变化和f的表示。简单地说,蝴蝶的I与编码后的目标函数相关联。但是,f是相对的,即应该由其他蝴蝶来感知。史蒂文斯幂定律中,为了将气味与其他形式区别开来,使用了C。现在,当I较少的蝴蝶向I较多的蝴蝶移动时,f比I增加得更快。因此,我们应该允许f随幂指数参数a实现的吸收程度而变化。在BOA中,香味被表示为刺激物的物理强度的函数,如下所示:

在这里插入图片描述

3 具体算法

为了用搜索算法演示上述讨论,将蝴蝶的上述特征理想化如下:

(1)所有的蝴蝶都可以发出气味,这使蝴蝶间相互吸引。

(2)每只蝴蝶都会随机移动或朝最好的蝴蝶移动,散发出更多的芳香。

(3)蝴蝶的刺激强度受目标函数的景观影响或决定。

该算法分为三个阶段:(1)初始化阶段、(2)迭代阶段和(3)结束阶段。

在BOA的每次运行中,首先执行初始化阶段,然后进行迭代搜索,最后在找到最优解时终止算法。BOA中使用的参数值也会被分配,设置这些值后,算法将继续创建初始蝴蝶种群以进行优化。由于在BOA的模拟过程中蝴蝶总数保持不变,分配了一个固定大小的内存来存储信息。蝴蝶的位置是在搜索空间中随机生成的,并计算和存储它们的香味和适应值。这样就完成了初始化阶段,算法开始了迭代阶段,该阶段使用创建的人工蝶形执行搜索。算法的第二阶段,即迭代阶段,由算法执行多次迭代。在每次迭代中,解空间中的所有蝶形都移到新位置,然后重新评估其适应性值。算法首先计算解空间中不同位置的所有蝴蝶的适应度值。那么这些蝴蝶就会利用式1在自己的位置产生香味。该算法有两个关键步骤,即全局搜索阶段和局部搜索阶段。在全局搜索阶段,蝴蝶向最合适的蝴蝶/解g?迈出一步,该蝴蝶/解g可以用公式(2)来表示。

在这里插入图片描述

这里,g?表示在当前迭代的所有解中找到的当前最佳解;fi表示第i只蝴蝶的香味,r是[0,1]中的随机数。局部搜索阶段可以表示为

在这里插入图片描述

其中,xjt和xkt是解空间中的第j个蝴蝶和第k个蝴蝶。
蝴蝶寻找食物、交配伙伴可以在局部和全局范围内发生。考虑到地理上的接近和各种其他因素,如雨、风等,在整个交配伙伴或蝴蝶的觅食活动中,寻找食物可能占很大比例。因此,在BOA中使用切换概率p来在普通全局搜索和密集局部搜索之间切换。
在未达到停止标准之前,一直进行迭代。迭代结束的标准可以有多个,如使用的最大CPU时间、达到的最大迭代次数、没有改进的最大迭代次数、达到错误率的特定值或任何其他适当的标准。当迭代阶段结束时,算法输出具有最佳适应度的最优解。

三、部分源代码

clc;clear
%% 下载数据
% 加工数据包括加工时间,加工机器,机器数,各机器权重,工件数,各工件对应的工序数
load data operation_time operation_machine num_machine machine_weight num_job num_op
%% 基本参数
MAXGEN=200;             % 最大迭代次数
sizepop=201;            % 种群规模
e=0.5;                  % 目标值权重
N_size=30;              % 邻域解数量
S_size=15;              % 共享解数量
G=5;                    % 巡回次数
G1=20;                  % 竞争机制1参数
G2=10;                  % 竞争机制2参数
trace=zeros(2,MAXGEN);
chrom_best=[];

%% ===========================种群初始化============================
total_op_num=sum(num_op);
chroms=initialization(num_op,num_job,total_op_num,sizepop,operation_machine,operation_time);
[Z,~,~,~,~]=fitness(chroms,num_machine,e,num_job,num_op);
% 将最好的解划分为领飞鸟
[Z_leader,ind]=min(Z);
leader=chroms(ind,:); 
% 从chroms中移出领飞鸟,然后划分左右两个跟飞鸟种群
chroms(ind,:)=[];
Z(ind)=[];
sp=(sizepop-1)/2;
lefts=chroms(1:sp,:);
Z_left=Z(1:sp);
rights=chroms(sp+1:end,:);
Z_right=Z(sp+1:end);
%% 候鸟算法中的交叉函数与遗传算法的不同
%% 候鸟算法输入两个染色体种群,分别来自左右队列
%--------------------------------------------------------------------------
function [lefts,Z_left,rights,Z_right]= crossover(lefts,rights,Z_left,Z_right,total_op_num,num_machine,e,num_job,num_op)
chroms1=lefts;
chroms2=rights;
for i=1:size(chroms1,1)
    %% 面向工序码的交叉操作
    % 父代染色体
    parent1=lefts(i,:);
    parent2=rights(i,:);
    Job=randperm(num_job);
    % 将工件随机分成两个集合
    J1=Job(1:round(num_job/2));
    J2=Job(length(J1)+1:end);
    op_p1=[];
    op_p2=[];
    for j=1:length(J2)
        %找出父代中J2片段对应的位置
        op_p1=[op_p1,find(parent1(1:total_op_num)==J2(j))];
        op_p2=[op_p2,find(parent2(1:total_op_num)==J2(j))];
    end
    op_s1=sort(op_p1);
    op_s2=sort(op_p2);
    % 子代1交换J2片段的基因,机器码对应位置的基因,工时码对应位置的基因
    chroms1(i,op_s1)=parent2(op_s2);
    chroms1(i,total_op_num+op_s1)=parent2(total_op_num+op_s2);
    chroms1(i,total_op_num*2+op_s1)=parent2(total_op_num*2+op_s2);
    % 子代2同理
    chroms2(i,op_s2)=parent1(op_s1);
    chroms2(i,total_op_num+op_s2)=parent1(total_op_num+op_s1);
    chroms2(i,total_op_num*2+op_s2)=parent1(total_op_num*2+op_s1);
    %% 面向机器码的交叉操作
    parent1=chroms1(i,:);
    parent2=chroms2(i,:);
    % 随机产生与染色体长度相等的0,1序列
    rand0_1=randi([0,1],1,total_op_num);
    for n=1:num_job
        ind_0=find(rand0_1(num_op(n)*(n-1)+1:num_op(n)*n)==0);
        if ~isempty(ind_0)
            ind1=find(parent1(1:total_op_num)==n);
            ind2=find(parent2(1:total_op_num)==n);
            chroms1(i,total_op_num+ind1(ind_0))=parent2(total_op_num+ind2(ind_0));
            chroms1(i,total_op_num*2+ind1(ind_0))=parent2(total_op_num*2+ind2(ind_0));
            chroms2(i,total_op_num+ind2(ind_0))=parent1(total_op_num+ind1(ind_0));
            chroms2(i,total_op_num*2+ind2(ind_0))=parent1(total_op_num*2+ind1(ind_0));
        end
    end
end
%% 判断个体是否可以更新
[Z1,~,~,~,~]=fitness(chroms1,num_machine,e,num_job,num_op);
[Z2,~,~,~,~]=fitness(chroms2,num_machine,e,num_job,num_op);
lefts(Z1<Z_left,:)=chroms1(Z1<Z_left,:);
Z_left(Z1<Z_left)=Z1(Z1<Z_left);
rights(Z2<Z_right,:)=chroms2(Z2<Z_right,:);
Z_right(Z2<Z_right)=Z2(Z2<Z_right);
function [Z,makespan,machine_load,machine_weight,pvals] = fitness(chroms,num_machine,e,num_job,num_op)
sizepop=size(chroms,1);
pvals=cell(1,sizepop);
makespan=zeros(1,sizepop);
machine_load=makespan;
total_op_num=sum(num_op);  % 总工序数
for k=1:sizepop
    chrom=chroms(k,:);
    machine=zeros(1,num_machine);  % 记录各机器变化时间
    job=zeros(1,num_job);  % 记录各工件变化时间
    machine_time=zeros(1,num_machine);  % 计算各机器的实际加工时间
    pval=zeros(2,total_op_num);  % 记录各工序开始和结束时间
    for i=1:total_op_num
        else
            pval(1,i)=job(chrom(i));
            job(chrom(i))=job(chrom(i))+chrom(total_op_num*2+i);
            machine(chrom(total_op_num+i))=job(chrom(i));
            pval(2,i)=job(chrom(i));
        end
        machine_time(chrom(total_op_num+i))=machine_time(chrom(total_op_num+i))+chrom(total_op_num*2+i);
    end
    makespan(k)=max(machine);
    machine_weight=machine_time;
    machine_load(k)=max(machine_weight)-min(machine_weight);
    pvals{k}=pval;
end
Z=e*makespan+(1-e)*machine_load;

四、运行结果

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例(第2版)[M].电子工业出版社,2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社,2017.
[3]蝴蝶优化算法

您可能感兴趣的文章:

相关文章