吸引-排斥优化算法

论文参考

Karol Cymerys, Mariusz Oszust,Attraction–Repulsion Optimization Algorithm for Global Optimization Problems,Swarm and Evolutionary Computation,Volume 84, 2024, 101459, https://doi.org/10.1016/j.swevo.2023.101459.

demoGlobalOptimization.m

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
% Attraction - Repulsion Optimization Algorithm (AROA) for Global Optimization Problems
% K. Cymerys, M. Oszust 

clear all
clc

SearchAgents_no = 40;    %   Number of search agents  
Max_FES = 3000;          % Maximum number of fitness function evaluations
NoRepeats = 3;           % 
  
fobj = @(x)F22(x);       %just a function

dim = 100;
lb = -5.12 * ones(1,dim);
ub = 5.12 * ones(1,dim);
%%fmin = 0; 

outputF=[];
outputX=[];
    for k = 1:NoRepeats

        tic;

            [Best_score, Best_pos, cg_curve] = AROA(SearchAgents_no,dim,lb,ub, fobj,Max_FES);

        toc 

       Best_score            
       outputF = [outputF;Best_score];
       outputX = [outputX;Best_pos];

    end
[Y,I] = min(outputF) ;
Y
meanV = mean(outputF)
stdV = std(outputF)

定义参数与设置：

SearchAgents_no = 40：设定搜索代理（即种群个体）的数量为40。

Max_FES = 3000：设置最大适应度函数评估次数为3000，作为算法的迭代终止条件。

NoRepeats = 3：设置算法重复运行次数为3次，以获取多次运行结果的统计信息。

定义目标函数与问题维度：

fobj = @(x)F22(x)：将目标函数F22封装为匿名函数fobj，用于计算给定解x的目标函数值。

dim = 100：设定问题的维度为100，即目标函数的自变量为100维向量。

lb = -5.12 * ones(1,dim)：设定搜索空间下界为-5.12的100维全1向量，即每个自变量的最小值均为-5.12。

ub = 5.12 * ones(1,dim)：设定搜索空间上界为5.12的100维全1向量，即每个自变量的最大值均为5.12。

%%fmin = 0;：注释行，未启用（可能原本打算设置目标函数最小值为0，但被注释掉了）。

初始化结果存储数组：

outputF=[]：创建空矩阵outputF，用于存储每次运行得到的最佳目标函数值。

outputX=[]：创建空矩阵outputX，用于存储每次运行得到的最佳解（即最优自变量值）。

主循环：重复运行AROA算法并记录结果

for k = 1:NoRepeats：对于设定的重复次数（3次），执行以下操作。 1. tic：开始计时，记录本次AROA算法运行时间。 2. [Best_score, Best_pos, cg_curve] = AROA(SearchAgents_no,dim,lb,ub, fobj,Max_FES);：调用AROA函数，输入参数包括搜索代理数量、问题维度、搜索空间上下界、目标函数以及最大适应度函数评估次数。函数返回本次运行的最佳目标函数值（Best_score）、最佳解（Best_pos）以及收敛曲线数据（cg_curve）。 3. toc：结束计时，输出本次AROA算法运行时间。 4. Best_score：显示本次运行的最佳目标函数值。 5. outputF = [outputF;Best_score];：将本次最佳目标函数值添加到outputF矩阵中。 6. outputX = [outputX;Best_pos];：将本次最佳解添加到outputX矩阵中。

统计分析多次运行结果：

[Y,I] = min(outputF)：从outputF矩阵中找出最小目标函数值（Y）及其对应的索引（I）。

Y：显示最小目标函数值。

meanV = mean(outputF)：计算outputF矩阵中所有最佳目标函数值的平均值（meanV）。

stdV = std(outputF)：计算outputF矩阵中所有最佳目标函数值的标准差（stdV）。

综上所述，该代码使用AROA算法对具有100维搜索空间的问题进行全局优化，重复运行3次，并记录每次运行的最佳解和目标函数值。最后，对多次运行结果进行统计分析，输出最小目标函数值、平均值及标准差。

F22.m

1
2
3
4
5
6
7
8
function  y  = F22(x)
 
