【基于遗传算法改进的粒子群GA-PSO算法优化shubert函数及MATLAB编程实现，测试函数shubert(十)

测试函数shubert（十）

shubert函数属于周期性多峰函数，如图1所示拥有多个全局最优值，如图二所示在一个周期内只有一个全局最优值，局部最优解较多，适合测试算法的收敛性能，粒子群算法是一种收敛速度较快的算法，运算速度快，但是粒子群算法容易陷入局部最优，有些时候会导致收敛慢，或者不收敛，本文用标准粒子群算法进行改进，利于遗传算法的变异算子，交叉算子，对标准粒子群进行改进，，用基于遗传算法改进的粒子群GA-PSO求解，如有疑问，欢迎大家留言交流！
函数图像如下


MATLAB编程shubert代码如下：

clc
clear
close all
x = -2:0.1:2;
y  = -2:0.1:2;
x = -10:0.1:10;
y  = -10:0.1:10;
[x,y] = meshgrid(x,y);
[m,n] = size(x);
z = zeros(m,n);
for ii = 1:m
    for jj = 1:n
        xx = [x(ii,jj) y(ii,jj)];
        z(ii,jj) = shubertfun(xx);
    end
end
figure
surf(x,y,z)
xlabel('x1')
ylabel('x2')
zlabel('z')
axis([-2 2 -2 2 -200 200])
% axis([-10 10 -10 10 -200 200])
% shading interp
title('Shubert Function')
set(gca,'fontsize',12)
colormap jet
function [ out] = shubertfun( x )
x1 = x(1);
x2 = x(2);
sum1 = 0;
sum2 = 0;

for ii = 1:5
	new1 = ii * cos((ii+1)*x1+ii);
	new2 = ii * cos((ii+1)*x2+ii);
	sum1 = sum1 + new1;
	sum2 = sum2 + new2;
end

out = sum1 * sum2;

end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

遗传算法改进粒子群GA-PSO算法的原理
粒子群优化算法（PSO）又翻译为粒子群算法、微粒群算法、或微粒群优化算法，PSO是由Kennedy和Eberhart共同提出，最初用于模拟社会行为，作为鸟群中有机体运动的形式化表示。自然界中各种生物体均具有一定的群体行为，Kennedy和Eberhart的主要研究方向之一是探索自然界生物的群体行为，从而在计算机上构建其群体模型。PSO是一种启发式算法，因为它很少或没有对被优化的问题作出假设，并且能够对非常大候选解决方案空间进行搜索。PSO算法初始化为一群随机的粒子，然后通过多次迭代找到最优解。每一次的迭代过程中，粒子通过本身所找到的最优解（被成为个体极值）和整个种群目前找到的最优解（被称为全局极值）来更新自己。也可以使用粒子本身的邻居（被称为局部极值）的极值来更新自己，粒子群算法收敛快，但是容易陷入局部最优，针对次问题，本文通过判断每个粒子群的适应度，对适应度低的粒子，利用遗传算法的交叉算子，变异算子进行交叉运算和变异运算，增强了算法的全局收敛能力，能防止种群陷入局部最优。

GA-PSO粒子群算法的主要参数
一、种群个数popsize，既算法中粒子的个数；
二、最大迭代次数gen，既算法迭代gen次后停止迭代；
三、种群维度dim，既需要优化的自变量个数；
四、种群位置pop，既每个粒子群的对应的自变量的值，一个粒子对应一组自变量，相当于一个解；
五、种群速度v，既粒子群每次迭代更新的飞行速度，粒子群位置更新的步长；
六、种群全局最优值gbest，既迭代过程中曾经出现的最优解，包括最优位置和对应的目标函数值；
七、个体最优，既每个粒子迭代过程中单个体曾经出现的个体最优解，，包括个体最优位置和对应的目标函数值；
八、个体学习因子c1，既个体最优解对粒子群飞行的影响能力；
九、全局学习因子c2，既全局最优值对粒子群飞行的影响能力；
十、惯性权重w，既个体位置所占的权重，权重越大，粒子群收敛越慢，全局搜索能力越强；
十一、变异概率pm，既遗传算子的变异算子的概率。
十二、交叉概率pc，既遗传算法的交叉算子的概率

遗传算法改进粒子群GA-PSO算法流程图

