RecursiveFrog 发表于 2022-1-7 13:14

【学界】进化算法简介(附代码)

文章作者：Poll
责任编辑：张帆（东北大学系统工程硕士生）
本篇文章转自【博客园】的优秀文章： 进化算法简介（原文链接网址：http://u6.gg/dFFyy）
欢迎原链接转发，转载请前往 @留德华叫兽 的主页获取信息，盗版必究。
敬请关注和扩散本专栏及同名公众号，会邀请全球知名学者陆续发布运筹学、人工智能中优化理论等相关干货、知乎Live及行业动态：『运筹OR帷幄』大数据人工智能时代的运筹学进化算法，也被成为是演化算法（evolutionary algorithms，简称EAs），它不是一个具体的算法，而是一个“算法簇”。进化算法的产生的灵感借鉴了大自然中生物的进化操作，它一般包括基因编码，种群初始化，交叉变异算子，经营保留机制等基本操作。与传统的基于微积分的方法和穷举方法等优化算法（具体介绍见博客 常见的几种最优化方法中的其他数学优化方法）相比，进化计算是一种成熟的具有高鲁棒性和广泛适用性的全局优化方法，具有自组织、自适应、自学习的特性，能够不受问题性质的限制，有效地处理传统优化算法难以解决的复杂问题（比如NP难优化问题）。
除了上述优点以外，进化算法还经常被用到多目标问题的优化求解中来，我们一般称这类进化算法为进化多目标优化算法（MOEAs）。目前进化计算的相关算法已经被广泛用于参数优化、工业调度、资源分配、复杂网络分析等领域。


1、（遗传算法（Genetic Algorithm，GA））

遗传算法（Genetic Algorithm，简称GA）是一种最基本的进化算法，它是模拟达尔文生物进化理论的一种优化模型，最早由J.Holland教授于1975年提出。遗传算法中种群分每个个体都是解空间上的一个可行解，通过模拟生物的进化过程，从而在解空间内搜索最优解。
遗传算法的基本操作可以用下图来描述：


个体的编码方式确定以后，针对上图操作的具体描述如下：
Step 1 种群初始化：根据问题特性设计合适的初始化操作（初始化操作应尽量简单，时间复杂度不易过高）对种群中的N个个体进行初始化操作；
Step 2 个体评价：根据优化的目标函数计算种群中个体的适应值（fitness value）；
Step 3 迭代设置：设置种群最大迭代次数gmax，并令当前迭代次数g=1；
Step 4 个体选择：设计合适的选择算子来对种群P(g)个体进行选择，被选择的个体将进入交配池中组成父代种群FP(g)，用于交叉变换以产生新的个体。选择策略要基于个体适应值来进行，假如要优化的问题为最小化问题，那么具有较小适应值的个体被选择的概率相应应该大一些。常用的选择策略有轮盘赌选择，锦标赛选择等；
Step 5 交叉算子：根据交叉概率pm（预先指定，一般为0.9）来判断父代个体是否需要进行交叉操作。交叉算子要根据被优化问题的特性来设计，它是整个遗传算法的核心，它被设计的好坏将直接决定整个算法性能的优劣；
Step 6 变异算子：根据变异概率pc（预先指定，一般为0.1）来判断父代个体是否需要进行变异操作。变异算子的主要作用是保持种群的多样性，防止种群陷入局部最优，所以其一般被设计为一种随机变换。
通过交叉变异操作以后父代种群FP(g)生成了新的子代种群P(g+1)，令种群迭代次数g=g+1，进行下一轮的迭代操作（跳转到Step 4），直至迭代次数达到最大的迭代次数。
为了更形象说明交叉操作的作用，我们以下图为例来深入理解一下交叉操作的作用：


