A list,
Grey Wolf Optimizer is a meta-heuristic algorithm proposed by Seyedali Mirjalili in 2014, which is inspired by the hunting strategy of the Big Grey Wolf. It mainly simulates searching, surrounding and attacking prey. After following the source code of the public account, it can reply to “Grey Wolf” or “GWO” for obtaining.
Inspired by the
Gray wolves belong to the canid family and are the top predators at the top of the food chain. They tend to live in groups with an average of 5 to 12 wolves per group. What is particularly interesting is that they have very strict social hierarchies, as shown below.
FIG. 1 Stratification of gray wolves
The alpha Wolf (the leader) is a male and a female, called alphas. Alpha wolves are responsible for making decisions about hunting, where to sleep, when to wake up, and so on, and then making decisions about the whole pack. If the wolves have their tails down, they all agree. However, there is also democratic behavior in the population, and alpha wolves obey other wolves in the population. Interestingly, the Alpha Wolf was not necessarily the strongest member, but was the best at managing the team. This shows that a pack’s organization and discipline are far more important than its strength.
The second layer is beta. Deta wolves are subordinate wolves, either male or female, used to assist alpha wolves in decision-making or other population activities. When one of the Alpha wolves dies or grows old, he is the best replacement for the Alpha Wolf. The beta Wolf should obey the alpha Wolf, but also give orders to other lower-level wolves, thus acting as a link between the preceding and the following.
At the bottom are the Omega wolves. They play the role of “scapegoat”, having to submit to the other leaders and come last to feed. It may seem that omega wolves are not an important individual in the pack, but without omega wolves, the whole pack will face internal strife and problems.
If a Wolf is neither alpha, beta or omega, it is called a subordinate Wolf or delta Wolf. Delta wolves must obey alpha and beta wolves, but will dominate omega wolves. Scout wolves, guard wolves, old wolves, predator wolves, and custodian wolves are all in this category. Scout wolves watch the boundaries of their territory and warn wolves of danger; To protect and ensure the safety of wolves; An old Wolf is an experienced Wolf who has retired from alpha or beta; Predator wolves help alpha and beta hunt and provide food for wolves; The Wolf warden is responsible for the care of the sick and old of the pack.
In addition, pack hunting is another interesting social behavior of gray wolves, which mainly includes the following stages:
To track, chase, and approach prey;
Chase, surround, and harass prey until it stops moving;
Attack the prey.
FIG. 2 predation behavior of gray wolves :(A) chasing, approaching and stalking prey. (B-d) Hunt, harass and surround prey. (E) Static states and attacks
Mathematical models and algorithms
Mathematical models of social hierarchy, tracking, encircling, and attacking prey are presented first, followed by the complete GWO algorithm.
Social hierarchy
In order to mathematically model the social hierarchy of wolves when designing GWO, it was thought that the most suitable solution was alpha(α \alphaα), Then the second and third best solutions are expressed as beta(β \betaβ) and delta(δ \deltaδ) respectively, while the remaining solutions are assumed to be omega(ω \omegaω). In GWO, omega omega listens to all three wolves through alpha \alpha, β \betaβ and δ \deltaδ to guide predation (optimization).
Surrounded by a predator
Gray wolves surround their prey when they hunt. This behavior can be expressed by using the following formula:
D ⃗ = ∣ ⃗ C ⋅ X p – (t) – X ⃗ (t) ∣ \ vec (1) {D} = | \ vec {C} \ cdot \ overrightarrow {X_ {p}} (t) – \ vec {X} (t) | \ tag {1} D = ∣ C ⋅ Xp (t) – X (t) ∣ (1)
⃗ X (t + 1) = p – > X (t) – A ⃗ ⃗ ⋅ D (2) \ vec {X} (t + 1) = \ overrightarrow {X_ {p}} (t) – \ vec {A} \ cdot \ \ vec {D} tag (t + 1) = {2} X Xp (t) – A ⋅ D (2)
Where t tt represents the current iteration number, A ⃗ \vec{A}A and C ⃗ \vec{C}C are coefficient vectors, Xp → \overrightarrow{X_{p}}Xp is the position vector of prey, X ⃗ \vec{X}X is the position vector of gray Wolf.
Vectors A ⃗ \vec{A}A and C ⃗ \vec{C}C are computed as follows:
A ⃗ = 2 A ⃗ ⋅ – A ⃗ r ⃗ 1 (3) \ vec \ {A} = 2 vec {A} \ cdot \ vec {r}{1} – \ vec {a} \ tag {3} a = 2 a ⋅ r1 – a (3)
C ⃗ =2 ⋅ r 2 → (4) \vec{C}=2 \cdot \overrightarrow{r{2}} \ tag {4} C ⋅ r2 = 2 (4)
Where a ⃗ \vec{a} the components of a decrease linearly from 2 to 0 in the iterative process, r ⃗ 1 \vec{r}R1 and r 2 → \overrightarrow{rR2 is the random vector between [0,1].
In order to clearly reflect the effects of equations (1) and (2), two-dimensional position vectors and possible neighborhoods are shown in FIG. 3(a). It can be seen that the positions of gray wolves (X,Y)(X,Y) can be based on the positions of prey (X ∗,Y ∗)(X *,Y)(X∗,Y∗) is updated, and by adjusting the values of A ⃗ \vec{A}A and C ⃗ \vec{C}C, different places can be reached around the optimal agent relative to the current position. For example, when A ⃗ = (0, 1) \ vec {A} = (0, 1) ⃗ = A = (0, 1) and C (1, 1) \ vec {C} = (1, 1) = C (1, 1), the gray Wolf to A new location for (∗ X – X, Y ∗) (X * – X, Y∗) (X – X, Y ∗). It’s similar in three dimensions. Note, here’s two pictures only shows A ⃗ = (0, 1) \ vec {A} = (0, 1) ⃗ = A = (0, 1) and C (1, 1) \ vec {C} = (1, 1) = C (1, 1) of this kind of circumstance, When the random vectors R 1, R_1R1 and R 2, r_2R2 take different values, the gray Wolf can reach the position between any two points. It is also noted that different values of A and AA also determine whether gray wolves approach or stay away from the prey, which will be explained in detail later.
FIG. 3 2D and 3D position vectors and their possible next positions
hunting
Gray wolves can recognize the location of prey and surround them. Hunting is usually led by alpha wolves, with beta and Delta wolves occasionally participating in hunting. However, in an abstract search space, we don’t know the best location. To mathematically model gray Wolf hunting behavior, we assume that alpha(best candidate), beta, and Delta wolves have a better understanding of the potential location of prey. Therefore, we save the first three best solutions obtained so far and force other search agents (including Omegas) to update their positions based on the location of the best search agent. In this regard, the following formula is proposed:
D alpha – = C ⋅ ⃗ 1 X ∣ alpha – – X ⃗ ∣, D beta – = ∣ C ⋅ ⃗ 2 X beta – – X ⃗ ∣, Delta – > = ∣ C 3 D – ⋅ X delta – – 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 alpha = {5} D ∣ ∣ ∣ C1 ⋅ alpha – X X ∣ ∣ ∣, D beta = ∣ ∣ ∣ C2 ⋅ beta – X X ∣ ∣ ∣, D the delta = ∣ C3 ⋅ delta – X X ∣ (5)
– > = X X 1 alpha – – A 1 – ⋅ alpha – (D), X – > = 2 X beta – 2 – A – ⋅ beta – (D), X 3 – > = X delta – – A 3 – ⋅ delta – (D) (6) (A) \ 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) the X1 = X alpha – A1 ⋅ alpha (D), X2 = X beta – A2 ⋅ beta (D), the X3 = X delta – A3 ⋅ delta (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)
Figure 4 shows how agent locations are updated in 2D space according to alpha, beta, and delta. As you can see, the final position will be a random position within the circle defined by alpha, beta, and Delta. In other words, alpha, beta, and Delta wolves estimate the location of prey, while other wolves randomly update their position around prey.
Figure 4. Location update schematic in GWO
To attack (use) prey.
As mentioned above, when the prey stops moving, the gray Wolf attacks it to complete the hunt. In order to model the approach of prey, the value of a ⃗ \vec AA needs to be constantly reduced, so the fluctuation range of a ⃗ \vec AA will also be reduced. When A ⃗ ∈[−1,1] \vec A\in[-1,1]A∈[−1,1], the next position of the search agent can be any position between the current position of the agent and the position of the prey. Figure 5 (a) shows that when ∣ a ∣ < 1 | a | < 1 ∣ ∣ a < 1, the gray Wolf to attack prey.
Figure 5. Prey attack vs hunt prey
Hunting (exploring) Gray wolves usually search by the location of alpha, beta, and Delta wolves, which disperse from each other in search of prey and then converge to attack it. To mathematically model dispersion, A ⃗ \vec AA with A random value greater than 1 or less than -1 is used to force the search agent to deviate from the prey, thus ensuring exploration. Figure 5 (b) shows that when ∣ A ∣ > 1 | A | 1 ∣ A ∣ > > 1, forced the wolves to leave their prey, hoping to find A more suitable prey. Another factor supporting GWO exploration is C ⃗ \vec CC. According to Formula (4), the value range of C ⃗ \vec CC is [0,2][0,2] [0,2]. This component provides random weight for prey, To emphasize (C>1) or to weaken (C<1) the role of prey in defining distance in equation (1). This helps GWO exhibit more random behavior throughout the optimization process, facilitating exploration and avoiding local optimizations. C CC does not decrease linearly like A AA, and is deliberately required to provide random values at all times in order to emphasize exploration not only in the initial iteration, but also in the final iteration.
The C CC vector can also be considered as the effect of obstacles on approaching prey in nature. Generally speaking, obstacles in nature appear in the hunting path of wolves and actually prevent them from approaching prey quickly and conveniently. That’s what the vector C, CC does. Depending on where the Wolf is, it can randomly assign a weight to its prey, making the Wolf’s hunt more difficult and distant, and vice versa.
GWO algorithm In summary, the search process begins by creating a random population of gray wolves (candidate solutions) in the GWO algorithm. During iteration, alpha, beta, and Delta wolves estimate the likely location of prey. Each candidate solution updates its distance from its prey. In order to emphasize exploration and utilization respectively, the parameter A aa is reduced from 2 to 0. When ∣ A ⃗ ∣ > 1 | | \ vec A > 1 ∣ A ∣ > 1, the candidate solutions have A tendency to deviate from its prey, when ∣ A ⃗ ∣ < 1 | | \ vec A < 1 ∣ ∣ A < 1, the candidate solution converge to its prey. Finally, the GWO algorithm is terminated when the end condition is satisfied. The pseudocode for GWO is as follows.
I (I =1,2… n) X_i(I =1,2… , n) Xi (I = 1, 2,… ,n)
Initialize a,A,C A, A,Ca, a, C
Calculate the fitness value for each search agent
Xα= X_{\alpha}=Xα= optimal search agent
Xβ= X_{\beta}=Xβ= second best search agent
Xδ= X_{\delta}=Xδ= third best search agent
While (t< maximum number of iterations)
For each search agent
Update the location of the current agent according to equation (7)
end for
Update a,A,C A, a, Ca, a, C
Calculates fitness values for all search agents
Update X alpha X_ {\ alpha} alpha, beta X_ X X {\ \ beta} beta, X X delta X_ {\} delta X delta
KaTeX Parse error: Expected ‘EOF’, got ‘&’ at position 1: &̲emsp; T = t +…
end while
Return the X alpha X_ {\ alpha} X alpha
Here’s how GWO theoretically solves optimization problems:
The proposed social hierarchy helps GWO to preserve the current optimal solution in the iterative process. The proposed enveloping mechanism defines a circular neighborhood around the solution, which can be extended to higher dimensions as a hypersphere. The auxiliary candidate solutions of random parameters A AA and C CC are hyperspheres with different random radii. The proposed hunting method allows candidate solutions to determine the possible location of prey. The adaptation of A AA and A AA ensures exploration and utilization; The adaptive values of parameters A AA and A AA make GWO achieve a smooth transition between exploration and utilization. With A decline in AA, half iteration is committed to explore (∣ A ∣ p 1 | | A ∣ A ∣ acuity 1 or higher), the other half is committed to using (∣ A ∣ < 1 | A | < 1 ∣ ∣ A < 1); GWO has only two major parameters to adjust (A AA and C CC).
Ii. Source code
%该程序用于解决柔性作业车间调度，m个工件，n道工序，其中n为最大工序数，工件的工序
%数可以少于n，加工机器数为M，每个工件的每道工序具有多个机器可以选择，对应的时间
%不同，其中初始种群的储存方式采用cell数据类型
%Version:1.3
%fileDescription:调度机器可选的柔性作业车间问题，甘特图已完善,GWO,8*8实例
%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,采用细胞结构，每个元素为8*4
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的机器号，大小为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的数据，第一行为可选机器数，第二行为对应的加工时间
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(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); %将机器赋予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(R); %向下取整得到各工件在第一道工序使用的机器号
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
Copy the code3. Operation results