sum = 0;
for ii = 1:length(x)	 
	sum = sum + x(ii)^2;
end

y = sum;

AROA.m

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
function [fbest, xbest, Convergence_curve] = AROA(N, dim, lb, ub, fobj, maxEvals)
    
    % Algorithm parameters definition
    c = 0.95;
    fr1 = 0.15;
    fr2 = 0.6;
    p1 = 0.2;  
    p2 = 0.8; 
    Ef = 0.4;
    tr1 = 0.9;
    tr2 = 0.85; 
    tr3 = 0.9;
    % Algorithm parameters definition
    
    tmax = ceil((maxEvals - N)/(2*N));
    evalCounter = 0;

    Convergence_curve = zeros(1,tmax);
    Xmin = repmat(ones(1,dim).*lb,N,1);
    Xmax = repmat(ones(1,dim).*ub,N,1);


    % random initialization - Eq (3)
    X = rand(N,dim).*(ub-lb) + lb;
    [X, F, evalCounter] = evaluate_population(X, fobj, ub, lb, evalCounter, maxEvals);
    [fbest, ibest] = min(F);
    xbest = X(ibest,:);
    % random initialization - Eq (3)

    X_memory = X;
    F_memory = F;

    % Main loop
    for t=1:tmax	
        D = squareform(pdist(X, 'squaredeuclidean'));  % Eq (4) 
        m = tanh(t, tmax, [-2, 7]);   % Eq (11)    

        for i=1:N
           Dimax = max(D(i,:));
           k = floor((1-t/tmax)*N)+1;  % Eq (9)
           [~, neighbors] = sort(D(i,:));
           
           % Attraction-Repulsion operator % Eq (6)
           delta_ni = zeros(1,dim);
           for j=neighbors(1:k)
                I = 1 - (D(i,j)/Dimax);  % Eq (7)
                s = sign(F(j)-F(i));  % Eq (8)
                delta_ni = delta_ni + c*(X_memory(i,:)-X_memory(j,:))*I*s;
           end
           ni = delta_ni/N;
           % Attraction-Repulsion operator % Eq (6)

           % Attraction to best solusion Eq (10)
           if rand < p1
               bi = m*c.*(rand(1,dim).*xbest - X_memory(i,:));
           else
               bi = m*c.*(xbest - X_memory(i,:));
           end
           % Attraction to best solusion Eq (10)

           % Local search operators Eq (15)
           if rand < p2
               if rand > 0.5*t/tmax + 0.25
                   u1 = rand(1, dim) > tr1;
                   ri = u1.*random('Normal', zeros(1,dim), fr1*(1-t/tmax)*(ub-lb));  % Eq (12)
               else
                   % Eq (13)
                   u2 = rand(1,dim) > tr2;
                   w = index_roulette_wheel_selection(F, k);
                   Xw = X_memory(w,:);
                   if rand < 0.5
                       ri = fr2*u2.*(1-t/tmax).*sin(2*pi*rand(1,dim)).*abs(rand(1,dim).*Xw-X_memory(i,:));
                   else
                       ri = fr2*u2.*(1-t/tmax).*cos(2*pi*rand(1,dim)).*abs(rand(1,dim).*Xw-X_memory(i,:));
                   end
                   % Eq (13)
               end
           else
               u3 = rand(1,dim) > tr3;
               ri = u3.*(2*rand(1,dim)-ones(1,dim)) .* (ub-lb);  % Eq (14)
           end
           % Local search operators Eq (15)

           X(i,:) = X(i,:) + ni + bi + ri;  % Eq(16)
        end
        
        [X, F, evalCounter] = evaluate_population(X, fobj, ub, lb, evalCounter, maxEvals);
        [fbest_candidate, ibest_candidate] = min(F);

        if fbest_candidate < fbest
            fbest = fbest_candidate;
            xbest = X(ibest_candidate, :);
        end

	    [X, F] = memory_operator(X, F, X_memory, F_memory);  % Eq (18)
        X_memory = X;
        F_memory = F; 
			
        % Eq (17)			 
        CF=(1-t/tmax)^3 ;			 
        if rand < Ef
            u4 = rand(N,dim) < Ef;                                                                                              
            X = X + CF*(u4.*(rand(N,dim).*(Xmax-Xmin) + Xmin));
        else
            r7 = rand();		
            X = X + (CF*(1-r7) + r7)*(X(randperm(N),:) - X(randperm(N),:));
        end
		% Eq (17)

        [X, F, evalCounter] = evaluate_population(X, fobj, ub, lb, evalCounter, maxEvals);
        [fbest_candidate, ibest_candidate] = min(F);

        if fbest_candidate < fbest
            fbest = fbest_candidate;
            xbest = X(ibest_candidate, :);
        end
	 
        [X, F] = memory_operator(X, F, X_memory, F_memory);  % Eq (18)
        X_memory = X;
        F_memory = F; 

  	    Convergence_curve(t) = fbest;
    end		
