【车间调度】灰狼优化算法求解柔性作业车间问题【含Matlab源码 661期】

⛄一、车间调度简介

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小时。

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

相比于传统作业车间调度，柔性作业车间调度放宽了对加工机器的约束，更符合现实生产情况，每个工序可选加工机器变成了多个，可以由多个加工机器中的一个加工。比如上表中的实例，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]的机器组合。调度的目的还是总的完工时间最短（也可以是其他目标，比如机器最大负荷最短、总的机器负荷最短）

⛄二、灰狼算法简介

Grey Wolf Optimizer是Seyedali Mirjalili受大灰狼捕食策略的启发，于2014年提出的一种元启发式算法，主要模拟了搜索猎物、包围猎物和攻击猎物，源代码关注公众号后，回复"灰狼"或"GWO"获取。
启发
灰狼属于犬科动物，是食物链顶端的顶级掠食者，它们大多喜欢群居生活，每个种群平均5~12不等。特别有趣的是，它们有非常严格的社会等级，如下图所示。

图1 灰狼分层

领头狼（领导者）是一公一母，称为alphas（个人理解，之所以这么叫，是因为alpha是希腊字幕中的第一个，用来表示最前面的）。alpha狼主要负责决策狩猎、睡觉地点、起床时间等等，然后再将决策下达至整个种群，如果狼群垂下尾巴说明它们都认可。然而种群中也会存在民主行为，alpha狼也会听从种群中的其他狼。有趣的是，alpha狼不一定是最强壮的成员，但在管理团队方面是最好的。这说明一个狼群的组织和纪律远比它的力量重要。

第二层是beta。deta狼是从属狼，可公可母，用于辅助alpha狼制定决策或其他种群活动。当其中一只alpha狼去世或变老时，它就是alpha狼的最佳替补。beta狼应该遵从alpha狼，但也会命令其他低级别的狼，因而起到承上启下的作用。
最底层的是omega狼。它们扮演了"替罪羊"的角色（好惨）,必须屈服于其他领头狼，进食时也是排在最后。看起来omega狼在狼群中并不是一个重要的个体，但是一旦失去omega狼，整个狼群就会面临内部争斗和问题。
如果一头狼既不是alpha，beta，也不是omega，就称为从属狼或delta狼。delta狼必须听从于alpha和beta狼，但会支配omega狼。侦察狼、守卫狼、老狼、捕食狼和看管狼都是这一类。侦察狼负责监视领地的边界，一旦有危险就向狼群发出警告；守卫狼保护和保证狼群的安全；老狼是从alpha或beta退下来的经验丰富的狼；捕食狼帮助alpha和beta捕猎并为狼群提供食物；看管狼负责照顾狼群中的老弱病残。
除此之外，群体狩猎是灰狼的另一个有趣的社会行为，主要包括以下几个阶段：
跟踪、追逐和接近猎物；
追捕、包围和骚扰猎物，直到它停止移动；
向猎物发起攻击。

图2 灰狼捕食行为：(A)追逐、接近和跟踪猎物。(B-D)追捕、骚扰和包围猎物。(E)静止状态与攻击
数学模型和算法
首先给出社会等级、跟踪、包围和攻击猎物的数学模型，然后再提供完整的GWO算法。
社会等级
为了在设计GWO时对狼的社会等级进行数学建模，认为最合适的解是alpha(α \alphaα)，那么第二和第三最优解分别表示为beta(β \betaβ)和delta(δ \deltaδ),而剩余其他解都假定为omega(ω \omegaω)。在GWO中，通过α \alphaα、β \betaβ和δ \deltaδ来导引捕食（优化），ω \omegaω听从于这三种狼。
包围捕食
灰狼在捕食时会将猎物包围，使用下式进行表达这种行为：
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表示当前迭代次数，A ⃗ \vec{A}A和C ⃗ \vec{C}C为系数向量，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]之间的随机向量。
为了清楚地反映等式(1)和(2)的效果，图3(a)中显示了二维的位置向量以及可能的邻域，可以看出，灰狼的位置( X , Y ) (X,Y)(X,Y)可以根据猎物的位置( X ∗ , Y ∗ ) (X^*,Y)(X∗,Y∗)进行更新，通过调整A ⃗ \vec{A}A和C ⃗ \vec{C}C的值，可以在最优代理周围到达相对于当前位置的不同地方。例如，当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)这一种情况，当随机向量r 1 r_1r1和r 2 r_2r2取不同的值时，灰狼可以到达任意两点之间的位置。同时还注意到，A AA的取值不同还会决定灰狼靠近还是远离猎物，后面再详细说明。

