matlab实现遗传算法——以Ras函数为例( 初学子 友好子 )

  • Post author:
  • Post category:其他




遗传算法



演化思想



遗传算法的本质

遗传算法的本质=

模拟

生物演化。具体来说,模拟对象是生物演化中的种种自然现象(如变异,交叉互换,交配,淘汰)。因此,一个优秀的遗传算法首先应该做到生物演化的模拟。但是并非仅仅对生物演化进行简单复现,生物演化中有种种细节可以忽略不计,如DNA非编码段变异,等等

那演化中最重要的部分是什么?


对下一代种群的定向选择

。这句话有三个要点,一个是下一代,另一个是定向,最后一个是选择。先说后两者,所谓定向就是人为设定衡量不同个体的可计算指标;而选择则是保留满足指标的个体;最后一点,下一代怎么产生呢?就是遗传,变异和交叉互换。

生物学有关名词这里就不解释了,假设高中生物老师都讲过。



Ras函数

Ras函数(Rastrigin’s Function)是一个典型的非线性

多峰

函数,具有多个局部的极大值点和极小值点,常规算法很容易陷入局部最优解中。

在n维空间中,Ras函数表达式为:





f

(

x

)

=

A

n

+

i

=

1

n

 

[

x

i

2

A

c

o

s

(

2

π

x

i

)

]

f(x)=An+\sum_{i=1}^{n}\ [x^{2}_{i}-Acos(2\pi x_{i})]






f


(


x


)




=








A


n




+

















i


=


1



















n




















[



x











i










2






























A


c


o


s


(


2


π



x











i



















)


]







其中,



A

A






A





是可以人为定义的常量,



n

n






n





是Ras函数的维数





n

=

2

n=2






n




=








2





时,Ras函数为二维形式,这个函数在



(

0

,

0

)

(0,0)






(


0


,




0


)





有全局最小值。

下图是



A

=

10

,

n

=

2

A=10,n=2






A




=








1


0


,




n




=








2





时,Ras函数在



x

,

y

[

5

,

5

]

x,y\in[-5,5]






x


,




y













[





5


,




5


]





上的三维图象

hello



遗传算法




目标函数

objective function

%%
% objective function
function ans = ras(A,x) 
    ans = 0;
    for i=1:2
        ans = ans+A+x(i)^2-A*cos(2*pi*x(i));
    end
end
%%

目标函数是Ras函数。计算这个函数需要三个参数,分别是



A

A






A





,



x

x






x





,



y

y






y





。其中



A

A






A





是我们自行设定的常数,而



x

y

x,y






x





y





则对应二维平面上的点

补充:


function

——->

matlab官方文档-function



fsurf

——->

matlab官方文档-fsurf




初始化种群

initialization population

% initialization population - binary form
function x = init (pop_size,chromolength)
    x = round(rand(pop_size*2,chromolength));
end
  • 初始种群要足够随机且多样:目的是尽可能扩大搜索范围,找到最优解
  • 染色体长度等于初始种群中每个个体的二进制编码长度




解码

decoding:from binary to decimal