遗传算法优化粒子群GA-PSO算法，MATLAB编程代码如下：

clc
clear
close all
warning off
set(0,'defaultfigurecolor','w') 
%x y
xmax = [5.12 5.12];
xmin = [-5.12 -5.12];
vmax = 0.2*xmax;
vmin = -vmax;
fun = @shubert;
m=2;
%程序初始化
% global popsize;     %种群规模
gen=20;     %设置进化代数
popsize=300;     %设置种群规模大小
best_in_history(gen)=inf;   %初始化全局历史最优解
best_in_history(:)=inf;   %初始化全局历史最优解
best_fitness=inf;
fz = zeros(gen,5);
%设置种群数量
pc=0.5;                       %交叉概率选择，0和1之间
pm=0.1;                    %变异概率选择，0和1之间
pop1 = zeros(popsize,m);
pop2 = zeros(popsize,m);
pop3 = zeros(popsize,m);
pop6 = zeros(gen,m);%存储解码后的每代最优粒子
pop7 = zeros(popsize,m);%存储更新解码后的粒子的位置
for ii1=1:popsize
    pop1(ii1,:)=funx(xmin,xmax,m);   %初始化种群中的粒子位置，
    pop3(ii1,:)=pop1(ii1,:);  %初始状态下个体最优值等于初始位置
    pop2(ii1,:)=funv(vmax,m);    %初始化种群微粒速度，
    pop4(ii1,1)=inf;
    pop5(ii1,1)=inf;
end
pop0=pop1;
xmax = [5.12 5.12];
xmin = [-5.12 -5.12];
c1=2;
c2=2;

gbest_x=pop1(end,:);
% pop1(1:size(num,1),:) = num;   %全局最优初始值为种群第一个粒子的位置

for exetime=1:gen
    ww = 0.7*(gen-exetime)/gen+0.2;
    for ii4=1:popsize
        
        pop2(ii4,:)=(ww*pop2(ii4,:)+c1*rand(1,m).*(pop3(ii4,:)-pop1(ii4,:))+c2*rand(1,m).*(gbest_x-pop1(ii4,:)));    %更新速度
        for jj = 1:m
            if  pop2(ii4,jj)<vmin(jj)
                pop2(ii4,jj)=vmin(jj);
            elseif pop2(ii4,jj)>vmax(jj)
                pop2(ii4,jj)=vmax(jj);
            end
        end
    end
    %更新粒子位置
    for ii5=1:popsize
        pop1(ii5,:)=pop1(ii5,:)+pop2(ii5,:);
        for jj2 = 1:m
            if  pop1(ii5,jj2)>xmax(jj2)
                pop1(ii5,jj2) = xmax(jj2);
            elseif pop1(ii5,jj2)<xmin(jj2)
                pop1(ii5,jj2)=xmin(jj2);
            end
        end
        