图3 2D和3D位置向量及其可能的下一位置
狩猎
灰狼能够识别猎物的位置并包围它们，狩猎通常是是由alpha狼领导的，beta和和delta狼偶尔也会参与狩猎。然而，在一个抽象的搜索空间中，我们不知道最佳(猎物)的位置。为了在数学上模拟灰狼的狩猎行为，我们假设alpha(最佳候选解)、beta和delta狼对猎物的潜在位置有更好的了解。因此，我们保存到目前为止获得的前三个最佳解，并迫使其他搜索代理(包括omegas)根据最佳搜索代理的位置更新其位置。对此，提出以下公式：
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空间中如何根据alpha，beta和delta进行代理位置的更新。可以看到，最终位置将是圆内的一个随机位置，该圆由alpha、beta和delta定义，换句话说，alpha、beta和delta狼对猎物的位置进行估计，而其他狼则再猎物周围随机更新它们的位置。

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

攻击猎物（利用）
正如上面提到的，当猎物停止移动时，灰狼就会攻击它来完成狩猎。为了对接近猎物进行建模，需要不断降低a ⃗ \vec aa的值，那么A ⃗ \vec AA的波动范围也会降低。当A ⃗ ∈ [ − 1 , 1 ] \vec A\in[-1,1]A∈[−1,1]，搜索代理的下一位置可以是代理当前位置和猎物位置之间的任意位置。图5(a)表明当∣ A ∣ < 1 |A|<1∣A∣<1时，灰狼向猎物发起攻击。

图5 猎物攻击vs搜索猎物

搜索猎物（探索）
灰狼通常根据alpha、beta和delta狼的位置进行搜索，它们彼此分散寻找猎物，然后汇聚攻击猎物。为了数学上对分散建模，利用随机值大于1或小于-1的A ⃗ \vec AA迫使搜索代理偏离猎物，从而保证了探索。图5(b)表明当∣ A ∣ > 1 |A|>1∣A∣>1时，迫使灰狼离开猎物，希望能找到更合适的猎物。另一个支持GWO进行探索的因素是C ⃗ \vec CC，根据公式(4)，C ⃗ \vec CC的取值范围为[ 0 , 2 ] [0,2][0,2]，该分量为猎物提供随机权重，以随机强调(C>1)或弱化(C<1)
猎物在定义等式(1)中的距离时的作用。这有助于GWO在整个优化过程中表现出更随机的行为，有利于探索和避免局部最优。C CC并不是和A AA一样线性递减，特意要求C CC在任何时候都提供随机值，以便不仅在初始迭代中强调探索，而且在最终迭代中也强调探索。

C CC向量也可以被认为是在自然界中障碍物对接近猎物的影响，一般来说，自然中的障碍出现在狼的捕猎路径上，实际上阻碍了它们快速、方便地接近猎物。这就是向量C CC的作用。根据狼所处的位置，它可以随机地给猎物一个权重，从而让狼的捕食变得更加困难和遥远，反之亦然。

GWO算法
总之，搜索过程从在GWO算法中创建一个随机的灰狼种群(候选解)开始。在迭代过程中，alpha、beta和delta狼估计猎物可能的位置。每一个候选解更新它与猎物的距离。为了分别强调探索和利用，将参数a aa从2降低到0。当∣ A ⃗ ∣ > 1 |\vec A|>1∣A∣>1时,候选解有偏离猎物的倾向，当∣ A ⃗ ∣ < 1 |\vec A|<1∣A∣<1时,候选解收敛于猎物。最后，当满足结束条件时终止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

计算每个搜索代理的适应度值

X α = X_{\alpha}=Xα=最优搜索代理

X β = X_{\beta}=Xβ=第二优搜索代理

X δ = X_{\delta}=Xδ=第三优搜索代理

while(t<最大迭代次数)

for 每个搜索代理

根据等式(7)更新当前代理的位置

end for

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

计算所有搜索代理的适应度值

更新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在理论上是如何解决优化问题的：

所提出的社会等级有助于GWO在迭代过程中保存目前的最优解；
所提出的包围机制在解周围定义了一个圆形邻域，该邻域可以作为超球体扩展到更高维度；
随机参数A AA和C CC辅助候选解具有不同随机半径的超球体；
所提出的狩猎方法允许候选解确定猎物的可能位置；
a aa和A AA的适配保证了探索和利用；
参数a aa和A AA的自适应值使GWO在探索和利用之间实现平稳过渡；
随着A AA的下降，一半迭代致力于探索(∣ A ∣ ≥ 1 |A|≥1∣A∣≥1)，另一半致力于利用(∣ A ∣ < 1 |A|<1∣A∣<1)；
GWO只有两个主要参数需要调整(a aa和C CC)。

