matlab練習程序（演化策略ES）

還是這本書上的內容，不過我看演化計算這一章是倒着看的，這里練習的算法正好和書中介紹的順序是相反的。

演化策略是最古老的的演化算法之一，和上一篇DE算法類似，都是基於種群的隨機演化產生最優解的算法。

算法步驟如下：

1.設定種群個體數和需要迭代的次數。
2.選擇父代中的個體按照公式z1=sqrt(-2*ln(u1))*sin(2*pi*u2)*m，z2=sqrt(-2*ln(u1))*cos(2*pi*u2)*m進行演化。

這里u1,u2都是隨機值，m是控制因子，演化次數越多m，m越小，父代通過與z1，z2相加得到后代。

3.計算后代的適應性。

4.選擇后代中最優的適應性作為全局最優適應性。

其實整個過程和DE非常類似。過程都是隨機變異，求適應性，再找最優。

我還試着將z1和z2橫設為1，竟也能得到非常好的解。

算法結果如下：

matlab代碼如下：

main.m

clear all;close all;clc;

[x y]=meshgrid(-100:100,-100:100);
sigma=50;
img = (1/(2*pi*sigma^2))*exp(-(x.^2+y.^2)/(2*sigma^2)); %目標函數，高斯函數
mesh(img);
hold on;
n=50;       %種群個體的數量
iter=100;   %迭代次數

%初始化種群，定義結構體
par=struct([]);
for i=1:n
    par(i).x=-100+200*rand();       %個體的x特征在[-100 100]隨機初始化
    par(i).y=-100+200*rand();       %個體的y特征在[-100 100]隨機初始化
    par(i).fit=compute_fit(par(i));       %個體在[x,y]處的適應度
end
par_best=par(1);    %初始化種群中最佳個體

for k=1:iter     %迭代次數
    plot3(par_best.x+100,par_best.y+100,par_best.fit,'g*'); %畫出最佳個體的位置，+100為相對偏移
    [par par_best]=select_and_recombin(par,par_best,n,k,iter);     %差異演化函數
end

select_and_recombin.m

function [next_par par_best]=select_and_recombin(par,par_best,n,k,iter)
    mul=(iter-k)/iter;  %限制進化因子，代數越高變異越小
    next_par=par;       %新種群
    for i=1:n
        
        %產生變異隨機數        
        u1=rand();  
        u2=rand();
        z1=sqrt(-2*log(u1))*sin(2*pi*u2)*mul;
        z2=sqrt(-2*log(u1))*cos(2*pi*u2)*mul;

        %變異
        next_par(i).x=par(i).x+z1;      
        next_par(i).y=par(i).y+z2;              
      
        %計算變異后個體的適應度
        next_par(i).fit=compute_fit(next_par(i));
        %如果新個體沒有變異前個體適應度高，新個體還原為舊個體
        if par(i).fit>next_par(i).fit
            next_par(i)=par(i);
        end
        %如果變異后適應度高於種群最高適應個體，則更新種群適應度最高個體
        if next_par(i).fit>par_best.fit
            par_best=next_par(i);
        end
    end    
end

compute_fit.m

function re=compute_fit(par)
    x=par.x;
    y=par.y;
    sigma=50;
    if x<-100 || x>100 || y<-100 || y>100 
        re=0;        %超出范圍適應度為0
    else            %否則適應度按目標函數求解
        re=(1/(2*pi*sigma^2))*exp(-(x.^2+y.^2)/(2*sigma^2)); 
    end
end