利用遗传算法求解函数极值

1.利用遗传算法求解函数极值

例1 利用遗传算法求函数
f(x) = 11sin(6x) + 7cos(5x)，x∈[- π,π]的最大值点。

解：在MATLAB中编制绘制函数曲线的代码，运行得到题中函数的曲线图如下图所示。从下图中可以看出,该函数有多个极值点，如果使用其他的搜寻方法,很容易陷人局部最小点,而不能搜寻到真正的全局最小点,但遗传算法可以较好地弥补这个缺陷。

遗传算法的具体实现如下。

1. 问题分析

在本题中,设定自变量x为个体的基因组,即用二进制编码表示x;设定函数值
f(x)为个体的适应度,函数值越小,适应度越高。关于二进制编码方式，在精度允许的范围内,可以将区间内的无穷多点用间隔足够小的有限点来代替,以降低计算量同时保证精度损失不大。

1. 编写代码方法

用长度为十的二进制编码串来分别表示两个决策变量U_max (2π)、U_min(0)。10位二进制编码串可以表示为0~1023之间的1024个不同的数,故将Umax、Umim的定义域离散为1023个均等的区域，包括两个端点在内共有1024个不同的离散点。

从离散点Umin到Umax ,依次让它们分别对应从0000000000(0)到1111111111(1023)之间的二进制编码。再将分别表示Umax、Umin的二进制编码串联在一起,组成一个20位长的二进制编码串,它构成了这个函数优化问题的染色体编码方法。使用这种编码方法,解空间和遗传算法的搜索空间具有一一对应的关系。

1. 确定个体评价方法

在遗传算法中,以个体适应度的大小来确定该个体被遗传到下一代群体中的概率。

基本遗传算法使用比例选择算子来确定群体中各个个体遗传到下一代群体中的数量。为正确计算不同情况下各个个体的遗传概率,要求所有个体的适应度必须为整数或者是零,不能是负数。

为满足适应度取非负值的要求,将目标函数f(x)变换为个体适应度F(x),一般采用以下两种方法。

• 对于求目标函数最大值的优化问题,变换公式为

F(x) =f(x) +C_min， (f(x)+ C_min) > 0
F(x)=0 其他
其中,Cmin为一个适当地相对比较小的数,它可用预先指定一个较小的数,或者进化到当前代为止的最小目标函数值、当前代或最近几代群体中的最小目标函数。

• 对于求目标函数最小值的优化问题,变换公式为

F(x) = C_max-f(x)， f(x) < C_max
F(x) = 0 其他
其中,Cmax为一个适当地相对比较大的数，它可用预先指定一个较大的数，或者进化到当前代为止的最大目标函数值、当前代或最近几代群体中的最大目标函数值。

根据适者生存原则选择下一代的个体。在选择时,以适应度为选择原则。适应度准则体现了适者生存,不适应者淘汰的自然法则。

1. 设计遗传算子

(1)使用赌轮选择算法,决定哪些个体可以进入下一代,其过程如下。
①在第t代,由p_i= f_i/Σf_i= f_i/f_sum ， 计算f_sum和p_i。
②产生{0,1}的随机数rand(),求s=rand()* f_sum.
③求Σf_i≥s中最小的k ,则第k个个体被选中。
④进行N次②,③操作,得到N个个体,成为第t=t+1代种群。
(2)分别求出20个初始种群中每个种群个体的适应度函数,并计算所有种群的和S.
(3)在区间(0,S)上随机的产生一个数r从某个基因开始,逐一取出基因来,把它的适应度加到s上去(s开始为0),如果s大于r,则停止循环并返回当前基因.

1. 交叉运算

使用单点交叉算子,对于选中用于繁殖下一代的个体,随机地选择两个个体的相同位置,按交叉概率P在选中的位置实行交换。这个过程反映了随机信息交换,目的在于产生新的基因组合也即产生新的个体。
假设2个父代个体x1、x2为.
x1 = 0100110
x2 = 1010001
从每个个体的第3位开始交叉,交叉后得到两个新的子代个体y1,y2分别为
y1 =0100001
y2 = 1010110
这样两个子代个体就分别具有了两个父代个体的某些特征。

1. 变异运算

在所有的个体一样时,交叉是无法产生新个体的，这时只能靠变异产生新的个体。

在算法中，可以变异概率分别对每个或者几个基因位做变异操作就可以生成新的种
群个体。

父代中的每个个体的每一位都可以以设定的概率进行翻转,即由1变为0,或由0变
为1。变异增加了全局优化的特质。

根据以上分析,编写MATLAB代码如下:

主程序

clc;
clear all;
popsize = 20;       %设置初始参数，群体大小
chromlength = 8;    %字符出长度（个体长度），染色体长度
pc = 0.7;           %设置交叉概率
pm = 0.02;          %设置变异概率
pop = initpop(popsize,chromlength);     %运行初始化函数，随机产生初始群体
for i = 1:20
    [objvalue] = calobjvalue(pop);      %计算目标函数
    fitvalue = calfitvalue(objvalue);   %计算群体中每个个体中的适应度
    [newpop] = selection(pop,fitvalue); %复制
    [newpop] = crossover(pop,pc);       %交叉
    [newpop] = mutation(pop,pc);        %变异
    [bestindividual,bestfit] = best(pop,fitvalue);     %求出群体中适应值最大的个体及其适应值
    y(i) = max(bestfit)
    n(i) = i;
    pop5 = bestindividual;
    x(i) = decodechrom(pop5,1,chromlength)*10/1023
    pop = newpop;
