I. Introduction to artificial ecosystem algorithms

Artificial ecosystem based optimization (AEO) is a new optimization method to solve optimization problems. Driven by the flow of energy in the earth’s ecosystems, SDO mimics three unique behaviors of living things, including production, consumption and decomposition. AEO algorithm is a new meta-heuristic optimization algorithm proposed by Zhao et al. [27] in 2019 by simulating energy flow in the earth’s ecosystem. The algorithm simulates production, consumption and decomposition behaviors in the ecosystem through production operators, consumption operators and decomposition operators to achieve the purpose of solving optimization problems. The production operator is designed to strengthen the balance between AEO algorithm exploration and development. The exploration ability of AEO algorithm improved by consumption operator; The decomposition operator aims to improve the performance of AEO algorithm development. Compared with the traditional swarm intelligence algorithm,AEO algorithm is not only simple to implement, but also does not need to adjust any other parameters except the population size and the maximum number of iterations, and has better optimization accuracy and global search ability.

1 AEO algorithm follows the following three criteria: ① The ecosystem as a population includes three kinds of organisms: producer, consumer and decomposer, and there is only one individual as producer and decomposer respectively in the population, and the other individuals as consumers; (2) Each individual has the same probability of being selected as a carnivore, herbivore or omnivore; ③ The energy level of each individual in the population was evaluated by fitness value, which was sorted in descending order. The higher the fitness value was, the higher the energy level of the minimization problem was.

2 Mathematical description of AEO algorithm** A. In ** ecosystems, producers can use CO2, water and sunlight, as well as nutrients from decomposers, to produce food energy. In AEO algorithm, the producer (worst individual) in the population updates by searching upper and lower limits of space and decomposer (optimal individual), and the updated individual guides other individuals in the population to search different regions. The mathematical model to simulate producer behavior is as follows:Where, x1 is the individual spatial position of the producer; Xn is the optimal individual spatial position in the current population. N is population size; T is the maximum number of iterations; T is the current iteration number; Xrand is the randomly generated individual spatial position in the search space; U and L are the upper and lower limits of space respectively. R and r1 are random numbers between [0,1].

** B. ** After the producer provides food energy, each consumer can randomly select the consumer or producer with a lower energy level or both to obtain food energy. If the consumer is randomly selected as an herbivore, it feeds only on the producer; If a consumer is randomly selected as a carnivore, it can only randomly select consumers with higher energy levels to eat; If a consumer is randomly selected as an omnivore, it can choose to feed on both consumers and producers with higher energy levels. The mathematical models simulating the consumption behavior of herbivores, carnivores and omnivores are respectivelyWhere, Xi is the spatial location of the ith individual consumer; C is the consumption factor with Levy flight characteristics; N(0,1) is a probability density function with normal distribution, mean value of 0 and standard deviation of 1. Xj is the consumer and producer with higher energy level; R2 is a random number in the range of [0,1].

** C. ** Decomposition is a very important process in terms of ecosystem function, providing essential nutrients for producers. In order to improve the development performance of the algorithm,AEO algorithm allows the next position of each individual to propagate around the best individual (decomposer), and updates the spatial position of the ith consumer in the population by adjusting decomposition factor D and weight coefficients E and H. The mathematical model simulating decomposition behavior is 2 combination growth modelA. Weibull model. Weibull model was first proposed by Swedish engineer Waloddi Weibull in 1951 and gradually developed. At present, Weibull model has been applied in fields such as oil and gas field output, foundation settlement, debris flow warning, drug dissolution curve evaluation and so on. Weibull model has various expressions. This paper uses the function model shown in Equation (6) to predict water demand.Where, W ‘1 is the predicted value of water demand in Weibull model; Z ˆ ‘o is the BTH normalized value of the water requirement prediction influence factor, that is, the O-dimensional input vector. O (o=1,2… M, there are 7 influencing factors of water demand prediction in this paper,m=7) is the parameter to be optimized.

B. Richards model. Richards model is a nonlinear regression equation describing biological growth, containing four parameters, which has been applied in carbon emission prediction, population prediction, land surface subsidence and other fields. Richards model has various expressions. This paper uses the function model shown in Equation (7) to predict water demand.Where, W ‘2 is the predicted water demand of Richards model; A, B, CO and D are parameters to be optimized.

C. Usher model. Usher model was first proposed by Usher, an American scholar, in 1980 to describe the change of growth information with time. It has been applied in surrounding rock deformation prediction, oilfield development, settlement prediction and other fields. Usher model is expressed in various forms. In this paper, the functional model shown in Equation (8) is used to predict water demand.Where, W ‘3 is the predicted value of water demand of Usher model; A, B, Co and D are parameters to be optimized. AEO algorithm was used to optimize parameters and weight coefficients of Weibull-Richards-Usher, Weibull-Richards, Weibull-Usher and Richards-Usher models, and four combined growth models to be optimized were obtained:

Two, some source code

% -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - % -------------------------------------------------------------------------- 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
 

Copy the code

3. Operation results

Matlab version and references

1 matlab version 2014A

[1] Yang Baoyang, YU Jizhou, Yang Shan. Intelligent Optimization Algorithm and Its MATLAB Example (2nd Edition) [M]. Publishing House of Electronics Industry, 2016. [2] ZHANG Yan, WU Shuigen. MATLAB Optimization Algorithm source code [M]. Tsinghua University Press, 2017.