end

function [X, F, evalCounter] = evaluate_population(X, fobj, ub, lb, evalCounter, maxEvals)
    N = size(X,1);
    F = Inf(N,1);
    X = max(lb, min(ub, X)); % Check space bounds
    
    for i=1:N
        if evalCounter >= maxEvals
            break
        end
        F(i) = fobj(X(i,:));
        evalCounter = evalCounter + 1;
    end
end

function [X, F] = memory_operator(X, F, X_memory, F_memory)
    dim = size(X, 2);
    Inx = F_memory < F;
    Indx = repmat(Inx,1,dim);
    X = Indx.*X_memory + ~Indx.*X;
    F = Inx.*F_memory + ~Inx.*F;
end

function [y] = tanh(t, tmax, range)
    z = 2*(t/tmax*(range(2)-range(1)) + range(1));
    y = 0.5*((exp(z)-1)/(exp(z)+1) + 1);
end


function [selected_index] = index_roulette_wheel_selection(F, k)
    fitness = F(1:k);
    weights = max(fitness) - fitness;
    weights = cumsum(weights/sum(weights));
    
    selected_index = roulette_wheel_selection(weights);
end

function [selected_index] = roulette_wheel_selection(weights)
    r = rand();
    selected_index = 1;
    for index=size(weights,1)
        if r <= weights(index)
            selected_index = index;
            break;
        end
    end
end

解释： - N: 搜索代理（种群个体）数量 - dim: 问题维度 - lb: 搜索空间下界向量 - ub: 搜索空间上界向量 - fobj: 目标函数（匿名函数） - maxEvals: 最大适应度函数评估次数

函数返回以下输出结果：

  • fbest: 找到的最小目标函数值
  • xbest: 对应最小目标函数值的最优解
  • Convergence_curve: 收敛曲线数据

算法参数定义 - 定义了AROA算法中涉及的多个常数参数，如c、fr1、fr2等，这些参数用于控制算法中不同操作的比例、强度等。

初始化 - 计算最大迭代次数tmax，初始化评价计数器evalCounter。 - 初始化用于记录全局最优解的fbest、xbest，以及用于记录收敛曲线的Convergence_curve。 - 初始化搜索空间边界向量Xmin、Xmax。 - 随机生成初始种群X，并对其进行适应度评估（通过evaluate_population函数），得到目标函数值F。根据F找到当前全局最优解fbest、xbest，以及保存当前种群状态的X_memory和F_memory。

主循环（tmax次迭代） - 计算种群间距离矩阵D（欧氏距离）。 - 根据当前迭代次数t计算参数m。 - 对每个搜索代理（i）执行以下操作： - 计算吸引力-排斥力操作： - 计算当前搜索代理与其他邻居间的吸引力-排斥力向量delta_ni。 - 合并所有邻居的贡献得到ni。 - 吸引至全局最优解操作： - 根据概率p1决定是否使用随机缩放的全局最优解。 - 计算吸引至全局最优解的向量bi。 - 局部搜索操作： - 根据概率p2决定采用哪种局部搜索策略（高斯扰动、邻域扰动或均匀扰动）。 - 根据所选策略计算局部搜索向量ri。 - 更新搜索代理位置：X(i,:) = X(i,:) + ni + bi + ri。 - 对更新后的种群X进行适应度评估（通过evaluate_population函数），并检查是否找到新的全局最优解。 - 执行记忆操作（通过memory_operator函数），更新种群X和目标函数值F，同时更新记忆中的种群状态X_memory和F_memory。 - 执行混沌扰动操作，对种群进行混沌移动，再次进行适应度评估并检查全局最优解。 - 再次执行记忆操作，更新相关变量。 - 记录当前迭代的全局最优解fbest至收敛曲线Convergence_curve。

