---

EP(Evolutionary Programming)

演化规划是由L.J.FogeI等在20世纪60年代提出的。当时演化规划的目标是通过模拟进化来获得智能行为。 他们将智能视为能够预测其所在环境的状态，并按照预定目标作出适当响应的能力。对环境的预测能力是智能行为的一个重要特征。

关键理念：

个体的变异是产生新个体的唯一方式

认为变异参数是演化的关键

强调用自适应方式不断优化变异参数

演化规划的基本框架

个体表示：
(x1,x2,…,xn,σ1,σ2,…,σn)，前面的x为一个个体n维，而后面的σ为此个体的各个维度变异强度。一个个体的强度可以相同

接下来介绍EP算法的主要部分——变异和生存选择：

EP算法的变异对于使用这个算法的人来说是非常友好的，因为它直接给出了公式。公式的具体内容就是对于一个个体的一维更新是在这个维度上直接加上一个以0为期望，对应变异强度为方差的正态分布随机数 而对于变异强度的更新则是对于一个个体的每一维的变异强度加上一个以0为期望，倍数于自身为方差的正态分布随机数，这个倍数可以是1。变异强度不能太小如果太小了，则需要将其拉回一个临界

当μ个父代个体的都产生了子代，那么就会存在μ+μ个个体，如何从这μ+μ个个体中选出下一代呢？这个就是生存选择索要解决的事情，我们可以采用(μ+μ)随机q竞争选择，也可以简单排序择优选择。 以下是(μ+μ)随机q竞争选择的步骤：

有了这些东西我们就可以编出相应的代码：

function [answer,best]=EPpro(fun)
global popnum ubound lbound vs dim cost times
best=inf;
answer=zeros(times,cost/popnum);
for i=1:times
    [X,D]=initpop();
    for j=1:(cost/popnum) 
        [X,D,good]=iteration(X,D);
        last=keepbest(good);
        answer(i,j)=last;
        if last<best
            best=last;
        end
    end
end


function [X,D]=initpop() %生成种群
X=unifrnd(lbound,ubound,popnum,dim);
D=unifrnd(0,ubound,popnum,1);
end

function [newX,newD,good]=iteration(X,D) %迭代
newX=zeros(popnum,dim);
newD=zeros(popnum,1);
for i1=1:popnum
    for j1=1:dim
        newX(i1,j1)=X(i1,j1)+D(i1)*randn();
    end
    newD(i1)=D(i1)+normrnd(0,sqrt(D(i1)));
    if newD(i1)<vs
        newD(i1)=vs;
    end
end
f=fit([X;newX]);
luck=selecsur(f);
good=f(luck);
x=[X;newX];
newX=x(luck,:);
d=[D;newD];
newD=d(luck);
end

function f=fit(X)
f=zeros(length(X),1);
for i2=1:length(X)
    f(i2)=fun(X(i2,:));
end
end

function luck=selecsur(X)
score=zeros(length(X),1);
for i3=1:length(X)
    select=randperm(length(X),round(length(X)*0.8)); %随机选取小组成员下标
    compare=X(select); %小组成员
    score(i3)=sum(compare<X(i3)); %求第i个个体在其小组得分
end
[~,b]=sort(score); %评选优秀个体
luck=b(1:length(b)/2); %存活者下标
end

function best=keepbest(X)
best=min(X);
end
end

测试branin函数:

global popnum ubound lbound vs dim cost times
popnum=50; %种群数量
ubound=15; %个体维度大小上限
lbound=-5; %个体维度大小下限
vs=0.3;  %变异强度临界值
dim=2;   %个体维度大小
cost=20000;  %计算代价
times=30;   %重复计算次数
fun=@(x)(x(2)-(5.1/(4*pi^2))*x(1)^2+5*x(1)/pi-6)^2+10*(1-1/(8*pi))*cos(x(1))+10;
[answer,best]=EPpro(fun) %best为EP找到的最优值，answer为每次计算得到的最好结果变化矩阵
plt(answer);
function plt(answer)  %绘制L曲线图
hisbest=reshape(answer',1,size(answer,1)*size(answer,2));
hisbest(1)=max(hisbest);
for k=2:length(hisbest)
    if hisbest(k)>hisbest(k-1)
        hisbest(k)=hisbest(k-1);
    end
end
for t=2:length(hisbest)
    hisbest(t)=(hisbest(t)+hisbest(t-1)*(t-1))/t;
end
plot(log(hisbest));
hold on
end

结果如下：