end

fplot(@(x)11.*sin(6.*x)+7.*cos(5.*x));
grid on;
hold on;
plot(x,y,'r* ');
xlabel('自变量');
ylabel('目标函数值');
title('种群中的个体数目为20，二进制编码长度为8,交叉概率为0.7,变异概率为0.02');
fmax = max(y)
hold off

子函数

%初始化
function pop = initpop( popsize, chromlength)
pop = round(rand(popsize, chromlength));
%rand随机产生每个单元为{0,1}行数为popsize,列数为chromlength的随机数矩阵,
%roud对矩阵的每个单元进行圆整，这样产生的初始种群，
end
%round函数用于舍入到最接近的整数。语法形式只有1种：Y = round(X)，这里的X可以是数，向量，矩阵，输出对应。

%将二进制编码转换成十进制
function pop2 = decodebinary(pop)
[px,py]= size(pop);
%求pop行和列数
for i= 1:py
    pop1(:,i)=2.^(py - i).* pop(:,i);
end
pop2 = sum(pop1,2);
%求popl的每行之和,1为列的和
end

%将二进制编码转换成十进制
function pop2 = decodechrom( pop, spoint, length)
    pop1 = pop(:, spoint:spoint + length- 1);
    pop2 = decodebinary(pop1);
end

%实现目标函数的计算
function  [abjvalue] = calobjvalue(pop)
    templ = decodechrom(pop,1,8); %将pop每行转化成十进制数
    x= templ * 10/1023;           %将二值域中的数转化为变量域的数
    abjvalue = 11* sin(6*x)+7* cos(5*x);%计算目标函数值
end

%计算个体的适应值
function fitvalue = calfitvalue(objvalue)
global Cmin;
Cmin= 0;
[px, py]= size(objvalue);
for i = 1:px
    if objvalue(i) + Cmin >0
        temp = Cmin + objvalue(i);
    else
        temp = 0.0;
    end
        fitvalue(i) = temp;
end
fitvalue = fitvalue';
%%使用golbal的优点：

% 1 传递大数据的参数
%如果通过函数传参数的方式的话，系统会浪费过多的时间在复制数据的时间上，如果采用global的方式共享数据的话代码的效率会大大提高
% 
% 2 过多的常量需要传递
%如果每个量都作为函数函数的参数传递的话，代码参数列表就很长，如果采用global的话代码的可读性提高，函数调用也方便

function [newpop] = selection(pop, fitvalue)
totalfit = sum(fitvalue);       %求适应值之和
fitvalue = fitvalue/ totalfit;  %单个个体被选择的概率
fitvalue = cumsum(fitvalue);    %如fitvalue=[1 2 3 4],则cumsum(fitvalue)=[1 3 6 10]
[px,py]= size(pop);
ms = sort(rand(px,1));         %从小到大排列
fitin = 1;
newin = 1;
while newin <= px
    if (ms(newin)) < fitvalue(fitin)
        
        newpop(newin) = pop(fitin);
        newin= newin + 1;
    else
        fitin = fitin+1;
    end
end
end

%交叉
function [newpop] = crossover(pop, pc)
[px, py] = size(pop);
newpop = ones(size(pop));
for i= 1:2:px- 1  %2为间隔
    if(rand <pc)
        cpoint = round(rand* py);
        newpop(i,:) = [pop(i,1:cpoint),pop(i+ 1,cpoint+1:py)];
        nevpop(i+1,:)= [pop(i + 1,1:cpoint),pop(i,cpoint+ 1:py)];
    else
        newpop(i,:) = pop(i);
        newpop(i+ 1,:)=pop(i+1);
    end
end

%变异
function [newpop] = mutation(pop,pm)
[px,py]= size(pop); 
newpop = ones(size(pop));
for i= 1:px
    if(rand< pm)
        mpoint = round(rand * py);
    if mpoint<= 0
        mpoint= 1;
    end
        newpop(i) = pop(1);
            if any( newpop( i, mpoint))==0
                newpop( i, mpoint) = 1;
            else
                newpop(i, mpoint) = 0;
            end
    else
        newpop(i) = pop(i);
    end
end

%计算群体中最大的适应值及其个体
function [ bestindividual, bestfit] = best(pop, fitvalue)
[px,py]= size(pop);
bestindividual = pop(1,:);
bestfit= fitvalue(1);
for i= 2:px
    if fitvalue(i)> bestfit
        bestindividual= pop(i,:);
        bestfit= fitvalue(i);
    end
end

• 本题中选择二进制编码,种群中的个体数目为20,二进制编码长度为8,交叉概率为0.7,变异概率为0.02。 运行以上程序得到结果如下：

结果如下



这是用常规求解函数最大值，可以看出最优解相差不大，

• 选择二进制编码,降低交叉概率，即种群中的个体数目为20,二进制编码长度为8,交叉概率为0.2,变异概率为0.02。运行主程序,得到结果：

由图可见,交叉概率较小时,全局最优解也减少。

• 选择二进制编码,增加变异概率，即种群中的个体数目为20,二进制编码长度为8, 交叉概率为0.7,变异概率为0.01。运行主程序,得到结果
  
  
  由图可见,增加变量概率并不一定使全局最优解增加。图中函数最大输出值为17.7533,这与变异概率为0.02时的值相同。

• 选择二进制编码,增加种群个体数目,即种群中的个体数目为30,二进制编码长度为 8,交叉概率为0.7,变异概率为0.02。运行主程序,得到结果如下:
  
  
  由图可见,增加种群个体数目后,函数最大值变小。