matlab鳥群算法求解車間調(diào)度問題詳解及實現(xiàn)源碼

更新時間：2022年02月07日 14:32:18 作者：紫極神光

這篇文章主要為大家介紹了matlab鳥群算法求解車間調(diào)度的問題分析及實現(xiàn)源碼，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進步

一、車間調(diào)度簡介

1 車間調(diào)度定義

車間調(diào)度是指根據(jù)產(chǎn)品制造的合理需求分配加工車間順序，從而達到合理利用產(chǎn)品制造資源、提高企業(yè)經(jīng)濟效益的目的。車間調(diào)度問題從數(shù)學上可以描述為有n個待加工的零件要在m臺機器上加工。問題需要滿足的條件包括每個零件的各道工序使用每臺機器不多于1次，每個零件都按照一定的順序進行加工。

2 傳統(tǒng)作業(yè)車間調(diào)度

傳統(tǒng)作業(yè)車間帶調(diào)度實例

在這里插入圖片描述

有若干工件，每個工件有若干工序，有多個加工機器，但是每道工序只能在一臺機器上加工。對應到上面表格中的實例就是，兩個工件，工件J1有三道工序，工序Q11只能在M3上加工，加工時間是5小時。
約束是對于一個工件來說，工序的相對順序不能變。O11->O12->O13。每時刻，每個工件只能在一臺機器上加工；每個機器上只能有一個工件。
調(diào)度的任務則是安排出工序的加工順序，加工順序確定了，因為每道工序只有一臺機器可用，加工的機器也就確定了。
調(diào)度的目的是總的完工時間最短（也可以是其他目標）。舉個例子，比如確定了O21->O22->O11->O23->O12->O13的加工順序之后，我們就可以根據(jù)加工機器的約束，計算出總的加工時間。
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 柔性作業(yè)車間調(diào)度

柔性作業(yè)車間帶調(diào)度實例（參考自高亮老師論文
《改進遺傳算法求解柔性作業(yè)車間調(diào)度問題》——機械工程學報）

在這里插入圖片描述

相比于傳統(tǒng)作業(yè)車間調(diào)度，柔性作業(yè)車間調(diào)度放寬了對加工機器的約束，更符合現(xiàn)實生產(chǎn)情況，每個工序可選加工機器變成了多個，可以由多個加工機器中的一個加工。比如上表中的實例，J1的O12工序可以選擇M2和M4加工，加工時間分別是8小時和4小時，但是并不一定選擇M4加工，最后得出來的總的完工時間就更短，所以，需要調(diào)度算法求解優(yōu)化。

相比于傳統(tǒng)作業(yè)車間，柔性車間作業(yè)調(diào)度的調(diào)度任務不僅要確定工序的加工順序，而且需要確定每道工序的機器分配。比如，確定了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]的機器組合。調(diào)度的目的還是總的完工時間最短（也可以是其他目標，比如機器最大負荷最短、總的機器負荷最短）

二、蝴蝶優(yōu)化算法（MBO）簡介

1 介紹

蝴蝶優(yōu)化算法(butterfly optimization algorithm, BOA)是Arora 等人于2019年提出的一種元啟發(fā)式智能算法。該算法受到了蝴蝶覓食和交配行為的啟發(fā)，蝴蝶接收/感知并分析空氣中的氣味，以確定食物來源/交配伙伴的潛在方向。
蝴蝶利用它們的嗅覺、視覺、味覺、觸覺和聽覺來尋找食物和伴侶，這些感覺也有助于它們從一個地方遷徙到另一個地方，逃離捕食者并在合適的地方產(chǎn)卵。在所有感覺中，嗅覺是最重要的，它幫助蝴蝶尋找食物（通常是花蜜）。蝴蝶的嗅覺感受器分散在蝴蝶的身體部位，如觸角、腿、觸須等。這些感受器實際上是蝴蝶體表的神經(jīng)細胞，被稱為化學感受器。它引導蝴蝶尋找最佳的交配對象，以延續(xù)強大的遺傳基因。雄性蝴蝶能夠通過信息素識別雌性蝴蝶，信息素是雌性蝴蝶發(fā)出的氣味分泌物，會引起特定的反應。
通過觀察，發(fā)現(xiàn)蝴蝶對這些來源的位置有非常準確的判斷。此外，它們可以辨識出不同的香味，并感知它們的強度。蝴蝶會產(chǎn)生與其適應度相關(guān)的某種強度的香味，即當蝴蝶從一個位置移動到另一個位置時，它的適應度會相應地變化。當蝴蝶感覺到另一只蝴蝶在這個區(qū)域散發(fā)出更多的香味時，就會去靠近，這個階段被稱為全局搜索。另外一種情況，當蝴蝶不能感知大于它自己的香味時，它會隨機移動，這個階段稱為局部搜索。

2 香味