通过交叉操作，原始的两个个体组合生成了两个新的个体组合，这就相当于在解空间进行搜索，每个个体都是解空间的一个可行解。
遗传算法matlab代码如下：
主函数main.m
% This example is used to explain the GA.
%max f(x1, x2) = 21.5 + x1·sin(4pi*x1) + x2·sin(20pi*x2)
%s.t.    -3.0 <=x1 <= 12.1
%         4.1 <=x2 <= 5.8
clc; clear all;
M = 40;N = 33;
generation = 1000;
it = 1;
pm = 0.05;%变异概率
pc = 0.8;%交叉概率
maxdata = -200;
pop = round(rand(M,N));%初始化矩阵，产生初始种群      
x1 = decode_x1(pop(:,1:18));%对x1 x2进行解码
x2 = decode_x2(pop(:,19:end));
fitness = 21.5 + x1.*sin(4*pi*x1) +x2.*sin(20*pi*x2);%适应度函数
while it < generation
    = seclection(fitness,pop);%轮盘赌选择
    = crossover(pc,B);%交叉
    = mutation(pm,newpop);%变异操作
    pop = B;
    x1 = decode_x1(pop(:,1:18));
    x2 = decode_x2(pop(:,19:end));
    fitness = 21.5 + x1.*sin(4*pi*x1) +x2.*sin(20*pi*x2);%计算适应度
    = max(fitness);%本代中的最优目标值
if fit_best >= maxdata
      maxdata = fit_best;
      verter = pop(index,:);
      x_1 = decode_x1(verter(1:18));
      x_2 = decode_x2(verter(19:end));
    end
    num(it) = maxdata;
    it = it + 1;
end
disp(sprintf('max_f=：%.6f',num(it-1)));%输出最优解
disp(sprintf('x_1=：%.5fx_2=：%.5f',x_1,x_2));%最优解对应的x1 x2的值
figure(1)   
plot(num,'k');%绘制图
变量x1的编码函数：decode_x1.m
function x1= encode(code)
= size(code);
sum = zeros(N);
for i=1:M
for j=1:N
      sum(i)=sum(i)+code(i,j)*2^(N-j);
    end
    x1(i) = -3.0 + sum(i)*(12.1-(-3.0))/(2^N - 1);
end
变量x2的编码函数：decode_x2.m
function x2=encode_x2(code)
= size(code);
sum = zeros(N);
for i=1:M
for j=1:N
      sum(i)=sum(i)+code(i,j)*2^(N-j);
    end
    x2(i) = 4.1 + sum(i)*(5.8 - 4.1)/(2^N - 1);
end
轮盘赌选择函数：seclection.m
%轮盘赌选择
function =seclection(fitness,A)
    =size(A);
    fit_sum = sum(fitness);
for i = 1:M
      probability(i) = fitness(i)/fit_sum;
    end
for i = 2:M
      probability(i) =probability(i) + probability(i-1);
    end
for j = 1:M
      p = rand();
      i = 1;
while p > probability(i)
            i = i+1;
      end
      B(j,:) = A(i,:);
    end
个体交叉算子：crossover.m
%交叉
function =crossover(pc,A)
newpop=A;
=size(A);
for i= 1:2:M-1
if rand < pc
            p = ceil(rand()*N);
for j= p:N
                ch = A(i,j);
                newpop(i,j) = A(i+1,j);
                newpop(i+1,j) = ch;
            end
      end
    end
个体变异算子：mutation.m
%变异
function =mutation(pm,A)
=size(A);
newpop=A;
W = rand(M,N)<pm;
for i=1:M
for j=1:N
if W(i,j) ==1
if A(i,j)~=1
                  newpop(i,j)= 1;
else
                  newpop(i,j) = 0;
               end
         end
       end
    end
2、文化基因算法（Memetic Algorithm，MA）

文化基因算法（Memetic Algorithm，简称MA），也被称为是“密母算法”，它是由Mpscato在1989年提出的。文化基因算法是一种基于种群的全局搜索和基于个体的局部启发式搜索的结合体，它的本质可以理解为：
Memetic = GA + Local Search
即memetic算法实质上为遗传算法加上一个局部搜索（Local Search）算子。局部搜索算子可以根据不同的策略进行设计，比如常用的爬山机制、模拟退火、贪婪机制、禁忌搜索等。

3、进化多目标优化算法（Multi-Objective Evolutionary Algorithm，MOEA）

