【優化演算法】人工生態系統優化演算法（AEO）【含Matlab原始碼 023期】
一、人工生態系統演算法簡介
基於人工生態系統的優化（AEO）是一種解決優化問題的新型優化方法。 SDO受地球生態系統中能量流的驅動，該演算法模仿了生物的三種獨特行為，包括生產，消耗和分解。
AEO演算法是Zhao等[27]於2019年通過模擬地球生態系統中能量流動而提出一種新型元啟發式優化演算法,該演算法通過生產運算元、消費運算元和分解運算元對生態系統中的生產、消費和分解行為進行模擬來達到求解優化問題的目的。生產運算元旨在加強AEO演算法勘探和開發之間的平衡能力;消費運算元用於改進AEO演算法的探索能力;分解運算元旨在提升AEO演算法的開發效能。與傳統群智慧演算法相比,AEO演算法不但實現簡單,除群體規模和最大迭代次數外,無需調整其他任何引數,且具有較好的尋優精度和全域性搜尋能力。
1 AEO演算法遵行以下3個準則
①生態系統作為種群包括3種生物:生產者、消費者和分解者,且種群中分別只有一個個體作為生產者和分解者,其他個體作為消費者;
②每個個體都具有相同的概率被選擇為食肉動物、食草動物或雜食動物;
③群體中每個個體的能量水平通過適應度值進行評價,適應度值按降序排序,適應度值越大表示最小化問題的能量水平越高。
2 AEO演算法數學描述
a. 生產者。生態系統中,生產者可以利用CO2、水和陽光以及分解者提供的營養來產生食物能量。在AEO演算法中,種群中的生產者(最差個體)通過搜尋空間上下限和分解者(最優個體)進行更新,更新後的個體將引導種群中的其他個體搜尋不同的區域。
模擬生產者行為的數學模型如下:
式中:x1為生產者個體空間位置;xn為當前群體中最佳個體空間位置;n為種群規模;T為最大迭代次數;t為當前迭代次數;xrand為搜尋空間中隨機生成的個體空間位置;U、L分別為空間上、下限;r、r1為[0,1]之間的隨機數。
b. 消費者。生產者提供食物能量後,每個消費者均可隨機選擇能量水平較低的消費者或生產者或兩者兼有獲得食物能量。如果消費者被隨機選擇為食草動物,它只以生產者為食;如果消費者被隨機選擇為食肉動物,它只能隨機選擇能量水平較高的消費者為食;如果消費者被隨機選擇為雜食動物,它可以同時選擇能量水平較高的消費者和生產者為食。模擬食草動物、食肉動物、雜食動物消費行為的數學模型分別為
式中:xi為第i個消費者個體空間位置;C為具有levy飛行特性的消費因子;N(0,1)為呈正態分佈、均值為0、標準差為1的概率密度函式;xj為具有較高能量水平的消費者和生產者;r2為[0,1]範圍內的隨機數。
c. 分解者。就生態系統功能而言,分解是一個非常重要的過程,它為生產者提供必要的養分。為提高演算法的開發效能,AEO演算法允許每個個體的下一個位置圍繞最佳個體(分解者)傳播,並通過調節分解因子D和權重係數e、h來更新群體中第i個消費者的空間位置。其模擬分解行為的數學模型為
2 組合生長模型
a. Weibull模型。Weibull模型最早由瑞典工程師Waloddi Weibull於1951年提出並逐漸發展而來,目前已在油氣田產量、地基沉降、泥石流預警、藥物溶出曲線評價等領域得到應用。Weibull模型表述形式多樣,本文利用式(6)所示函式模型進行需水預測。
式中:W′1為Weibull模型需水預測值;zˆ′o為第b個需水預測影響因子歸一化值,即第o維輸入向量;α、γ、βo(o=1,2,…,m,本文需水預測影響因子有7個,m=7)為待優化引數。
b. Richards模型。Richards模型是描述生物生長的非線性迴歸方程,含有4個引數,目前已在碳排放量預測、人口預測、地表沉降等領域得到應用。Richards模型表述形式多樣,本文利用式(7)所示函式模型進行需水預測。
式中:W′2為Richards模型需水預測值;a、b、co、d為待優化引數。
c. Usher模型。Usher模型最早由美國學者Usher於1980年提出用於描述增長資訊隨時間變化的數學模型,目前已在圍巖變形預測、油田開發、沉降預測等領域得到應用。Usher模型表述形式多樣,本文利用式(8)所示函式模型進行需水預測。
式中:W′3為Usher模型需水預測值;A、B、Co、D為待優化引數。
利用AEO演算法對Weibull-Richards-Usher、Weibull-Richards、Weibull-Usher、Richards-Usher模型引數和權重係數進行優化,得到待優化的4種組合生長模型:
二、部分原始碼
%--------------------------------------------------------------------------
% --------------------------------------------------------------------------
clc;
clear;
MaxIteration=500;
PopSize=50;
FunIndex=1;
[BestX,BestF,HisBestF]=AEO(FunIndex,MaxIteration,PopSize);
display(['F_index=', num2str(FunIndex)]);
display(['The best fitness is: ', num2str(BestF)]);
%display(['The best solution is: ', num2str(BestX)]);
if BestF>0
semilogy(HisBestF,'r','LineWidth',2);
else
plot(HisBestF,'r','LineWidth',2);
end
xlabel('Iterations');
ylabel('Fitness');
title(['F',num2str(FunIndex)]);
%--------------------------------------------------------------------------
% --------------------------------------------------------------------------
% Artificial ecosystem-based optimization (AEO)
function [BestX,BestF,HisBestFit]=AEO(F_index,MaxIt,nPop)
% FunIndex: Index of function.
% MaxIt: The maximum number of iterations.
% PopSize: The size of population.
% PopPos: The position of population.
% PopFit: The fitness of population.
% Dim: The dimensionality of prloblem.
% C: The consumption factor.
% D: The decomposition factor.
% BestX: The best solution found so far.
% BestF: The best fitness corresponding to BestX.
% HisBestFit: History best fitness over iterations.
% Low: The low bound of search space.
% Up: The up bound of search space.
[Low,Up,Dim]=FunRange(F_index);
for i=1:nPop
PopPos(i,:)=rand(1,Dim).*(Up-Low)+Low;
PopFit(i)=BenFunctions(PopPos(i,:),F_index,Dim);
end
BestF=inf;
BestX=[];
[NFit, indF]=sort(PopFit,'descend');
PopPos=PopPos(indF,:);
PopFit=PopFit(indF);
BestF=PopFit(end);
BestX=PopPos(end,:);
HisBestFit=zeros(MaxIt,1);
Matr=[1,Dim];
for It=1:MaxIt
r1=rand;
a=(1-It/MaxIt)*r1;
xrand=rand(1,Dim).*(Up-Low)+Low;
newPopPos(1,:)=(1-a)*PopPos(nPop,:)+a*xrand; %equation (1)
u=randn(1,Dim);
v=randn(1,Dim);
C=1/2*u./abs(v); %equation (4)
newPopPos(2,:)=PopPos(2,:)+C.*(PopPos(2,:)-newPopPos(1,:)); %equation (6)
for i=3:nPop
u=randn(1,Dim);
v=randn(1,Dim);
C=1/2*u./abs(v);
r=rand;
if r<1/3
newPopPos(i,:)=PopPos(i,:)+C.*(PopPos(i,:)-newPopPos(1,:)); %equation (6)
else
if 1/3<r<2/3
newPopPos(i,:)=PopPos(i,:)+C.*(PopPos(i,:)- PopPos(randi([2 i-1]),:)); %equation (7)
else
r2=rand;
newPopPos(i,:)=PopPos(i,:)+C.*(r2.*(PopPos(i,:)- newPopPos(1,:))+(1-r2).*(PopPos(i,:)-PopPos(randi([2 i-1]),:))); %equation (8)
end
end
end
for i=1:nPop
newPopPos(i,:)=SpaceBound(newPopPos(i,:),Up,Low);
newPopFit(i)=BenFunctions(newPopPos(i,:),F_index,Dim);
if newPopFit(i)<PopFit(i)
PopFit(i)=newPopFit(i);
PopPos(i,:)=newPopPos(i,:);
end
end
[BestOne indOne]=min(PopFit);
for i=1:nPop
r3=rand; Ind=round(rand)+1;
newPopPos(i,:)=PopPos(indOne,:)+3*randn(1,Matr(Ind)).*((r3*randi([1 2])-1)*PopPos(indOne,:)-(2*r3-1)*PopPos(i,:)); %equation (9)
end
for i=1:nPop
newPopPos(i,:)=SpaceBound(newPopPos(i,:),Up,Low);
newPopFit(i)=BenFunctions(newPopPos(i,:),F_index,Dim);
if newPopFit(i)<PopFit(i)
PopPos(i,:)=newPopPos(i,:);
PopFit(i)=newPopFit(i);
end
end
[NFit,indF]=sort(PopFit,'descend');
PopPos=PopPos(indF,:);
PopFit=PopFit(indF);
for i=1:nPop
if PopFit(end)<BestF
BestF=PopFit(end);
BestX=PopPos(end,:);
end
end
HisBestFit(It)=BestF;
end
function Fit=BenFunctions(X,FunIndex,Dim)
switch FunIndex
%%%%%%%%%%%%%%%%%%%%%%%%%%unimodal function%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Sphere
case 1
Fit=sum(X.^2);
%Schwefel 2.22
case 2
Fit=sum(abs(X))+prod(abs(X));
%Schwefel 1.2
case 3
Fit=0;
for i=1:Dim
Fit=Fit+sum(X(1:i))^2;
end
%Schwefel 2.21
case 4
Fit=max(abs(X));
%Rosenbrock
case 5
Fit=sum(100*(X(2:Dim)-(X(1:Dim-1).^2)).^2+(X(1:Dim-1)-1).^2);
%Step
case 6
Fit=sum(floor((X+.5)).^2);
%Quartic
case 7
Fit=sum([1:Dim].*(X.^4))+rand;
%%%%%%%%%%%%%%%%%%%%%%%%%%multimodal function%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Schwefel
case 8
Fit=sum(-X.*sin(sqrt(abs(X))));
%Rastrigin
case 9
Fit=sum(X.^2-10*cos(2*pi.*X))+10*Dim;
%Ackley
case 10
Fit=-20*exp(-.2*sqrt(sum(X.^2)/Dim))-exp(sum(cos(2*pi.*X))/Dim)+20+exp(1);
%Griewank
case 11
Fit=sum(X.^2)/4000-prod(cos(X./sqrt([1:Dim])))+1;
%Penalized
case 12
a=10;k=100;m=4;
Dim=length(X);
Fit=(pi/Dim)*(10*((sin(pi*(1+(X(1)+1)/4)))^2)+sum((((X(1:Dim-1)+1)./4).^2).*...
(1+10.*((sin(pi.*(1+(X(2:Dim)+1)./4)))).^2))+((X(Dim)+1)/4)^2)+sum(k.*...
((X-a).^m).*(X>a)+k.*((-X-a).^m).*(X<(-a)));
%Penalized2
case 13
a=10;k=100;m=4;
Dim=length(X);
Fit=.1*((sin(3*pi*X(1)))^2+sum((X(1:Dim-1)-1).^2.*(1+(sin(3.*pi.*X(2:Dim))).^2))+...
((X(Dim)-1)^2)*(1+(sin(2*pi*X(Dim)))^2))+sum(k.*...
((X-a).^m).*(X>a)+k.*((-X-a).^m).*(X<(-a)));
%%%%%%%%%%%%%%%%%%%%%%%%%%fixed-dimensionalmultimodalfunction%%%%%%%%%%%%%%
%Foxholes
case 14
a=[-32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32;,...
-32 -32 -32 -32 -32 -16 -16 -16 -16 -16 0 0 0 0 0 16 16 16 16 16 32 32 32 32 32];
for j=1:25
b(j)=sum((X'-a(:,j)).^6);
end
Fit=(1/500+sum(1./([1:25]+b))).^(-1);
%Kowalik
case 15
a=[.1957 .1947 .1735 .16 .0844 .0627 .0456 .0342 .0323 .0235 .0246];
b=[.25 .5 1 2 4 6 8 10 12 14 16];b=1./b;
Fit=sum((a-((X(1).*(b.^2+X(2).*b))./(b.^2+X(3).*b+X(4)))).^2);
%Six Hump Camel
case 16
Fit=4*(X(1)^2)-2.1*(X(1)^4)+(X(1)^6)/3+X(1)*X(2)-4*(X(2)^2)+4*(X(2)^4);
%Branin
case 17
Fit=(X(2)-(X(1)^2)*5.1/(4*(pi^2))+5/pi*X(1)-6)^2+10*(1-1/(8*pi))*cos(X(1))+10;
%GoldStein-Price
case 18
Fit=(1+(X(1)+X(2)+1)^2*(19-14*X(1)+3*(X(1)^2)-14*X(2)+6*X(1)*X(2)+3*X(2)^2))*...
(30+(2*X(1)-3*X(2))^2*(18-32*X(1)+12*(X(1)^2)+48*X(2)-36*X(1)*X(2)+27*(X(2)^2)));
%Hartman 3
case 19
a=[3 10 30;.1 10 35;3 10 30;.1 10 35];c=[1 1.2 3 3.2];
p=[.3689 .117 .2673;.4699 .4387 .747;.1091 .8732 .5547;.03815 .5743 .8828];
Fit=0;
for i=1:4
end
end
三、執行結果
四、matlab版本及參考文獻
1 matlab版本
2014a
2 參考文獻
[1] 包子陽,餘繼周,楊杉.智慧優化演算法及其MATLAB例項（第2版）[M].電子工業出版社，2016.
[2]張巖,吳水根.MATLAB優化演算法原始碼[M].清華大學出版社，2017.