一、簡介

Grey Wolf Optimizer是Seyedali Mirjalili受大灰狼捕食策略的啟發(fā)，于2014年提出的一種元啟發(fā)式算法，主要模擬了搜索獵物、包圍獵物和攻擊獵物，源代碼關(guān)注公眾號后，回復(fù)"灰狼"或"GWO"獲取。
啟發(fā)
灰狼屬于犬科動物，是食物鏈頂端的頂級掠食者，它們大多喜歡群居生活，每個種群平均5~12不等。特別有趣的是，它們有非常嚴(yán)格的社會等級，如下圖所示。

圖1 灰狼分層

領(lǐng)頭狼（領(lǐng)導(dǎo)者）是一公一母，稱為alphas（個人理解，之所以這么叫，是因為alpha是希臘字幕中的第一個，用來表示最前面的）。alpha狼主要負(fù)責(zé)決策狩獵、睡覺地點、起床時間等等，然后再將決策下達(dá)至整個種群，如果狼群垂下尾巴說明它們都認(rèn)可。然而種群中也會存在民主行為，alpha狼也會聽從種群中的其他狼。有趣的是，alpha狼不一定是最強(qiáng)壯的成員，但在管理團(tuán)隊方面是最好的。這說明一個狼群的組織和紀(jì)律遠(yuǎn)比它的力量重要。

第二層是beta。deta狼是從屬狼，可公可母，用于輔助alpha狼制定決策或其他種群活動。當(dāng)其中一只alpha狼去世或變老時，它就是alpha狼的最佳替補。beta狼應(yīng)該遵從alpha狼，但也會命令其他低級別的狼，因而起到承上啟下的作用。
最底層的是omega狼。它們扮演了"替罪羊"的角色（好慘）,必須屈服于其他領(lǐng)頭狼，進(jìn)食時也是排在最后?？雌饋韔mega狼在狼群中并不是一個重要的個體，但是一旦失去omega狼，整個狼群就會面臨內(nèi)部爭斗和問題。
如果一頭狼既不是alpha，beta，也不是omega，就稱為從屬狼或delta狼。delta狼必須聽從于alpha和beta狼，但會支配omega狼。偵察狼、守衛(wèi)狼、老狼、捕食狼和看管狼都是這一類。偵察狼負(fù)責(zé)監(jiān)視領(lǐng)地的邊界，一旦有危險就向狼群發(fā)出警告；守衛(wèi)狼保護(hù)和保證狼群的安全；老狼是從alpha或beta退下來的經(jīng)驗豐富的狼；捕食狼幫助alpha和beta捕獵并為狼群提供食物；看管狼負(fù)責(zé)照顧狼群中的老弱病殘。
除此之外，群體狩獵是灰狼的另一個有趣的社會行為，主要包括以下幾個階段：
跟蹤、追逐和接近獵物；
追捕、包圍和騷擾獵物，直到它停止移動；
向獵物發(fā)起攻擊。