% coding - binary to decimal floating point
function x = decoding(pop)
    [r,c] = size(pop);
    vector = ones(16,1); % 0~65535
    for i=1:c/2
        for j=1:(c/2-i)
            vector(i) = vector(i)*2;
        end
    end
    pop_xy = mat2cell(pop',[c/2 c/2]);
    pop_x = pop_xy{1,1}'*vector;  %1*16
    pop_y = pop_xy{2,1}'*vector; %1*16
    x=[pop_x pop_y]/65535;
end

每个个体由32位二进制数字组成。32位数字从中间划分成两组,前16位表示x,后16位表示y,组合起来代表二维平面上的点。

  • mat2cell的时候注意元胞数组的索引方式

    每个个体与向量



    v

    e

    c

    t

    o

    r

    vector






    v


    e


    c


    t


    o


    r





    相乘,可以将16位二进制数转化成十进制数




    v

    e

    c

    t

    o

    r

    =

    [

    2

    15

    2

    14

    ,

    ,

    2

    ,

    1

    ]

    vector=[2^{15} ,2^{14},…,2,1]






    v


    e


    c


    t


    o


    r




    =








    [



    2











    1


    5














    2











    1


    4










    ,









    ,




    2


    ,




    1


    ]







    得到的乘积是一个范围落在在



    (

    0

    ,

    65535

    )

    (0,65535)






    (


    0


    ,




    6


    5


    5


    3


    5


    )





    内的整数。为了将这个范围缩小到10附近,再将结果除以6553。

    函数最后输出一个50*2的矩阵,储存50个个体的二维坐标值

    补充:


    mat2cell

    函数——->

    matlab官方文档-mat2cell



    ones

    函数——->

    matlab官方文档-ones




计算目标函数

calculate the real value

% calculate the real value
function obj = calobjvalue(pop2)
    [r,c] = size(pop2);
    obj = [ ];
    for i=1:r
        for j=1:c
            obj(i) = ras(pop2(i,:));
        end
    end
    obj =obj';  % 1*50 objective value in floating point
end




计算适应度

calculate the fittness value

% calculate the fittness value; range in 0 and 1
function fit = calfitvalue(pop_obj)
    [r,c] = size(pop_obj);
    % calculate the max & min value
    for i=1:r/2 % divide the pop_obj into two sets % 注意对矩阵行列的索引是不是弄反了
        if pop_obj(i) > pop_obj(r+1-i)
            maxSet(i) = pop_obj(i);
            minSet(i) = pop_obj(r+1-i);
        else
            maxSet(i) = pop_obj(r+1-i);
            minSet(i) = pop_obj(i);
        end
    end
    pop_max = maxSet(1);
    pop_min = minSet(1);
    for i=2:r/2 %compare two sets separately
        if maxSet(i) > pop_max
            pop_max = maxSet(i);
        elseif minSet(i) < pop_min
            pop_min = minSet(i);
        end
    end       
    for i=1:r
        fit(i) = 1-(pop_obj(i)-pop_min)/(pop_max-pop_min);
    end
    fit = fit';
end

衡量适应度采用常见归一化方法





X

=

x

x

m

i

n

x

m

a

x

x

m

i

n

X = \frac{x-x_{min}}{x_{max}-x_{min}}






X




=




















x











m


a


x



























x











m


i


n































x










x











m


i


n










































如果求解最大值,则



X

X






X





越大,个体适应度越高;

如果求解最小值时,则



(

1

X

)

(1-X)






(


1













X


)





的结果越大,个体适应度越高。




选择最优个体

choose the best individual

% best
function best = bestIndividual(pop)
    popDe = decoding(pop);
    obj = calobjvalue(popDe);
    fit = calfitvalue(obj);
    [r,c] = max(fit);
    best = popDe(r,:);
end




选择复制

selection copy

% evolution
function pop = evolutionpop(pop,fitvalue)
    % roulette wheel selection
    fitValue = fitvalue/sum(fitvalue);
    cumFit = cumsum(fitValue); % 工作是不是遗漏了一步
    [r,c] = size(pop);
    for i=1:r
        dartFather = rand; % select father 缺少交配这一步,仅仅是完全保留优秀个体
       % dartMother = rand; % select mother
        j = 1;
        while dartFather > cumFit(j)
            j = j+1;
        end
%         while dartMother >cumFit(k)
%             k = k+1;
%         end    
        popNew(i,:) = pop(j,:);
    end
    pop = popNew;
end

选择复制的思想是每次都保留种群中的优秀个体,逐代提高下一代种群适应度。特别要注意的是,种群适应度的提高的过程一定要 “

逐步

” ,提高速度太快,种群容易过早陷入局部最优解;提高速度过慢,则迭代次数过多,计算速度慢。

那么,怎样才能体现

逐步

这二字?

答案是,在每一代种群中选择优秀个体的同时,保留一些不那么优秀的个体,为种群提供多样性和可能性。

或者降低交叉互换,基因变异的概率

补充:

轮盘赌——–>

轮盘赌算法



cumsum

——–>

matlab帮助文档-cumsum


进行到这一步,结合计算适应度、计算目标函数、选择最优个体和选择复制,已经可以得到多次演化大致的图象

在这里插入图片描述

图象在前10次迭代里快速向高适应度方向演化,但是在演化10次左右后最小值不再改变,呈现为一条平滑的直线,陷入局部最优解。

用这种不进行交叉互换,不产生基因突变,仅仅依靠选择复制的算法,计算得最优个体为



(

3.3357

,

3.5710

)

(3.3357 ,3.5710)






(


3


.


3


3


5


7


,




3


.


5


7


1


0


)








交叉互换

crossover

function popNew = crossover(pop,pc,fitvalue)
    [r,c] = size(pop);
    fitValue = fitvalue/sum(fitvalue);
    cumFit = cumsum(fitValue);
    for i=1:r
        dart = rand;
        popNew(i,:) = pop(i,:);
        if dart < pc % 选择是否进行交叉互换
            % 选择pop(i,:)的交叉互换对象 pop(j,:)
            dart2 = rand;
            j = 1;
            while dart2 > cumFit(j)
               j = j+1;
            end
             % 确定彼此从哪一位开始交叉互换
            dart3 = (c/2)*rand;
            k = 0;
            while k < dart3
                k = k+1;
            end
            % 进行交叉互换
            for n=1:(k-1)
                popNew(i,n) = pop(i,n);
                popNew(i,n+c/2) = pop(i,n+c/2);
            end
            for m=k:c/2
                popNew(i,m) = pop(j,m);
                popNew(i,m+c/2) = pop(j,m+c/2);
            end
        end
    end
end

从任意位置开始随机选择两个个体交换染色体,进行交叉互换

在这里插入图片描述

交叉互换后得到的解更加逼近于全局最优解,得到的最优个体为



(

2.5010

,

1.5626

)

( 2.5010 ,1.5626 )






(


2


.


5


0


1


0


,




1


.


5


6


2


6


)







随机变异

mutation

% mutation
function popNew = mutation(pop,pm)
    % 同比例放大pm与dart
    [r,c] = size(pop);
    dart = rand;
    for i=1:r
        popNew(i,:) = pop(i,:);
        if dart*10 < pm*10 %判断基因是否变异
            % 选择某位进行变异
            dart2 = rand*c;
            k = 0;
            while k < dart2
                k = k+1;
            end
            if popNew(i,k) == 0
                popNew(i,k) = 1;
            else
                popNew(i,k) = 0;
            end
        end
    end
end

在这里插入图片描述

随机选择染色体中某位1变成0,0变成1,得到的最优个体为



0.0012

0.9962

)