对于一个优化问题而言，如果其只有一个目标方程，那么我们称之为单目标优化问题；而一旦方程个数达到两个或者两个以上，那么它被相应的称之为多目标优化问题（Multi-objective Optimization Problems，简称为MOPs）。
对于一个多目标优化问题而言，问题的最优解可能不止一个，而应该是一组，我们通常称这组最优解为相应多目标优化问题的一个非支配解集，或者称为是Pareto解集，其中的每一个解我们称之为Pareto解（Pareto是引自一个经济学的术语）。求解多目标优化问题的解法有很多，比如常见的目标规划方法，目标分解方法，目标化多为少方法（将多个目标表示为一个）等。进化算法在解决多目标问题上有着天然的优势，对于一个进化多目标优化算法而言，它可以对多个目标函数同时进行优化，而且输出一组非支配的Pareto解集，从而可以有效地求解多目标问题。
下图为一个具有两目标的多目标优化问题的Pareto解集示意图：


常用的进化多目标优化算法有Deb提出基于拥挤度距离度量的NSGA-II，张青富老师提出的基于分解思想的MOEA/D算法，以及西电公老师提出的基于免疫克隆机制的NNIA算法等。最近，Deb等人在NSGA-II的基础上又提出一种针对2~16个目标函数同时优化的高维问题的NSGA-III算法，对进化多目标优化感兴趣的博友可以深入的看看相关的参考文献，这方面的内容还是很有意思的。
NSGA-II算法的matlab代码参考如下：
主函数main.m
clear all
clc
tic
N=50;%the number of the population
numvariate=1;%the number of the variate of each individual
numfun=2;
minx=-10^3;
maxx=10^3;
iteration=1;
Pt0=initializationPt0(minx,maxx,N,numvariate,numfun);

while(iteration<500)
fprintf('-->iteration= %d\n',iteration);
=Genetic_Operators( Pt0,minx,maxx,numvariate,numfun);%%执行交叉变异后，子代个体的支配面号，拥挤度距离，这两项还没计算
[ Fi] = fast_nondominated_sort( Rt,numvariate,numfun);%Fi存储占优面上相应每个个体在种群Rt中的编号，即行号
for i=1:2*N
    numFi=Fi(i,1);
if (numFi~=0)
         =crowdiing_distance_assignment(Fi(i,:),Rt,numvariate,numfun);
for j=1:numFi
            Rt(Fi(i,j+1),numvariate+numfun+1)=i;
             Rt(Fi(i,j+1),numvariate+numfun+2)=Idistance(j);
         end
    end
end
Rt;

Pt1=zeros(1,N+1);
i=1;
while((Pt1(1)+Fi(i,1))<=N)
   Idistance1=crowdiing_distance_assignment(Fi(i,:),Rt,numvariate,numfun);
for j=1:Fi(i,1)
      Pt1(1)=Pt1(1)+1;
       Pt1(Pt1(1)+1)=Fi(i,j+1);
   end
   i=i+1;
end
if Pt1(1)<N
      Idistance=crowdiing_distance_assignment(Fi(i,:),Rt,numvariate,numfun);
   =sort(Idistance);
   NumFi=Fi(i,1);
for j=(Pt1(1)+2):(N+1)
      k=j-Pt1(1);
      Pt1(j)=Fi(i,rankIdis(NumFi-k+2)+1);
    end
end
b=Pt1(2:N+1);
Pt0=Rt(b,:);
iteration=iteration+1;
end
f1=Pt0(:,(numvariate+1));
f2=Pt0(:,(numvariate+2));
figure;
plot(f1,f2,'bo');
xlabel('f1');
ylabel('f2');
toc