圖2 灰狼捕食行為：(A)追逐、接近和跟蹤獵物。(B-D)追捕、騷擾和包圍獵物。(E)靜止?fàn)顟B(tài)與攻擊
數(shù)學(xué)模型和算法
首先給出社會等級、跟蹤、包圍和攻擊獵物的數(shù)學(xué)模型，然后再提供完整的GWO算法。
社會等級
為了在設(shè)計GWO時對狼的社會等級進(jìn)行數(shù)學(xué)建模，認(rèn)為最合適的解是alpha(α \alphaα)，那么第二和第三最優(yōu)解分別表示為beta(β \betaβ)和delta(δ \deltaδ),而剩余其他解都假定為omega(ω \omegaω)。在GWO中，通過α \alphaα、β \betaβ和δ \deltaδ來導(dǎo)引捕食（優(yōu)化），ω \omegaω聽從于這三種狼。
包圍捕食
灰狼在捕食時會將獵物包圍，使用下式進(jìn)行表達(dá)這種行為：
D ? = ∣ C ? ? X p → ( t ) ? X ? ( t ) ∣ (1) \vec{D}=|\vec{C} \cdot \overrightarrow{X_{p}}(t)-\vec{X}(t)|\tag{1}D=∣C?Xp(t)?X(t)∣(1)
X ? ( t + 1 ) = X p → ( t ) ? A ? ? D ? (2) \vec{X}(t+1)=\overrightarrow{X_{p}}(t)-\vec{A} \cdot \vec{D}\tag{2}X(t+1)=Xp(t)?A?D(2)
其中t tt表示當(dāng)前迭代次數(shù)，A ? \vec{A}A和C ? \vec{C}C為系數(shù)向量，X p → \overrightarrow{X_{p}}Xp是獵物的位置向量，X ? \vec{X}X是灰狼的位置向量。
向量A ? \vec{A}A和C ? \vec{C}C的計算如下：
A ? = 2 a ? ? r ? 1 ? a ? (3) \vec{A}=2 \vec{a} \cdot \vec{r}{1}-\vec{a}\tag{3}A=2a?r1?a(3)
C ? = 2 ? r 2 → (4) \vec{C}=2 \cdot \overrightarrow{r{2}}\tag{4}C=2?r2(4)
其中a ? \vec{a}a的各個分量在迭代過程中線性地從2減少到0，r ? 1 \vec{r}{1}r1和r 2 → \overrightarrow{r{2}}r2為[0,1]之間的隨機(jī)向量。
為了清楚地反映等式(1)和(2)的效果，圖3(a)中顯示了二維的位置向量以及可能的鄰域，可以看出，灰狼的位置( X , Y ) (X,Y)(X,Y)可以根據(jù)獵物的位置( X ? , Y ? ) (X*,Y)(X?,Y?)進(jìn)行更新，通過調(diào)整A ? \vec{A}A和C ? \vec{C}C的值，可以在最優(yōu)代理周圍到達(dá)相對于當(dāng)前位置的不同地方。例如，當(dāng)A ? = ( 0 , 1 ) \vec{A}=(0,1)A=(0,1)和C ? = ( 1 , 1 ) \vec{C}=(1,1)C=(1,1)時，灰狼的新位置為( X ? ? X , Y ? ) (X*-X,Y)(X??X,Y?)。三維空間中也是類似。注意，此處的兩張圖僅僅展示了A ? = ( 0 , 1 ) \vec{A}=(0,1)A=(0,1)和C ? = ( 1 , 1 ) \vec{C}=(1,1)C=(1,1)這一種情況，當(dāng)隨機(jī)向量r 1 r_1r1和r 2 r_2r2取不同的值時，灰狼可以到達(dá)任意兩點之間的位置。同時還注意到，A AA的取值不同還會決定灰狼靠近還是遠(yuǎn)離獵物，后面再詳細(xì)說明。

圖3 2D和3D位置向量及其可能的下一位置
狩獵
灰狼能夠識別獵物的位置并包圍它們，狩獵通常是是由alpha狼領(lǐng)導(dǎo)的，beta和和delta狼偶爾也會參與狩獵。然而，在一個抽象的搜索空間中，我們不知道最佳(獵物)的位置。為了在數(shù)學(xué)上模擬灰狼的狩獵行為，我們假設(shè)alpha(最佳候選解)、beta和delta狼對獵物的潛在位置有更好的了解。因此，我們保存到目前為止獲得的前三個最佳解，并迫使其他搜索代理(包括omegas)根據(jù)最佳搜索代理的位置更新其位置。對此，提出以下公式：
D α → = ∣ C ? 1 ? X α → ? X ? ∣ , D β → = ∣ C ? 2 ? X β → ? X ? ∣ , D δ → = ∣ C 3 → ? X δ → ? X ? ∣ (5) \overrightarrow{D_{\alpha}}=\left|\vec{C}{1} \cdot \overrightarrow{X{\alpha}}-\vec{X}\right|, \overrightarrow{D_{\beta}}=\left|\vec{C}{2} \cdot \overrightarrow{X{\beta}}-\vec{X}\right|, \overrightarrow{D_{\delta}}=|\overrightarrow{C_{3}} \cdot \overrightarrow{X_{\delta}}-\vec{X}|\tag{5}Dα=∣∣∣C1?Xα?X∣∣∣,Dβ=∣∣∣C2?Xβ?X∣∣∣,Dδ=∣C3?Xδ?X∣(5)