為了理解BOA中的香味是如何計算的，首先需要理解，像氣味、聲音、光、溫度等這樣的模態(tài)是如何計算的。感知、處理這些模態(tài)需要知道三個重要的術(shù)語：感覺模態(tài)C、刺激強度I和冪指數(shù)a。在感覺模態(tài)中，感覺意味著測量能量的形式并以類似方式對其進行處理，而模態(tài)是指傳感器使用的原始輸入。不同的形態(tài)可以是氣味，聲音，光線，溫度，在BOA中，模態(tài)是香味。I是物理刺激的大小。在BOA中，I與蝴蝶/解決方案的適應度相關(guān)。這意味著，當一只蝴蝶散發(fā)出更多的香味時，周圍的其他蝴蝶可以感知到并被吸引。冪是強度增加的指數(shù)。參數(shù)a允許正則表達式、線性響應和響應壓縮。響應擴展是當I增加時，香味（f）比I增長更快。響應壓縮是當I增加時，f比I增長慢。線性響應是當I增加時，f成比例地增加。經(jīng)實驗證明，有時隨著刺激的增強，昆蟲對刺激變化的敏感性變得越來越低。因此在BOA中，為了估計I的大小，使用了響應壓縮。
蝴蝶的自然現(xiàn)象基于兩個重要問題：I的變化和f的表示。簡單地說，蝴蝶的I與編碼后的目標函數(shù)相關(guān)聯(lián)。但是，f是相對的，即應該由其他蝴蝶來感知。史蒂文斯冪定律中，為了將氣味與其他形式區(qū)別開來，使用了C?，F(xiàn)在，當I較少的蝴蝶向I較多的蝴蝶移動時，f比I增加得更快。因此，我們應該允許f隨冪指數(shù)參數(shù)a實現(xiàn)的吸收程度而變化。在BOA中，香味被表示為刺激物的物理強度的函數(shù)，如下所示：

在這里插入圖片描述

3 具體算法

為了用搜索算法演示上述討論，將蝴蝶的上述特征理想化如下：

（1）所有的蝴蝶都可以發(fā)出氣味，這使蝴蝶間相互吸引。

（2）每只蝴蝶都會隨機移動或朝最好的蝴蝶移動，散發(fā)出更多的芳香。

（3）蝴蝶的刺激強度受目標函數(shù)的景觀影響或決定。

該算法分為三個階段：(1)初始化階段、(2)迭代階段和(3)結(jié)束階段。

在BOA的每次運行中，首先執(zhí)行初始化階段，然后進行迭代搜索，最后在找到最優(yōu)解時終止算法。BOA中使用的參數(shù)值也會被分配，設(shè)置這些值后，算法將繼續(xù)創(chuàng)建初始蝴蝶種群以進行優(yōu)化。由于在BOA的模擬過程中蝴蝶總數(shù)保持不變，分配了一個固定大小的內(nèi)存來存儲信息。蝴蝶的位置是在搜索空間中隨機生成的，并計算和存儲它們的香味和適應值。這樣就完成了初始化階段，算法開始了迭代階段，該階段使用創(chuàng)建的人工蝶形執(zhí)行搜索。算法的第二階段，即迭代階段，由算法執(zhí)行多次迭代。在每次迭代中，解空間中的所有蝶形都移到新位置，然后重新評估其適應性值。算法首先計算解空間中不同位置的所有蝴蝶的適應度值。那么這些蝴蝶就會利用式1在自己的位置產(chǎn)生香味。該算法有兩個關(guān)鍵步驟，即全局搜索階段和局部搜索階段。在全局搜索階段，蝴蝶向最合適的蝴蝶/解g∗邁出一步，該蝴蝶/解g可以用公式(2)來表示。

在這里插入圖片描述

這里，g∗表示在當前迭代的所有解中找到的當前最佳解；fi表示第i只蝴蝶的香味，r是[0，1]中的隨機數(shù)。局部搜索階段可以表示為

在這里插入圖片描述

其中，xjt和xkt是解空間中的第j個蝴蝶和第k個蝴蝶。
蝴蝶尋找食物、交配伙伴可以在局部和全局范圍內(nèi)發(fā)生?？紤]到地理上的接近和各種其他因素，如雨、風等，在整個交配伙伴或蝴蝶的覓食活動中，尋找食物可能占很大比例。因此，在BOA中使用切換概率p來在普通全局搜索和密集局部搜索之間切換。
在未達到停止標準之前，一直進行迭代。迭代結(jié)束的標準可以有多個，如使用的最大CPU時間、達到的最大迭代次數(shù)、沒有改進的最大迭代次數(shù)、達到錯誤率的特定值或任何其他適當?shù)臉藴?。當?shù)A段結(jié)束時，算法輸出具有最佳適應度的最優(yōu)解。

三、部分源代碼

clc;clear
%% 下載數(shù)據(jù)
% 加工數(shù)據(jù)包括加工時間，加工機器，機器數(shù)，各機器權(quán)重，工件數(shù)，各工件對應的工序數(shù)
load data operation_time operation_machine num_machine machine_weight num_job num_op
%% 基本參數(shù)
MAXGEN=200;             % 最大迭代次數(shù)
sizepop=201;            % 種群規(guī)模
e=0.5;                  % 目標值權(quán)重
N_size=30;              % 鄰域解數(shù)量
S_size=15;              % 共享解數(shù)量
G=5;                    % 巡回次數(shù)
G1=20;                  % 競爭機制1參數(shù)
G2=10;                  % 競爭機制2參數(shù)
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);
% 將最好的解劃分為領(lǐng)飛鳥
[Z_leader,ind]=min(Z);
leader=chroms(ind,:); 
% 從chroms中移出領(lǐng)飛鳥，然后劃分左右兩個跟飛鳥種群
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);
%% 候鳥算法中的交叉函數(shù)與遺傳算法的不同
%% 候鳥算法輸入兩個染色體種群，分別來自左右隊列
%--------------------------------------------------------------------------
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,:);
    % 隨機產(chǎn)生與染色體長度相等的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);  % 總工序數(shù)
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);  % 記錄各工序開始和結(jié)束時間
    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;

四、運行結(jié)果

在這里插入圖片描述