⛄三、部分源代码

%该程序用于解决柔性作业车间调度，m个工件，n道工序，其中n为最大工序数，工件的工序
%数可以少于n，加工机器数为M，每个工件的每道工序具有多个机器可以选择，对应的时间
%不同，其中初始种群的储存方式采用cell数据类型
%Version:1.3
%fileDescription:调度机器可选的柔性作业车间问题，甘特图已完善,GWO,88实例
%last edit time:2019-6-7
function GWO_Model_FJSP_1_3_8_8()
count = 600; %迭代次数
N = 50; %种群规模
m = 8; %工件数
n = 4; %工序数
M = 8; %机器数
a =2; %计算A/C协同系数的
plotif = 1; %控制程序是否进行绘图
s = input(m,n); %数据输入
[p,TN] = initial_p(m,n,N,s,M); %生成初始种群50,采用细胞结构，每个元素为84
P = machine(n,M);
FIT = zeros(count,1);
aveFIT = zeros(count,1);
X1=randperm(count); %收敛图形的横坐标X
X=sort(X1);
%------------------------输出最优解的时有用------------------------------
best_fit = 1000; %改变模型需要修改此参数
best_p = zeros(m,n);
best_TN = zeros(m,n);
Y1p = zeros(m,1);
Y2p = zeros(m,1);
Y3p = zeros(m,1);
minfit3 = 1000000000;
%-------------------------进行迭代--------------------------------------
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; %用于适应度函数的
aveFIT(i) = mean(fit); %用于适应度函数的
end
%------------------投射最佳方案数据--------------------------------------

fprintf('最优解：%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(‘最优解’,‘平均值’);
%-------------------------甘特图-----------------------------------------
figure;
w=0.5; %横条宽度
set(gcf,‘color’,‘w’); %图的背景设为白色
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]); %设置小图框四个点的x坐标
y=Y3p(i,j)+[-w/2 w/2 w/2 -w/2]; %设置小图框四个点的y坐标
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)]); %显示小图框里的数字位置和数值
end
end
xlabel(‘加工时间/s’); %横坐标名称
ylabel(‘机器’); %纵坐标名称
title({[num2str(m),‘*’,num2str(M),‘的一个最佳调度（最短完工时间为’,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’);’'}); %显示机器号
hold on;
end
%--------------------------输入数据---------------------------------
function s = input(m,n) %输入数据
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为初始解集的机器集
TN = cell(N,1); %TN为初始解集的时间集
for i = 1:N %产生N个初始解
store_m = zeros(M,1); %用于储存生成初始方案时的各机器数量
pz = zeros(m,n); %pz为中间储存量，用于储存解i的机器号，大小为mn
tz = zeros(m,n); %tz为中间储存量，用于储存解i的加工时间，大小为mn
for j = 1:m
for k = 1:n
sle = s(j,k); %sle为工件j的工序k的数据，第一行为可选机器数，第二行为对应的加工时间
sle2 = cell2mat(sle); %sle为cell结构，需要将sle用cell2mat函数转换为double类型
b = size(sle2,2); %数据中有0数组，所以需要判断
if b == 0
pz(j,k) = 0;
tz(j,k) = 0;
else
c = randperm(b,1); %产生一个1到b的随机数，用于选择机器
if store_m© >= (mn)/M
c = randperm(b,1);
if store_m© >= (mn)/M
c = randperm(b,1);
if store_m© >= (m*n)/M
c = randperm(b,1);
end
end
end
store_m© = store_m©+1;
pz(j,k) = sle2(1,c); %将机器赋予pz(j,k)
tz(j,k) = sle2(2,c); %将加工时间赋予tz(j,k)
end
end
end
p{i} = pz;
TN{i} = tz;
end
%---------------------------输入各工序机器数量-----------------------
function P = machine(n,M)
P = zeros(n,1);
for i= 1:n
P(i) = M; %每道工序的可选机器数设为M
end
%-------------------------计算各染色体的适应度-----------------------
function [fit,Y1,Y2,Y3] = object(p,TN,N,P,m,n) %计算各染色体的适应度
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); %取出第一道工序的机器号
Q3 = floor®; %向下取整得到各工件在第一道工序使用的机器号
for k =1:P(1) %第一道工序的时间安排，k为机器号
pos = find(Q3 == k); %在Q3中取出用机器k加工的工件编号
lenpos = length(pos); %使用机器k的工件数量
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

⛄四、运行结果

⛄五、matlab版本及参考文献

1 matlab版本
2014a

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

3 备注
简介此部分摘自互联网，仅供参考，若侵权，联系删除