%         if  rand>0.85
%             k=ceil(m*rand);
%             pop1(ii5,k) = (xmax( k)-xmin(k)).*rand(1,1)+xmin(k);
%         end
%         if  pop5(ii5)>sum(pop5)/popsize
%             pop1(ii5,:) = (xmax(1,m)-xmin(1,m)).*rand(1,m)+xmin(1,m);
%         end
    end
    for jj2 = 1:m
        if  pop1(ii5,jj2)>xmax(jj2)
            pop1(ii5,jj2) = xmax(jj2);
        elseif pop1(ii5,jj2)<xmin(jj2)
            pop1(ii5,jj2)=xmin(jj2);
        end
    end
    if  exetime-1>0
        plot(1:length(best_in_history(1:exetime-1)),best_in_history(1:exetime-1));
        xlabel('迭代次数')
        ylabel('适应度')
        title('遗传算法粒子群GA-PSO算法')
        hold on;
        pause(0.1)
    end
    
    pop1(end,:) = gbest_x;
    %计算适应值并赋值
    for ii3=1:popsize
        [my,mx] = fun2(pop1(ii3,:));
        %          [my,mx] = fun2(gbest_x,num,xmax,xmin);
        pop5(ii3,1)=my;
        pop7(ii3,:) = mx;
        if  pop4(ii3,1)>pop5(ii3,1)    %若当前适应值优于个体最优值，则进行个体最优信息的更新
            pop4(ii3,1)=pop5(ii3,1);          %适值更新
            pop3(ii3,:)=pop1(ii3,:);      %位置坐标更新
        end
    end
    
    %计算完适应值后寻找当前全局最优位置并记录其坐标
    if  best_fitness>min(pop4(:,1))
        best_fitness=min(pop4(:,1))  ;  %全局最优值
        ag = [];
        ag =find(pop4(:,1)==min(pop4(:,1)));
        
        gbest_x(1,:)=(pop1(ag(1),:));    %全局最优粒子的位置

        pop6(exetime,:) = pop7(ag(1),:);
    else
        fz(exetime,:) = fz(exetime-1,:);
        if exetime>1
            pop6(exetime,:) = pop6(exetime-1,:);
        end
    end
    best_in_history(exetime)=best_fitness;  %记录当前全局最优
    
    meany = mean(pop5);
    %变异
    for hh = 1:popsize
        if pop5(hh)>meany
            if  rand()>pm
                k=ceil(rand());
                pop1(hh,k) = (xmax( k)-xmin(k)).*rand(1,1)+xmin(k);
            end
        end
    end
    %交叉
    pany = [];
    for hh = 1:popsize
        if pop5(hh)>meany
            if isempty(pany)
                pany =[pop1(hh,:) hh];
            else
                panx = [pop1(hh,:) hh];
                if  rand()>pc
                    k=ceil(m*rand());
                    pop1(pany(end),k) = panx(k);
                    pop1(panx(end),k) =pany(k);
                    pany = [];
                    
                end
            end
        end
    end
end

x = -2:0.1:2;
y  = -2:0.1:2;
x = -10:0.1:10;
y  = -10:0.1:10;
[x,y] = meshgrid(x,y);
[m,n] = size(x);
z = zeros(m,n);
for ii = 1:m
    for jj = 1:n
        xx = [x(ii,jj) y(ii,jj)];
        z(ii,jj) = shubertfun(xx);
    end
end
figure
surf(x,y,z)
hold on
xlabel('x1')
ylabel('x2')
zlabel('z')
axis([-2 2 -2 2 -200 200])
% axis([-10 10 -10 10 -200 200])
% shading interp
title('Shubert Function')
set(gca,'fontsize',12)
colormap jet
plot(pop0(:,1),pop0(:,2),'ro','MarkerFaceColor','r')
xlabel('X')
ylabel('Y')
title('初始种群')
set(gca,'fontsize',12)
view([-130 40])

figure
surf(x,y,z)
hold on
xlabel('x1')
ylabel('x2')
zlabel('z')
axis([-2 2 -2 2 -200 200])
% axis([-10 10 -10 10 -200 200])
% shading interp
title('Shubert Function')
set(gca,'fontsize',12)
colormap jet
plot(pop1(:,1),pop1(:,2),'ro','MarkerFaceColor','r')
xlabel('X')
ylabel('Y')
title('收敛后的种群')
set(gca,'fontsize',12)
view([-130 40])
function [out,xy]= fun2(xy)  
xmax = [5.12 5.12 ];
xmin = [-5.12 -5.12 ];
for ii = 1:2
    if xy(ii)>xmax(ii)
        xy(ii)=xmax(ii);
    elseif xy(ii)<xmin(ii);
        xy(ii)=xmin(ii);
    end
end


x1 = xy(1);
x2 = xy(2);
sum1 = 0;
sum2 = 0;

for ii = 1:5
	new1 = ii * cos((ii+1)*x1+ii);
	new2 = ii * cos((ii+1)*x2+ii);
	sum1 = sum1 + new1;
	sum2 = sum2 + new2;
end

out = sum1 * sum2;

end
    
function x=funv(vmax,m)

    x = 2.*vmax(1:m).*rand(1,m)-vmax(1:m);
    
end

function x=funx(xmin,xmax,m)

    x  = (xmax(1:m)-xmin(1:m)).*rand(1,m)+xmin(1:m);
    
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241

遗传算法改进粒子群GA-PSO算法优化shubert函数，效果图如下：

从效果图看，收敛图比较温和，收敛后种群没有完全聚集，有一部分粒子在全局寻优，大部分粒子在最优解附近寻优，既保证了对最优解的精益求精，又可以避免陷入局部最优，全局搜索能力不降低。