( 0.0012 ,0.9962)









0


.


0


0


1


2





0


.


9


9


6


2


)





,更加逼近全局最优解。




总代码

clear all;
clc;
pc = 0.8; 				% 交叉互换概率
pm = 0.05;				% 基因变异概率
pop_size = 50; 			% 种群大小
chromolength = 32; 		% 染色体长度
times = 100;			% 迭代次数

pop{1,1} = init(pop_size,chromolength); %产生初始种群
for i=1:times %进行一百次种群演化
    pop2 = decoding(pop{i,1});
    pop_obj = calobjvalue(pop2);
    total(i) = sum(pop_obj);
    pop_fit = calfitvalue(pop_obj);
    pop{i+1,1} = evolutionpop(pop{i,1},pop_fit);
    pop{i+1,1} = crossover(pop{i+1,1},pc,pop_fit);
    pop{i+1,1} = mutation(pop{i+1,1},pm);
end
ret = choosebest(pop{times+1,1})
plot(total/pop_size);
%%
% objective function
function ans = ras(x)
    ans = 0;
    A = 10;
    for i=1:2
        ans = ans+A+x(i)^2-A*cos(2*pi*x(i));
    end
end
%%

%%
% initialization population - binary
function x = init (pop_size,chromolength)
    x = round(rand(pop_size,chromolength)); %产生随机01矩阵
end