目标函数：funfvalue.m
%%%input :the population Rt of size N*m,output:the value of the population
%%%Rt of size N*numf,numf if the number of the objective
function [ funvalueRt ] = funfvalue( Rt)
%% 测试数据SCH
funvalueRt(:,1)=Rt.^2;
funvalueRt(:,2)=(Rt-2).^2;
种群初始化函数：GeneratePt0.m
% %The interval is divided into N/10 equal parts,
function [ Pt0 ] = GeneratePt0(minx,maxx,N,numvariate,numfun)
Pt0=unifrnd(minx,maxx,N,numvariate);%种群数量为100
fPt0=zeros(N,numfun+2);
Pt0=cat(2,Pt0,fPt0);
end
基因操作（交叉和变异）：Genetic_Operators.m
% % input Pt ,size of N*m;output Rt,size of 2N*m
function [ Rt ] = Genetic_Operators( Pt,minx,maxx,numvariate,numfun)
%% 不交叉则变异，避免种群中数值相同的个体存在，这样可以保持种群的多样性
pc=0.9;
pm=0.1;
cetac=20;
cetam=20;
maxminx=maxx-minx;
N=size(Pt,1);
numvariate;
Qt=Pt;
crossPt=zeros(N,numvariate);
for i=1:numvariate
crossPt(:,i)=randperm(N);
end
detamax=10;%这个参数有待调整
for j=1:numvariate %%先交叉
for i=1:ceil(N/2)
      pccross=rand(1);
if(pccross<pc)
             ui=rand(1);
if ui<=0.5
               betau=(2*ui)^(1/(cetac+1));
else
if ui==1
                     ui=1-0.005;
               end
               betau=(1/(2*(1-ui)))^(1/(cetac+1));
             end
            Qt1=0.5*((1-betau)*Pt(crossPt(i,j),j)+(1+betau)*Pt(N-crossPt(i,j)+1,j));
            Qt2=0.5*((1+betau)*Pt(crossPt(i,j),j)+(1-betau)*Pt(N-crossPt(i,j)+1,j));
            Qt(2*i-1,j)=mod(Qt1-minx,maxminx)+minx;
            Qt(2*i,j)=mod(Qt2-minx,maxminx)+minx;
else
             pmu1=rand(1);
if pmu1<0.5
                deta1=(2*pmu1)^(1/(cetam+1))-1;
else deta1=1-(2*(1-pmu1))^(cetam+1);
             end
             pmu2=rand(1);
if pmu2<0.5
                deta2=(2*pmu2)^(1/(cetam+1))-1;
else deta2=1-(2*(1-pmu2))^(cetam+1);
             end
            Qt1=Pt(crossPt(i,j),j)+deta1*detamax;
            Qt2=Pt(N-crossPt(i,j)+1,j)+deta2*detamax;
            Qt(2*i-1,j)=mod(Qt1-minx,maxminx)+minx;
            Qt(2*i,j)=mod(Qt2-minx,maxminx)+minx;
          end
   end
end

funvalueQt=funfvalue(Qt(:,1:numvariate));
Qt(:,(numvariate+1):(numvariate+numfun))=funvalueQt;
Rt=cat(1,Pt,Qt);
end
快速的非支配排序函数：fast_nondominated_sort.m
% %%Fii存储占优面上相应每个个体在种群Rt中的编号，即行号
function = fast_nondominated_sort( Rt,numvariate,numfun)
N=size(Rt,1);
Sp=zeros(N,N+1);%Sp的第一列表示每个Spi所拥有的个数，即p所支配的个数，从第二列开始表示被p支配的个体
np=zeros(N,1);
Finum=zeros(N,N+1);
% Fi=zeros(N,N+1,numvariate);
prank=zeros(N,1);
funvalueRt=Rt(:,(numvariate+1):(numvariate+numfun));
for p=1:N
for q=1:N
if q~=p
%%% 这里还有另外一种情况，即Rt(p)不支配Rt(q)，也不被Rt(q)支配
if dominate(funvalueRt(p,:),funvalueRt(q,:))==1
                Sp(p,1)=Sp(p,1)+1;
                Sp(p,Sp(p,1)+1)=q;
            elseif dominate(funvalueRt(q,:),funvalueRt(p,:))==1
                  np(p)=np(p)+1;
            end
      end
    end

if np(p)==0
      prank(p)=1;
      Finum(1,1)=Finum(1,1)+1;
      Finum(1,Finum(1,1)+1)=p;
    end