X 1 → = X α → ? A 1 → ? ( D α → ) , X 2 → = X β → ? A 2 → ? ( D β → ) , X 3 → = X δ → ? A 3 → ? ( D δ → ) 6 ) (() \overrightarrow{X_{1}}=\overrightarrow{X_{\alpha}}-\overrightarrow{A_{1}} \cdot(\overrightarrow{D_{\alpha}}), \overrightarrow{X_{2}}=\overrightarrow{X_{\beta}}-\overrightarrow{A_{2}} \cdot(\overrightarrow{D_{\beta}}), \overrightarrow{X_{3}}=\overrightarrow{X_{\delta}}-\overrightarrow{A_{3}} \cdot(\overrightarrow{D_{\delta}})\tag(6)X1=Xα?A1?(Dα),X2=Xβ?A2?(Dβ),X3=Xδ?A3?(Dδ)6)(()

X ? ( t + 1 ) = X 1 → + X 2 → + X 3 → 3 (7) \vec{X}(t+1)=\frac{\overrightarrow{X_{1}}+\overrightarrow{X_{2}}+\overrightarrow{X_{3}}}{3}\tag{7}X(t+1)=3X1+X2+X3(7)

圖4展示了2D空間中如何根據(jù)alpha，beta和delta進(jìn)行代理位置的更新?？梢钥吹?，最終位置將是圓內(nèi)的一個隨機(jī)位置，該圓由alpha、beta和delta定義，換句話說，alpha、beta和delta狼對獵物的位置進(jìn)行估計，而其他狼則再獵物周圍隨機(jī)更新它們的位置。

圖4 GWO中的位置更新示意圖

攻擊獵物（利用）
正如上面提到的，當(dāng)獵物停止移動時，灰狼就會攻擊它來完成狩獵。為了對接近獵物進(jìn)行建模，需要不斷降低a ? \vec aa的值，那么A ? \vec AA的波動范圍也會降低。當(dāng)A ? ∈ [ ? 1 , 1 ] \vec A\in[-1,1]A∈[?1,1]，搜索代理的下一位置可以是代理當(dāng)前位置和獵物位置之間的任意位置。圖5(a)表明當(dāng)∣ A ∣ < 1 |A|<1∣A∣<1時，灰狼向獵物發(fā)起攻擊。

圖5 獵物攻擊vs搜索獵物

搜索獵物（探索）
灰狼通常根據(jù)alpha、beta和delta狼的位置進(jìn)行搜索，它們彼此分散尋找獵物，然后匯聚攻擊獵物。為了數(shù)學(xué)上對分散建模，利用隨機(jī)值大于1或小于-1的A ? \vec AA迫使搜索代理偏離獵物，從而保證了探索。圖5(b)表明當(dāng)∣ A ∣ > 1 |A|>1∣A∣>1時，迫使灰狼離開獵物，希望能找到更合適的獵物。另一個支持GWO進(jìn)行探索的因素是C ? \vec CC，根據(jù)公式(4)，C ? \vec CC的取值范圍為[ 0 , 2 ] [0,2][0,2]，該分量為獵物提供隨機(jī)權(quán)重，以隨機(jī)強(qiáng)調(diào)(C>1)或弱化(C<1)
獵物在定義等式(1)中的距離時的作用。這有助于GWO在整個優(yōu)化過程中表現(xiàn)出更隨機(jī)的行為，有利于探索和避免局部最優(yōu)。C CC并不是和A AA一樣線性遞減，特意要求C CC在任何時候都提供隨機(jī)值，以便不僅在初始迭代中強(qiáng)調(diào)探索，而且在最終迭代中也強(qiáng)調(diào)探索。

C CC向量也可以被認(rèn)為是在自然界中障礙物對接近獵物的影響，一般來說，自然中的障礙出現(xiàn)在狼的捕獵路徑上，實際上阻礙了它們快速、方便地接近獵物。這就是向量C CC的作用。根據(jù)狼所處的位置，它可以隨機(jī)地給獵物一個權(quán)重，從而讓狼的捕食變得更加困難和遙遠(yuǎn)，反之亦然。

GWO算法
總之，搜索過程從在GWO算法中創(chuàng)建一個隨機(jī)的灰狼種群(候選解)開始。在迭代過程中，alpha、beta和delta狼估計獵物可能的位置。每一個候選解更新它與獵物的距離。為了分別強(qiáng)調(diào)探索和利用，將參數(shù)a aa從2降低到0。當(dāng)∣ A ? ∣ > 1 |\vec A|>1∣A∣>1時,候選解有偏離獵物的傾向，當(dāng)∣ A ? ∣ < 1 |\vec A|<1∣A∣<1時,候選解收斂于獵物。最后，當(dāng)滿足結(jié)束條件時終止GWO算法。GWO的偽代碼如下。

初始化灰狼種群X i ( i = 1 , 2 , . . . , n ) X_i(i=1,2,…,n)Xi(i=1,2,…,n)

初始化a , A , C a,A,Ca,A,C

計算每個搜索代理的適應(yīng)度值

X α = X_{\alpha}=Xα=最優(yōu)搜索代理

X β = X_{\beta}=Xβ=第二優(yōu)搜索代理

X δ = X_{\delta}=Xδ=第三優(yōu)搜索代理

while(t<最大迭代次數(shù))

for 每個搜索代理

根據(jù)等式(7)更新當(dāng)前代理的位置

end for

更新a , A , C a,A,Ca,A,C

計算所有搜索代理的適應(yīng)度值

更新X α X_{\alpha}Xα、X β X_{\beta}Xβ、X δ X_{\delta}Xδ

KaTeX parse error: Expected ‘EOF’, got ‘&’ at position 1: &?emsp;?t=t+…

end while

return X α X_{\alpha}Xα

通過以下幾點可以了解GWO在理論上是如何解決優(yōu)化問題的：

所提出的社會等級有助于GWO在迭代過程中保存目前的最優(yōu)解；
所提出的包圍機(jī)制在解周圍定義了一個圓形鄰域，該鄰域可以作為超球體擴(kuò)展到更高維度；
隨機(jī)參數(shù)A AA和C CC輔助候選解具有不同隨機(jī)半徑的超球體；
所提出的狩獵方法允許候選解確定獵物的可能位置；
a aa和A AA的適配保證了探索和利用；
參數(shù)a aa和A AA的自適應(yīng)值使GWO在探索和利用之間實現(xiàn)平穩(wěn)過渡；
隨著A AA的下降，一半迭代致力于探索(∣ A ∣ ≥ 1 |A|≥1∣A∣≥1)，另一半致力于利用(∣ A ∣ < 1 |A|<1∣A∣<1)；
GWO只有兩個主要參數(shù)需要調(diào)整(a aa和C CC)。

二、源代碼

%該程序用于解決柔性作業(yè)車間調(diào)度，m個工件，n道工序，其中n為最大工序數(shù)，工件的工序 %數(shù)可以少于n，加工機(jī)器數(shù)為M，每個工件的每道工序具有多個機(jī)器可以選擇，對應(yīng)的時間 %不同，其中初始種群的儲存方式采用cell數(shù)據(jù)類型 %Version:1.3 %fileDescription:調(diào)度機(jī)器可選的柔性作業(yè)車間問題，甘特圖已完善,GWO,8*8實例 %last edit time:2019-6-7 function GWO_Model_FJSP_1_3_8_8() count = 600; ? ? %迭代次數(shù) N = 50; ? ? ? ? ?%種群規(guī)模 m = 8; ? ? ? ? ? ? %工件數(shù) n = 4; ? ? ? ? ? ? %工序數(shù) M = 8; ? ? ? ? ? ? %機(jī)器數(shù) a =2; ? ? ? ? ? ? ?%計算A/C協(xié)同系數(shù)的 plotif = 1; ? ? ? ?%控制程序是否進(jìn)行繪圖 s = input(m,n); ? ?%數(shù)據(jù)輸入 [p,TN] = initial_p(m,n,N,s,M); ? ?%生成初始種群50,采用細(xì)胞結(jié)構(gòu)，每個元素為8*4 P = machine(n,M); FIT = zeros(count,1); aveFIT = zeros(count,1); X1=randperm(count); ? ? ? %收斂圖形的橫坐標(biāo)X X=sort(X1); %------------------------輸出最優(yōu)解的時有用------------------------------ best_fit = 1000; ? ? ? ? ? ?%改變模型需要修改此參數(shù) best_p = zeros(m,n); best_TN = zeros(m,n); Y1p = zeros(m,1); Y2p = zeros(m,1); Y3p = zeros(m,1); minfit3 ?= ?1000000000; %-------------------------進(jìn)行迭代-------------------------------------- for i = 1:count ? ?[fit,Y1,Y2,Y3] = object(p,TN,N,P,m,n); ? ? ?[newp,newTN] = GWO(fit,p,TN,N,m,n,s,a); ? ?a = a-2/(count-1); ? ? ? ?%a的值會線性下降 ? ?if best_fit > min(fit) ? ? ? ?[best_p,best_TN,best_fit,Y1p,Y2p,Y3p]=best(best_fit,best_p,fit,best_TN,Y1p,Y2p,Y3p,p,TN,Y1,Y2,Y3); ? ?end ? ?p = newp; ? ?TN = newTN; ? ?minfit = min(fit); ? ?if minfit3>minfit ? ? ? ?minfit3 = minfit; ? ?end ? ?FIT(i) = minfit3; ? ?%用于適應(yīng)度函數(shù)的 ? ?aveFIT(i) = mean(fit); ? ? ?%用于適應(yīng)度函數(shù)的 end %------------------投射最佳方案數(shù)據(jù)-------------------------------------- ? ? ?fprintf('最優(yōu)解：%d\n',best_fit); ? ?fprintf('工序1 工序2 工序3 工序4\n'); ? ?best_p ? ?fprintf('時間1 時間2 時間3 時間4\n'); ? ?best_TN %------------------------收斂曲線---------------------------------------- ? ?if plotif == 1 ? ?figure; ? ?plot(X,FIT,'r'); ? ?hold on; ? ?plot(X,aveFIT,'b'); ? ?title('收斂曲線'); ? ?hold on; ? ?legend('最優(yōu)解','平均值'); %-------------------------甘特圖----------------------------------------- figure; w=0.5; ? ? ? %橫條寬度 set(gcf,'color','w'); ? ? ?%圖的背景設(shè)為白色 for i = 1:m ? ?for j = 1:n ? ? ? ?color=[1,0.98,0.98;1,0.89,0.71;0.86,0.86,0.86;0.38,0.72,1;1,0,1;0,1,1;0,1,0.49;1,0.87,0.67;0.39,0.58,0.92;0.56,0.73,0.56]; ? ? ? ?a = [Y1p(i,j),Y2p(i,j)]; ? ? ? ?x=a(1,[1 1 2 2]); ? ? ?%設(shè)置小圖框四個點的x坐標(biāo) ? ? ? ?y=Y3p(i,j)+[-w/2 w/2 w/2 -w/2]; ? %設(shè)置小圖框四個點的y坐標(biāo) ? ? ? ?color = [color(i,1),color(i,2),color(i,3)]; ? ? ? ?p=patch('xdata',x,'ydata',y,'facecolor',color,'edgecolor','k'); ? ?%facecolor為填充顏色，edgecolor為圖框顏色 ? ? ? ? ? ?text(a(1,1)+0.5,Y3p(i,j),[num2str(i),'-',num2str(j)]); ? ?%顯示小圖框里的數(shù)字位置和數(shù)值 ? ?end end xlabel('加工時間/s'); ? ? ?%橫坐標(biāo)名稱 ylabel('機(jī)器'); ? ? ? ? ? ?%縱坐標(biāo)名稱 title({[num2str(m),'*',num2str(M),'的一個最佳調(diào)度（最短完工時間為',num2str(best_fit),')']}); ? ? ?%圖形名稱 axis([0,best_fit+2,0,M+1]); ? ? ? ? %x軸，y軸的范圍 set(gca,'Box','on'); ? ? ? %顯示圖形邊框 set(gca,'YTick',0:M+1); ? ? %y軸的增長幅度 set(gca,'YTickLabel',{'';num2str((1:M)','M%d');''}); ?%顯示機(jī)器號 hold on; ? ?end %--------------------------輸入數(shù)據(jù)--------------------------------- function s = input(m,n) ? ? ?%輸入數(shù)據(jù) s = cell(m,n); s{1,1}=[1 2 3 4 5 7 8;5 3 5 3 3 10 9]; s{1,2}=[1 3 4 5 6 7 8;10 5 8 3 9 9 6]; s{1,3}=[2 4 5 6 7 8;10 5 6 2 4 5]; s{1,4}=[0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0]; s{2,1}=[1 2 3 4 5 7;5 7 3 9 8 9]; s{2,2}=[2 3 4 5 6 7 8;8 5 2 6 7 10 9]; s{2,3}=[2 4 5 6 7 8;10 5 6 4 1 7]; s{2,4}=[1 2 3 4 5 6;10 8 9 6 4 7]; s{3,1}=[1 4 5 6 7 8;10 7 6 5 2 4]; s{3,2}=[2 3 4 5 6 7;10 6 4 8 9 10]; s{3,3}=[1 2 3 4 6 8;1 4 5 6 10 7]; s{3,4}=[0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0]; s{4,1}=[1 2 3 4 5 6 7 8;3 1 6 5 9 7 8 4]; s{4,2}=[1 2 3 4 5 6 7 8;12 11 7 8 10 5 6 9]; s{4,3}=[1 2 3 4 5 6 7 8;4 6 2 10 3 9 5 7]; s{4,4}=[0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0]; s{5,1}=[1 2 3 4 5 7;3 6 7 8 9 10]; s{5,2}=[1 3 4 5 6 7;10 7 4 9 8 6]; s{5,3}=[2 3 4 5 6 7;9 8 7 4 2 7]; s{5,4}=[1 2 4 5 6 7 8;11 9 6 7 5 3 6]; s{6,1}=[1 2 3 4 5 6 8;6 7 1 4 6 9 10]; s{6,2}=[1 3 4 5 6 7 8;11 9 9 9 7 6 4]; s{6,3}=[1 2 3 4 5 7;10 5 9 10 11 10]; s{6,4}=[0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0]; s{7,1}=[1 2 3 4 5 7;5 4 2 6 7 10]; s{7,2}=[2 4 5 6 7 8;9 9 11 9 10 5]; s{7,3}=[2 3 4 5 6 8;8 9 3 8 6 10]; s{7,4}=[0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0]; s{8,1}=[1 2 3 4 6 8;2 8 5 9 4 10]; s{8,2}=[1 2 3 4 5 7;7 4 7 8 9 10]; s{8,3}=[1 2 4 5 6 7 8;9 9 8 5 6 7 1]; s{8,4}=[1 3 4 5 6 7;9 3 7 1 5 8]; ? %---------------------------建立初始種群----------------------------- function [p,TN] = initial_p(m,n,N,s,M) ? ? %建立初始種群 p = cell(N,1); ? ? ? ? ? ?%p為初始解集的機(jī)器集 TN = cell(N,1); ? ? ? ? ? ?%TN為初始解集的時間集 for i = 1:N ? ? ? ? ? ? ? ? ?%產(chǎn)生N個初始解 ? ?store_m = zeros(M,1); ? ?%用于儲存生成初始方案時的各機(jī)器數(shù)量 ? ?pz = zeros(m,n); ? ? ? ? %pz為中間儲存量，用于儲存解i的機(jī)器號，大小為m*n ? ?tz = zeros(m,n); ? ? ? ? %tz為中間儲存量，用于儲存解i的加工時間，大小為m*n ? ?for j = 1:m ? ? ? ?for k = 1:n ? ? ? ? ? ?sle = s(j,k); ? ? ? %sle為工件j的工序k的數(shù)據(jù)，第一行為可選機(jī)器數(shù)，第二行為對應(yīng)的加工時間 ? ? ? ? ? ?sle2 = cell2mat(sle); ? ?%sle為cell結(jié)構(gòu)，需要將sle用cell2mat函數(shù)轉(zhuǎn)換為double類型 ? ? ? ? ? ?b = size(sle2,2); ? ? ? %數(shù)據(jù)中有0數(shù)組，所以需要判斷 ? ? ? ? ? ?if b == 0 ? ? ? ? ? ? ? ?pz(j,k) = 0; ? ? ? ? ? ? ? ?tz(j,k) = 0; ? ? ? ? ? ?else ? ? ? ? ? ?c = randperm(b,1); ? %產(chǎn)生一個1到b的隨機(jī)數(shù)，用于選擇機(jī)器 ? ? ? ? ? ? ? ?if store_m(c) >= (m*n)/M ? ? ? ? ? ? ? ? ? ?c = randperm(b,1); ? ? ? ? ? ? ? ? ? ? ? ?if store_m(c) >= (m*n)/M ? ? ? ? ? ? ? ? ? ? ? ? ? ? c = randperm(b,1); ? ? ? ? ? ? ? ? ? ? ? ? ? ? if store_m(c) >= (m*n)/M ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?c = randperm(b,1); ? ? ? ? ? ? ? ? ? ? ? ? ? ? end ? ? ? ? ? ? ? ? ? ? ? ?end ? ? ? ? ? ? ? ?end ? ? ? ? ? ?store_m(c) = store_m(c)+1; ? ? ? ? ? ?pz(j,k) = sle2(1,c); ? ? %將機(jī)器賦予pz(j,k) ? ? ? ? ? ?tz(j,k) = sle2(2,c); ? ? %將加工時間賦予tz(j,k) ? ? ? ? ? ?end ? ? ? ?end ? ?end ? ?p{i} = pz; ? ?TN{i} = tz; end %---------------------------輸入各工序機(jī)器數(shù)量----------------------- function P = machine(n,M) P = zeros(n,1); for i= 1:n ? ?P(i) = M; ? ? ?%每道工序的可選機(jī)器數(shù)設(shè)為M end %-------------------------計算各染色體的適應(yīng)度----------------------- function [fit,Y1,Y2,Y3] = object(p,TN,N,P,m,n) ?%計算各染色體的適應(yīng)度 fit = zeros(N,1); Y1 = cell(N,1); Y2 = cell(N,1); Y3 = cell(N,1); ? ?for j = 1:N ? ? ? ?Y1{j} = zeros(m,n); ? ? ? ?Y2{j} = zeros(m,n); ? ? ? ?Y3{j} = zeros(m,n); ? ?end for w = 1:N ? ?X = p{w}; ? ? ? ? ? ? ? ? ?%變量初始化 ? ?T = TN{w}; ? ?[m,n] = size(X); ? ?Y1p = zeros(m,n); ? ?Y2p = zeros(m,n); ? ?Y3p = zeros(m,n); ? ?Q1 = zeros(m,1); ? ? ? ? %計算第一道工序的安排 ? ?Q2 = zeros(m,1); ? ?R = X(:,1); ? ? ? ? ? ? %取出第一道工序的機(jī)器號 ? ?Q3 = floor(R); ? ? ? ? ?%向下取整得到各工件在第一道工序使用的機(jī)器號 ? ?for k =1:P(1) ? ? ? ? ? %第一道工序的時間安排，k為機(jī)器號 ? ? ? ?pos = find(Q3 == k); ? ? %在Q3中取出用機(jī)器k加工的工件編號 ? ? ? ?lenpos = length(pos); ? ?%使用機(jī)器k的工件數(shù)量 ? ? ? ?if lenpos == 0 ? ? ? ?end ? ? ? ?if lenpos >= 1 ? ? ? ? ? ?Q1(pos(1)) = 0; ? ? ? ? ? ?Q2(pos(1)) = Q1(pos(1)) + T(pos(1),1); ? ? ? ? ? ?if lenpos >= 2 ? ? ? ? ? ? ? ?for j = 2:lenpos ? ? ? ? ? ? ? ? ? ?Q1(pos(j)) = Q2(pos(j-1)); ? ? ? ? ? ? ? ? ? ?Q2(pos(j)) = Q1(pos(j)) + T(pos(j),1); ? ? ? ? ? ? ? ?end ? ? ? ? ? ?end ? ? ? ?end ? ?end

三、運行結(jié)果

?