辅助函数 - evaluate_population(X, fobj, ub, lb, evalCounter, maxEvals): 评估种群X在目标函数fobj下的适应度，确保解位于搜索空间内，并更新评价计数器。 - memory_operator(X, F, X_memory, F_memory): 根据记忆中的最优解和当前种群的最优解，选择保留更好的解。 - tanh(t, tmax, range): 计算双曲正切函数（tanh）的值，用于调整算法参数。 - index_roulette_wheel_selection(F, k): 使用轮盘赌选择法在前k个个体中按适应度比例选取一个个体的索引。 - roulette_wheel_selection(weights): 实现基本的轮盘赌选择算法，根据权重向量weights随机选择一个索引。

选一个函数做测试（2维0-10范围的最小值为:-18.5547） 

1
2
3
4
5
6
7
function z = F22(x)
    X = x(:, 1);
    Y = x(:, 2);    
    %z = X.^2 + Y.^2 + (25 * (sin(X).^2 + sin(Y).^2));
    z=X.*sin(4.*X)+ 1.1.*Y.*sin(2.*Y);
    y = z;
end 
1
2
3
4
5
6
7
8
9
10
11
12
13
...
Y =

  -18.5547

meanV =

  -18.5547

stdV =

   1.4862e-14
>>

测试函数： ## Rastrigin 

1
2
3
4
5
6
function z = F22(x)
% F22- Rastrigin  lower=[-5.12], upper=[5.12] gminimum=[0], dim=30
n = size(x, 2);
A = 10;
z = (A * n) + (sum(x .^2 - A * cos(2 * pi .* x), 2));
end
1
2
3
4
5
6
7
8
9
10
11
Y =

   1.0630e-11

meanV =

   24.0180

stdV =

   23.2159
Rastrigin函数是一个典型的非线性多峰函数，在搜索区域内存在许多极大值和极小值，导致寻找全局最小值比较困难，常用来测试寻优算法的性能。Rastrigin函数表达式和函数图像如下： alt text alt text

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
clc
clear
close all
[x,y]= meshgrid(-5:0.1:5);

z = x.^2+y.^2+(0.5.*i.*x).^2+(0.5.*i*y).^2+(0.5.*i.*x).^4+(0.5.*i.*y).^4;

xx=15*rand(30,1)-5;
yy=15*rand(30,1)-5;
zz = xx.^2+yy.^2+(0.5.*i.*xx).^2+(0.5.*i*yy).^2+(0.5.*i.*xx).^4+(0.5.*i.*yy).^4;

[nx,ny] = min(z);
[mx,my] = min(nx);
figure
subplot(1,2,1)
surf(x,y,z)
hold on
plot3(xx,yy,zz,'k*')%,'MarkerFaceColor','r'
shading interp
xlabel('x')
ylabel('y')
zlabel('f')
colormap(jet)
view([-107 49])
set(gca,'fontsize',12)
subplot(1,2,2)
surf(x,y,z)
hold on
plot3(x(ny(1),my),y(ny(1),my),mx,'ro','MarkerFaceColor','r','MarkerSize',10)%,'MarkerFaceColor','r'
shading interp
xlabel('x')
ylabel('y')
zlabel('f')
view([-107 49])
colormap(jet)
set(gca,'fontsize',12)

Schwefel’s Problem 1.2

Schwefel’s Problem 1.2

https://al-roomi.org/component/content/article?id=188:schwefel-s-function-no-1-2-double-sum-or-rotated-hyper-ellipsoid-function

https://www.oreilly.com/library/view/evolutionary-computation-with/9781848218079/16_appendix01.xhtml