end
Sp;
% np;
i=1;
while Finum(i,1)~=0
    Q=[];
%   fprintf('-->i= %d\n',i);
    Finum(i,:);
for ip=2:(Finum(i,1)+1)
      p=Finum(i,ip);
for iq=2:(Sp(p,1)+1)
            q=Sp(p,iq);
            np(q)=np(q)-1;
if np(q)==0
                prank(q)=i+1;
                Q=;
            end
      end
    end
    i=i+1;
    nun2Q=size(Q,2);
    Finum(i,1)=nun2Q;
for j=1:nun2Q
      Finum(i,j+1)=Q(j);
    end
end
Fii=Finum;
end
拥挤度距离计算函数：
%%% input:Fi of size 1*(2*N+1),output:Idistance of size 1*(FiRt(1,1))
function = crowdiing_distance_assignment(Fi,Rt,numvariate,numfun)
NFi=Fi(1,1);
Idistance=zeros(1,NFi);
funvalueI1=zeros(NFi,numfun);
if NFi<=2
    Idistance=Inf*ones(1,NFi);
else
    b=Fi(2:(Fi(1,1)+1));
    funvalueI1=Rt(b,(numvariate+1):(numvariate+numfun));
%   b=Fi(2:(Fi(1,1)+1));
%   for i=1:NFi
%         funvalueI1(i,:)=Rt(b(i),(numvariate+1):(numvariate+numfun));
%   end

%   fprintf('-->Idistancel= %d\n',l);
for i=1:numfun
      =sort(funvalueI1(:,i));
      Idistance(rankI1(1))=Inf;
      Idistance(rankI1(NFi))=Inf;
for j=2:(NFi-1)
            Idistance(rankI1(j))=Idistance(rankI1(j))+(valueI1(j+1)-valueI1(j-1))/(valueI1(NFi)-valueI1(1));
      end
    end
end
end

支配函数：dominate.m
%% 如果Rtp支配Rtq，返回值value=1，否则value=0；
function = dominate(Rtp,Rtq)
funvaluepq=Rtp-Rtq;
value=0;
if funvaluepq<0
    value=1;
elseif funvaluepq<=0
if(sum(funvaluepq)<0)
      value=1;
else value=0;
    end
end
end

参考文献

进化多目标优化算法研究，软件学报，2009.
Messy genetic algirithms: Motivation, analysis, and first results. Complex systems, 1989.
Genetic algorithms in search, optimization, and mechine learning. Addion wesley, 1989.
A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transaction on Evolutionary Computation, 2002.
<hr/>如果你是运筹学/人工智能硕博或在读，请在下图的公众号后台留言：“加微信群”。系统会自动辨认你的关键字，并提示您进一步的加群要求和步骤，邀请您进全球运筹或AI学者群（群内学界、业界大佬云集）。
同时我们有：【运筹学|优化爱好者】【供应链|物流】【人工智能】【数据科学|分析】千人QQ群，想入群的小伙伴可以关注下方公众号点击“加入社区”按钮，获得入群传送门。
学术界|工业界招聘、征稿等信息免费发布，请见下图：

APSchmidt 发表于 2022-1-7 13:19

有关于GA的书籍推荐吗

unityloverz 发表于 2022-1-7 13:28

关键词回复里貌似有几本书的，也欢迎加学术群交流

unityloverz 发表于 2022-1-7 13:36

请问编码函数中的sum=zeros(N)是不是错了，表示个体的应该是M吧

JamesB 发表于 2022-1-7 13:39

could be

zt3ff3n 发表于 2022-1-7 13:43

大佬，有文化基因算法的代码吗？
[害羞]

Zephus 发表于 2022-1-7 13:45

我不是做这块的，欢迎加学术群交流哈~

Mecanim 发表于 2022-1-7 13:46

谢谢

Mecanim 发表于 2022-1-7 13:47

能不能推荐下进化算法的书籍

RecursiveFrog 发表于 2022-1-7 13:51

欢迎进 @运筹OR帷幄 学术群交流
页: [1] 2
查看完整版本: 【学界】进化算法简介(附代码)