matlab練習程序（差異演化DE）

這兩天在看M.Tim Jones的《人工智能》，書中不只介紹原理，而且都有相應的c代碼實現。

雖然代碼不完全，不過缺少的部分完全可以自己補完。

差異演化和昨天實現的PSO很類似，都屬於優化算法。

算法步驟：

1.設定種群個體個數和需要迭代的次數。當然也可以設定條件，然后判斷是否停止迭代。

2.定義交叉概率CR，個體有一定概率進行變異，如果變異則進行第3步，如果不變異則下一代個體和當前個體一樣。

3.在種群中隨機選出三個互不相同的個體進行變異，變異公式如下（有博士論文總結了一大堆變異公式）：

X_i,g+1=X_r1,g+F(X_r2,g-X_r3,g)

其中X是種群中個體特征向量，這里就是x和y坐標。g是種群的代數。i代表當前個體，r1,r2,r3是和i不同，並且也互不相同的個體。

4.計算變異后個體的適應度，如果變異后適應度不如變異前，那么將變異后個體重新恢復為變異前個體。

5.比較變異后個體適應度和種群最優個體適應度，將適應度高的個體賦給種群最優個體。

最后當然還是需要目標函數才能計算適應度。

算法結果如下圖，小綠點代表種群最優個體適應度：

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=20;   %種群個體的數量

%初始化種群，定義結構體
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=0;               %個體適應度為0初始化
end
par_best=par(1);    %初始化種群中最佳個體

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

select_and_recombine.m

function [next_par par_best]=select_and_recombine(par,par_best,n)
    F=0.5;          %加速因子
    CR=0.8;         %變異率
    next_par=par;   %新種群
    for i=1:n
        
        while 1         %在原種群中任選三個互不相同的個體進行交叉變異
            r1=floor(1+20*rand());
            r2=floor(1+20*rand());
            r3=floor(1+20*rand());
            if i~=r1 && i~=r2 && i~=r3 &&...
               r1~=r2 && r1~=r3 && r2~=r3
                break;
            end
        end
                     
        if rand()<CR    %變異率，可以對每一個特征分別設置，我這里要變一起變了
            next_par(i).x=par(r1).x+F*(par(r2).x-par(r3).x);    %交叉變異准則
            next_par(i).y=par(r1).y+F*(par(r2).y-par(r3).y);              
        end
        
        %計算變異后個體的適應度
        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