Circle-Inspired Optimization Algorithm (CIOA) 圆启发优化算法

Circle-Inspired Optimization Algorithm

参考文献：CIOA: Circle-Inspired Optimization Algorithm, an algorithm for engineering optimization
该博客主要分享CIOA算法的代码，并带有中文注释！
分享不易，喜欢的请大家点赞加收藏，谢谢！
(1) CIOA主程序

% CIOA算法
% 参考文献：CIOA: Circle-Inspired Optimization Algorithm, an algorithm for engineering optimization
% Code by Luzhenhui Yangtze University
clc
clear 
format long
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输入参数
Nvar = 6;                   % 变量的数目
Nag = 250;                  % 种群个体数目
Nit = 400;                  % 迭代次数
Ub = [1 1 1 1 16 16];       % 设计变量的上界
Lb = [0 0 0 0 1e-5 1e-5];   % 设计变量的下界
ThetaAngle = 17;            % Theta 角度 (推荐值:13,17或19)
GlobIt = 0.85;              % 全局搜索占算法总的迭代次数的比例（推荐取值区间[0.75 0.95]）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 算法的参数初始化
Angle = (ThetaAngle*pi/180);
Convergence = zeros(1,Nit);             % 用于记录每一次迭代的最优值
BestVar = zeros(Nit,Nvar);              % 用于记录每一次迭代的最优解（及最优个体）
% x = ones(1,Nvar);
Bestobj = zeros(Nit,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 半径计算
raux = 1*(sqrt(max(Ub)-min(Lb)))/Nag;   % 半径辅助变量Cr 如公式（2）
r = zeros(1,Nag);                       % 初始化半径向量
h = Nag;
while h >= 1                            % 计算半径向量 如公式（1）
    r(h) = raux*(h)^2/Nag;              
    h = h-1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 算法种群个体初始化
i0 = 1;
Var0 = zeros(Nag,Nvar);
obj0 = zeros(Nag,1);
while i0 <= Nag
    Var0(i0,:) = Lb+(Ub-Lb).*rand(1,Nvar); 
    i0 = i0+1;
end
% 计算初始种群个体的适应度值
Sol0 = zeros(1,Nag);
i1 = 1;
while i1 <= Nag
    x = Var0(i1,:);                     
    [sol] = Objective_Function_CIOA(x); 
    Sol0(i1) = sol;                      
    i1 = i1+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 全局搜索
It = 1;
while It <= floor(Nit*GlobIt)
    [class,ord] = sort(Sol0);           % 种群个体按其适应度值进行排序（由小到大）
    i2 = 1;
    while i2 <= Nag
        i3 = 1;
        while i3 <= Nvar                % i3（某个种群个体的第i3个值），根据其奇偶选择不同的公式进行更新 如文章公式（3）和（4）
            if mod(i3,2)==0
                Var0(i2,i3) = (Var0(ord(i2),i3))-((0.0+(1-0.0)*rand())*r(i2)*sin((It*Angle)-Angle))+((0.0+(1-0.0)*rand())*r(i2)*sin(It*Angle));
            else
                Var0(i2,i3) = (Var0(ord(i2),i3))-((0.0+(1-0.0)*rand())*r(i2)*cos((It*Angle)-Angle))+((0.0+(1-0.0)*rand())*r(i2)*cos(It*Angle));
            end
            i3 = i3+1;
        end
        i2 = i2+1;
    end
    % 边界处理
    i4 = 1;
    while i4 <= Nag
        i5 = 1;
        while i5 <= Nvar 
            if Var0(1,i5) > Ub(i5)
               Var0(1,i5) = Ub(i5);
            elseif Var0(1,i5) < Lb(i5)
                   Var0(1,i5) = Lb(i5);                 
            end
            if Var0(i4,i5) > Ub(i5)
               Var0(i4,i5) = Var0(1,i5);
            elseif Var0(i4,i5) < Lb(i5)
                   Var0(i4,i5) = Var0(1,i5);                
            end
            i5 = i5+1;
        end
        i4 = i4+1;
    end
    % 计算更新后种群个体的适应度值
    i6 = 1;
    while i6 <= Nag
        x = Var0(i6,:);                          
        [sol, obj] = Objective_Function_CIOA(x); 
        Sol0(i6) = sol;                          
        obj0(i6) = obj;
        i6 = i6+1;
    end
    % 确定全局最优值及其相对应的最优个体（这里针对最小化问题，可根据自己的需求进行修改）
    [Best,Position] = min(Sol0);        
    BestSolution = Best;                
    BestVar(It,:) = Var0(Position,:);   
    Bestobj(It) = obj0(Position);
    Convergence(It) = BestSolution;
    % 经过一定迭代次数后，若r超过360，那么就需要下面的公式对r进行更新，以加快算法的收敛，如文中的公式（5）
    if mod(It,floor(360/ThetaAngle))==0
        r = r.*0.99;                    
    end
    It = It+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 局部搜索
Lb1 = BestVar(It-1,:)-(ones(1,Nvar).*((Ub-Lb)/10000)); % 公式（6）
Ub1 = BestVar(It-1,:)+(ones(1,Nvar).*((Ub-Lb)/10000)); % 公式（7）
i7 = 1;
Var1 = zeros(Nag,Nvar);
while i7 <= Nag
    Var1(i7,:) = BestVar(It-1,:);    
    i7 = i7+1;
end
x = Var1(1,:);
[sol] = Objective_Function_CIOA(x);  
Sol1 = zeros(1, Nag);
i8 = 1;
while i8 <= Nag
    Sol1(i8) = sol;
    i8 = i8+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 局部搜索的主要步骤
while It <= Nit
    [class,ord] = sort(Sol1);           
    i9 = 1;
    while i9 <= Nag
        i10 = 1;
        while i10 <= Nvar
            if mod(i10,2)==0
                Var1(i9,i10) = (Var1(ord(i9),i10))-((0.0+(1-0.0)*rand())*r(i9)*sin((It*Angle)-Angle))+((0.0+(1-0.0)*rand())*r(i9)*sin(It*Angle));
            else
                Var1(i9,i10) = (Var1(ord(i9),i10))-((0.0+(1-0.0)*rand())*r(i9)*cos((It*Angle)-Angle))+((0.0+(1-0.0)*rand())*r(i9)*cos(It*Angle));
            end
            i10 = i10+1;
        end
        i9 = i9+1;
    end
    i11 = 1;
    while i11 <= Nag
        i12 = 1;
        while i12 <= Nvar
            if Var1(1,i12) > Ub(i12)
               Var1(1,i12) = Ub(i12);
            elseif Var1(1,i12) < Lb(i12)
                   Var1(1,i12) = Lb(i12);
            end
            if Var1(i11,i12) > Ub1
               Var1(i11,i12) = Var1(1,i12);
            elseif Var1(i11,i12) > Ub(i12)
                   Var1(i11,i12) = Ub(i12);
            end
            if Var1(i11,i12) < Lb1
               Var1(i11,i12) = Var1(1,i12);
            elseif Var1(i11,i12) < Lb(i12)
                   Var1(i11,i12) = Lb(i12);
            end            
            i12 = i12+1;
        end        
        i11 = i11+1;
    end
    % 计算种群个体的适应度值并找出最优值及其所对应的最优个体
    i13 = 1;
    while i13 <= Nag
        x = Var1(i13,:);
        [sol, obj] = Objective_Function_CIOA(x); 
        Sol1(i13) = sol;
        obj0(i13) = obj;
        i13 = i13+1;
    end
    [Best,Position] = min(Sol1);
    BestSolution = Best;
    BestVar(It,:) = Var1(Position,:);
    Bestobj(It) = obj0(Position);
    Convergence(It) = BestSolution;

    if mod(It,floor(360/ThetaAngle))==0
        r = r.*0.99;
    end
    It = It+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 结果统计：全局最优值及其相对应的最优个体，绘制收敛曲线
ConvergenceCurve = zeros(1,Nit);
ConvergencePlot = zeros(1,Nit);
ConvergenceCurve(1) = Convergence(1);
ConvergencePlot(1) = Bestobj(1);
i14 = 2;
while i14 <= Nit
    if Convergence(i14) < ConvergenceCurve(i14-1)
        ConvergenceCurve(i14) = Convergence(i14);
        ConvergencePlot(i14) = Bestobj(i14);
    else
        ConvergenceCurve(i14) = ConvergenceCurve(i14-1);
        ConvergencePlot(i14) = ConvergencePlot(i14-1);
    end
    i14 = i14+1;
end
plot(ConvergencePlot','b-','LineWidth',2)  
xlabel('\fontsize{11}\bf Iteration Number')
ylabel('\fontsize{11}\bf Objective Function Value')
grid on
[BestSol,PositionSol] = min(Convergence);
BestSolution=ConvergencePlot(Nit);        
BestVariables=BestVar(PositionSol,:);    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 约束条件处理
x=BestVariables;
[sol,obj,rg,rh,nug,nuh]=Objective_Function_CIOA(x);
% 不等式约束条件
Violg = 0;
for icong=1:nug                        
    Violg = Violg+max(rg(icong),0);
end
% 等式约束条件
Violh = 0;
for iconh=1:nuh                      
    Violh = Violh+(max(abs(rh(iconh))-0.0001,0));
end
% 约束汇总
if (nug+nuh)==0
    Violation = (Violg+Violh)/1;        
else
    Violation = (Violg+Violh)/(nug+nuh);
end

(2) 测试函数

% Code by Luzhenhui Yangtze University
function [sol,obj,rg,rh,nug,nuh]=Objective_Function_CIOA(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 目标方程
k1 = 0.09755988;
k2 = 0.99*k1;
k3 = 0.0391908;
k4 = 0.9*k3;
obj = -x(4);       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 约束
nug = 1;           % 不等式约束数目
nuh = 4;           % 等式约束数目
rg = zeros(nug,1);
rh = zeros(nuh,1);
% 不等式约束：
rg(1) = x(5)^0.5+x(6)^0.5-4;
% 等式约束：
rh(1) = x(1)+k1*x(2)*x(5)-1;
rh(2) = x(2)-x(1)+k2*x(2)*x(6);
rh(3) = x(3)+x(1)+k3*x(3)*x(5)-1;
rh(4) = x(4)-x(3)+x(2)-x(1)+k4*x(4)*x(6);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 约束条件处理，惩罚函数
Weight = 10^5;     % 惩罚系数
% 惩罚-不等式约束
Pen = 0;
for icong=1:nug
    Pen = Pen+(max(rg(icong),0));              
end
% 惩罚-等式约束
for iconh = 1:nuh
    Pen = Pen+(max(abs(rh(iconh))-0.0001,0));  
end
Pen = Weight*Pen;
% 惩罚目标函数
sol = obj + Pen;   
end