https://blog.csdn.net/abc991835105/article/details/127936204

1
2
3
4
5
6
7
8
9
10
11
12
13
function z = F2(x)
% F22- Schwefel 1.2 lower=[-100], upper=[100], gminimum=[0], dim=30
z = 0;
v1 = 0;
dim=size(x,2);
for i = 1:dim
    v1=0;
    for j = 1: i
        v1 = v1 + x(:,j);
    end
    z = z + v1.^2;
end
end

alt text alt text alt text alt text alt text

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
% Schwefel's Problem # 1.2 (Double-Sum or Rotated Hyper-Ellipsoid Function)
% Range of initial points: -100 <= xj <= 100 , j=1,2,...,n
% Global maxima: (x1,x2,...,xn)=0
% f(X)=0
% Coded by: Ali R. Alroomi | Last Update: 07 June 2015 | www.al-roomi.org
    
clear
clc
warning off
    
x1min=-100;
x1max=100;
x2min=-100;
x2max=100;
R=1500; % steps resolution
x1=x1min:(x1max-x1min)/R:x1max;
x2=x2min:(x2max-x2min)/R:x2max;
 
for j=1:length(x1)
    
    % For 1-dimensional plotting
    f1(j)=x1(j)^2;
    
    % For 2-dimensional plotting
    for i=1:length(x2)
        fn(i)=x1(j)^2+(x1(j)+x2(i))^2;
    end
    
    fn_tot(j,:)=fn;
 
end

figure(1)
plot(x1,f1);set(gca,'FontSize',12);
xlabel('x','FontName','Times','FontSize',20,'FontAngle','italic');
ylabel('f(x)','FontName','Times','FontSize',20,'FontAngle','italic');
title('2D View','FontName','Times','FontSize',24,'FontWeight','bold');
 
figure(2)
meshc(x1,x2,fn_tot);colorbar;set(gca,'FontSize',12);
xlabel('x_2','FontName','Times','FontSize',20,'FontAngle','italic');
set(get(gca,'xlabel'),'rotation',25,'VerticalAlignment','bottom');
ylabel('x_1','FontName','Times','FontSize',20,'FontAngle','italic');
set(get(gca,'ylabel'),'rotation',-25,'VerticalAlignment','bottom');
zlabel('f(X)','FontName','Times','FontSize',20,'FontAngle','italic');
title('3D View','FontName','Times','FontSize',24,'FontWeight','bold');
 
figure(3)
mesh(x1,x2,fn_tot);view(0,90);colorbar;set(gca,'FontSize',12);
xlabel('x_2','FontName','Times','FontSize',20,'FontAngle','italic');
ylabel('x_1','FontName','Times','FontSize',20,'FontAngle','italic');
zlabel('f(X)','FontName','Times','FontSize',20,'FontAngle','italic');
title('X-Y Plane View','FontName','Times','FontSize',24,'FontWeight','bold');
 
figure(4)
mesh(x1,x2,fn_tot);view(90,0);colorbar;set(gca,'FontSize',12);
xlabel('x_2','FontName','Times','FontSize',20,'FontAngle','italic');
ylabel('x_1','FontName','Times','FontSize',20,'FontAngle','italic');
zlabel('f(X)','FontName','Times','FontSize',20,'FontAngle','italic');
title('X-Z Plane View','FontName','Times','FontSize',24,'FontWeight','bold');
 
figure(5)
mesh(x1,x2,fn_tot);view(0,0);colorbar;set(gca,'FontSize',12);
xlabel('x_2','FontName','Times','FontSize',20,'FontAngle','italic');
ylabel('x_1','FontName','Times','FontSize',20,'FontAngle','italic');
zlabel('f(X)','FontName','Times','FontSize',20,'FontAngle','italic');
title('Y-Z Plane View','FontName','Times','FontSize',24,'FontWeight','bold');
算法优化
吸引-排斥优化算法
http://jrhu0048.github.io/2024/04/26/xi-yin-pai-chi-you-hua-suan-fa/
作者
JR.HU
发布于
2024年4月26日
更新于
2024年4月27日
许可协议