% coding - binary to decimal floating point
function x = decoding(pop)
    [r,c] = size(pop);
    vector = ones(c/2,1); % 0~65535
    for i=1:c/2
        for j=1:(c/2-i)
            vector(i) = vector(i)*2;
        end
    end
    pop_xy = mat2cell(pop',[c/2 c/2]);
    pop_x = pop_xy{1,1}'*vector;  %1*16
    pop_y = pop_xy{2,1}'*vector; %1*16
    x=[pop_x pop_y]/6553;;
end
%%

% calculate the real value
function obj = calobjvalue(pop2)
    [r,c] = size(pop2);
    obj = [ ];
    for i=1:r
        for j=1:c
            obj(i) = ras(pop2(i,:));
        end
    end
    obj =obj';  % 1*50 objective value in floating point
end

% calculate the fittness value; range in 0 and 1
function fit = calfitvalue(pop_obj)
    [r,c] = size(pop_obj);
    % calculate the max & min value
    for i=1:r/2 % divide the pop_obj into two sets % 注意对矩阵行列的索引是不是弄反了
        if pop_obj(i) > pop_obj(r+1-i)
            maxSet(i) = pop_obj(i);
            minSet(i) = pop_obj(r+1-i);
        else
            maxSet(i) = pop_obj(r+1-i);
            minSet(i) = pop_obj(i);
        end
    end
    pop_max = maxSet(1);
    pop_min = minSet(1);
    for i=2:r/2 %compare two sets separately
        if maxSet(i) > pop_max
            pop_max = maxSet(i);
        elseif minSet(i) < pop_min
            pop_min = minSet(i);
        end
    end       
    for i=1:r
        fit(i) = 1-(pop_obj(i)-pop_min)/(pop_max-pop_min);
    end
    fit = fit';
end

% evolution
function popNew = evolutionpop(pop,fitvalue)
    % roulette wheel selection
    fitValue = fitvalue/sum(fitvalue);
    cumFit = cumsum(fitValue);
    [r,c] = size(pop);
    for i=1:r
        dartFather = rand; % select father 缺少交配这一步,仅仅是完全保留优秀个体
       % dartMother = rand; % select mother
        j = 1;
        while dartFather > cumFit(j)
            j = j+1;
        end
%         while dartMother >cumFit(k) %其实这里是想模拟母本的……
%             k = k+1;
%         end    
        popNew(i,:) = pop(j,:);
    end
    
end

function popNew = crossover(pop,pc,fitvalue)
    [r,c] = size(pop);
    fitValue = fitvalue/sum(fitvalue);
    cumFit = cumsum(fitValue);
    for i=1:r
        dart = rand;
        popNew(i,:) = pop(i,:);
        if dart < pc % 选择是否进行交叉互换
            % 选择pop(i,:)的交叉互换对象 pop(j,:)
            dart2 = rand;
            j = 1;
            while dart2 > cumFit(j)
               j = j+1;
            end
             % 确定彼此从哪一位开始交叉互换
            dart3 = (c/2)*rand;
            k = 0;
            while k < dart3
                k = k+1;
            end
            % 进行交叉互换
            for n=1:(k-1)
                popNew(i,n) = pop(i,n);
                popNew(i,n+c/2) = pop(i,n+c/2);
            end
            for m=k:c/2
                popNew(i,m) = pop(j,m);
                popNew(i,m+c/2) = pop(j,m+c/2);
            end
        end
    end
end

% mutation
function popNew = mutation(pop,pm)
    % 同比例放大pm与dart
    [r,c] = size(pop);
    dart = rand;
    for i=1:r
        popNew(i,:) = pop(i,:);
        if dart*10 < pm*10 %判断基因是否变异
            % 选择某位进行变异
            dart2 = rand*c;
            k = 0;
            while k < dart2
                k = k+1;
            end
            if popNew(i,k) == 0
                popNew(i,k) = 1;
            else
                popNew(i,k) = 0;
            end
        end
    end
end

% best
function best = bestIndividual(pop)
    obj = calobjvalue(pop);
    best = decoding(pop(1,:));
end




未回答的问题

算法虽好,仍有遗憾

  • 陷入局部最优无法跳出窘境

在这里插入图片描述

上图是迭代一千次的结果,每一个细小尖微的波动都是突变在演化中产生的影响,但是仍没有跳出局部最优的窘境。



沿着哲理的山路



局部最优与全局最优



随机



\ne









































=






均匀



退化反而是进化的开始



版权声明:本文为weixin_44634406原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。