如何使用matlab鳥群算法求解車間調度問題

這篇文章主要介紹了如何使用matlab鳥群算法求解車間調度問題，具有一定借鑒價值，感興趣的朋友可以參考下，希望大家閱讀完這篇文章之后大有收獲，下面讓小編帶著大家一起了解一下。

一、車間調度簡介

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鳥群算法求解車間調度問題”這篇文章對大家有幫助，同時也希望大家多多支持億速云，關注億速云行業資訊頻道，更多相關知識等著你來學習!