基于群智能算法的函数最值优化问题

基于群智能算法的函数最值优化问题

摘要
针对求解函数的极值问题的群智能算法，大部分函数在定义域内都可以通过微分的方式求得极值点和找到最值。但是存在一些测试函数，他们的最值用求解的微分的方式只能使得计算量大大增加。随着计算机技术的发展，可以通过matble编写基于群智能算法的随机搜索算法，不断去逼近极值，不断地优化，改变控制参数，来找到全局最优解。

关键词：遗传算法；粒子群优化算法；全局最优解；控制参数选择；
一、遗传算法(genetic algorithms，GA)
遗传算法一类借鉴生物界自然选择和自然遗传机制的随机搜索算法。鉴于生物在自然界的生存环境中适者生存、不适者被淘汰的过程。具有遗传特征，子代总是与亲代(父代)相似。
所以遗传算法非常适用于处理传统搜索方法难以解决的复杂和非线性优化问题。
在本题目中:使用遗传算法求解给定区间内Sphere函数的最小值问题。
当Sphere的变量为两个时，其在空间中的形状为圆形抛物面，使用matble作图可得，其最小值的位置在原点处。

（一）、遗传算法的基本思想
适者生存，最适合的的个体生存和繁殖后代的可能比较大，在求解的过程中先随机设置一个种群，多个个体，即多个个体，后根据适应程度进行生存和繁殖后代。

1、 编码
在这里采用二进制编码的方式，每一个基因由两个dna13行1列矩阵拼接而成，前十位后后十位分别表示函数中的点的横纵坐标，采用这种编码方式可以在基因的遗传和变异中按位进行操作。更能体现其的随机性。:类似生物染色体的组成，算法易于用生物遗传理论解释，遗传操作如交叉、变异等易实现;算法处理的模式数最多。
其中13基因中每一位作用：
1：符号位，1表示负数，0表示正数
2-8：横纵坐标的整数部分，其中7位二进制数最大可表示到128，对于本题目中要求的定义域【-100，100】以可以满足
9-13：横纵坐标的小数部分，
使用randi函数产生两个，其中5位二进制数最小可表示到0.03125，对于一些函数，极值点少，且非常明显的已经可以满足要求

    dna1 = randi([0, 1], [N, L]);   % 基因100hang 13lie的均匀矩阵
    dna2=  randi([0, 1], [N, L]) ;    % 基因100hang 13lie的均匀矩阵

2、解码decode
解码就是将基因中的二进制数字转化成为实数，以便可以带入具体的函数中
解码方式依据二进制转十进制的规则，一次在每一位上加上对应的权重

%解码
function numbers = decode(t, row)
randshu = [];
for i = 1: row
% 整数部分
judgement = t(i, 1);
bit1 = 7;
bit_integer = t(i,2:bit1+1);
% 小数部分
bit2 = 5;
bit_decimal = t(i,bit1+2:13);
integer = 0;
decimal = 0;
% 计算整数部分
for i = 1 : bit1
    integer = integer + bit_integer(i) * (2^(bit1 -i));
end
% 计算小数部分
for j = 1 : bit2
    decimal = decimal + bit_decimal(j) * (2^(-j));
end
% 整合
num = integer + decimal;

num = sum(num);
if judgement == 1
    num = -num;
end
randshu = [randshu num];

end

numbers = randshu;
end

3、初始种群的产生
给定x，y的范围-100-100，使用linspace函数从-128-128之间均匀的取100个数，然后使用surface函数画出函数的等高线图，以便于可以在标出初始种群和子代种群

    x=linspace(xRange(1), xRange(2), 100); %一列数100个从-128-128
    y=linspace(yRange(1), yRange(2), 100); %一列数100个从-128-128
    [xx, yy]=meshgrid(x, y); %2-D and 3-D grids
    surface( xx , yy ,Fx(xx, yy ) )  %Create surface object

然后随机产生两个矩阵，每个矩阵都是100行13列的，这就表示100个初始个体，作为最基本的种群，使用随机函数可以使最初的种群充满不确定性，可以为以后的寻找全局最优解提供初始条件。

4、适应度函数设计
适应度函数：将目标函数映射成适应度函数(fitness function)的方法:
若目标函数f(x)为最大化问题，则适应度函数取为：Fit( f(x)) = f(x)
若目标函数f(x)为最小化问题，则适应度函数取为：Fit( f()x)= - f(x)
在这里我们求解的是最小值问题，既可以转化成Fit( f()x)= - f(x)
判断个体优良与否的准则是各个个体的适应度值:个体适应度越高，其被选择的机会就越多。将适应度值放到矩阵的最后一列，就可以对种群的基因进行排序，同样依据适应度比例方法

%适应度函数
function result=Fx(x , y)

    result  = x.^2 + y.^2; %最小值
     
end

5、选择个体（轮盘赌选择）

按照个体顺序求出每个个体的累积概率。
产生一个随机数，它落入累积概率的哪个区域就选择相应的个体交叉。
6、交叉产生子代个体

一点交叉:在个体串中随机设定一个交叉点，实行交叉时，该点前或后的两个个体的部分结构进行互换，并生成两个新的个体。
在这里我选用一点交叉，可以减少计算的复杂程度
每一个基因dna （x,y）是一个1行26列的行矩阵，每次都有一个交叉概率，产生一个随机数，若其小于这个概率，就进行交叉。
交叉中先产生0-26的一个随机数，在随机选取一个一个基因与其进行一点交叉，充分的体现了交叉的随机性与灵活性。

%交叉函数，一条dna，产生两条交叉后的dna
function [x1,x2]=genetic(gene, pc)
    global N L num % n = 100, l  = 13, num = 2
    x1 = zeros(N, num*L);                % 初始化子代基因，提速用
    x2 = x1;                         % 同上
    for i = 1: N                 % 交叉操作
        if rand < pc % rand返回一个0-1的随机数
           d = randi(N);    % rand返回一个0-N的随机数
           m = gene(d,:);            % 随机抽选一个交叉的个体，取一行dna
           d = randi(num*L-1);          % 确定交叉断点 1-26随机去一个数字
           x1(i,:) = [gene(i,1:d), m(d+1: num*L)];  % 新个体 1   前半部分是源基因的，后半部分是随机抽取的基因后半部分     
           x2(i,:) = [m(1:d), gene(i, d+1: num*L)];  % 新个体 2  前半部分是随机抽取基因的，后半部分源基因后半部分 
        end
    end
end

7、变异
位点变异:对群体中的个体码串随机挑选一个或多个基因座，并对这些基因座的基因值以变异概率作变动。
若某位原为“0”，通过变异操作变成了“1”也可以变异还是“0”。
因为变异的概率较小，同样也采用遍历的方式，产生一个随机数，如小于变异概率就执行变异过程。

%变异
function result=variation(gene, pm)
    global N L num
    for i = 1: N                           % 变异操作
        if rand < pm
            gene(i,randi(num*L)) = randi([0, 1]);
            %产生一个i行num*L列的变异基因,只替换其中的一个基因
        end
    end
    result=gene;

end

8、子代种群自然选择
对适应度进行降序排名，进行种群的自然选择，根据自然选择与种群生存环境的限制，只能保留50条基因，同样也会淘汰，模仿自然现象，也可能会淘汰优质基因。

%选择
function result=choice(gene, fi)
    % 自然选择--排名法
    global L N num
    dna = flipud(sortrows(gene, L*num + 1));   % 对适应度进行降序排名，
    %sortrows(A,column) 基于向量 column 中指定的列对 A 进行排序。
    % flipud()上下翻转矩阵
    
    while size(dna, 1) > N                 % 自然选择,只保留50条基因
        d = randi(size(dna, 1));           % 排名法 size（dna，1）表示行数，介于0-行数之间的一个随机数
        %当d非常大时，则容易被淘汰
        if rand < (d - 1) / size(dna, 1)  %根据概率大小，淘汰排名靠后的小概率基因 
            dna(d,:) = [];
            fi(d, :) = [];
        end
    end
    result=dna ; 
end

（二）、解决情况与种群分布
其中黑色点表示自带种群，蓝色点表示初始种群
N = 2 n = 3

N = 5 n = 10

（三）、总结

当求解函数比较简单时，可以很快的求解出最优解，且繁殖次数也比较少，可以很容易得出结论。不断去逼近极值，不断地优化，改变控制参数，来找到全局最优解。
在全局最优解和局部最优解的的判断上，基本将种群的规模扩大时，可以得到合适的全局最优解。
1、模型评价优化
总之，遗传算法是一个比较黑的黑箱模型，通过不断地计算，优胜略汰。但是大致可以说，注意，“大致可以说”，只有那些稳定性和单个区间的实力兼备的基因，最后才有可能在残酷的竞争中占据主导地位，稳定性和实力必须兼顾。
关于遗传运算的终止进化代数，设定一个计数器，如果连续N代出现的最优个体的适应度都一样时，（严格的说应该是，连续N代子代种群的最优个体适应度都<=父代最优个性的适应度）可以终止运算。不用重复多次无用的计算。

同时用其他测试函数进行测试时
2、easom函数

查资料可得esaom函数的极值点为（π，π）

很多代后可以基本确定，子代在极值点附近活动，但是每次距离极值点都有一些距离
分析可得
由于采用二进制编码的方式，最小的精度为
0.03125不能很好的逼近π，所以对于这种函数采用二进制编码就感觉有写吃力。
如果小数位的精度选取更多的话就会，大大增加计算量，使得函数的左右有写得不尝试

3、 himmelblaus函数

容易得知，此函数存在四个极值点，故使用遗传算法存在局部最优解问题

经过多次迭代可以看出最小值在（3， 2）处取得，但是 需要次数比较多，总之多极值点问题也基本可以解决，改变同的交叉概率和变异概率可以得到不同效果
4、 其他函数
z = xcos(4x)^2 + y*sin(y) 基本也是可以解决

二、粒子群优化(particle swarm optimization，PSO)算法
1、基本思想
将群体中的每个个体看作n维飞行通过群体有质量的粒子，在搜索空间中以一定的速度飞行，通过群体粒子间的合作与竞争产生的群体智能指导优化搜索。每个粒子有一个由适应度函数决定的适应值。
2、基本原理:
PSO初始化为一群随机粒子，然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个“极值”来更新自己。第一个就是粒子本身所找到的最优解，这个解称为个体极值。另一个是整个种群目前找到的最优解，这个解称为全局极值。

（一）、遗传算法的基本思想

上式右边的第1部分称为“记忆项”，表示粒子在前一时刻的速度对下一时刻速度的影响;
第2部分为个体“认知(cognition)”分量，表示粒子本身的思考，将现在的位置和曾经经历过的最优位置相比。
第3部分是群体“社会(social)”分量,表示粒子间的信息共享与相互合作。
和p2分别控制个体认知分量和群体社会分量相对贡献的学习率(学习因子)。
随机系数增加认知和社会搜索方向的随机性和算法多样性
（二）、解决情况与种群分布
N = 20 n = 400

（三）、总结
PSO算法的搜索性能取决于其全局探索和局部细化的平衡，这在很大程度上依赖于算法的控制参数，包括粒子群初始化、惯性因子w、最大速度和加速常数 等。
1、PSO算法优点：

1)需要调整的参数少，原理简单，容易实现。
2)协同搜索，同时利用个体局部信息和群体全局信息指导搜索。
3)对于目标函数仅能提供极少搜索最优值的信息， PSO算法的粒子具有飞越性的特点使其能够跨过搜索平面上信息严重不足的障碍，飞抵全局最优目标值。
2、PSO算法的缺点：

1)算法局部搜索能力较差，搜索精度不够高。

2)算法不能绝对保证搜索到全局最优解，主要有两方面的原因：

1有时粒子群在俯冲过程中会错失全局最优解。粒子飞翔过程中的俯冲动作使搜索行为不够精细，不容易发现全局最优目标值，所以对粒子的最大飞翔速度进行限制既是为了使粒子不要冲出搜索区域的边界，同时也是为了使搜索行为不至于太粗糙。
2算法搜索性能对参数具有一定的依赖性。对于特定的优化问题，如果用户经验不足，参数调整的确是个棘手的问题。参数值的大小直接影响到算法是否收敛以及求解结果的精度。
三、遗传算法与粒子群算法的比较
用遗传算法求解函数比较简单时，可以很快的求解出最优解，且繁殖次数也比较少，可以很容易得出结论。不断去逼近极值，不断地优化，改变控制参数，来找到全局最优解。在全局最优解和局部最优解的的判断上，基本将种群的规模扩大时，可以得到合适的全局最优解。
但是算法的实现比较复杂，但是解决问题的效率和准确性都比较高。

PSO算法局部搜索能力较差，搜索精度不够高。算法不能绝对保证搜索到全局最优解， 有时粒子群在俯冲过程中会错失全局最优解。
粒子飞翔过程中的俯冲动作使搜索行为不够精细，不容易发现全局最优目标值，所以对粒子的最大飞翔速度进行限制既是为了使粒子不要冲出搜索区域的边界，同时也是为了使搜索行为不至于太粗糙。
同时算法搜索性能对参数具有一定的依赖性。对于特定的优化问题，如果用户经验不足，参数调整的确是个棘手的问题。参数值的大小直接影响到算法是否收敛以及求解结果的精度。
但是粒子群算法原理简单，容易实现。

四、参考文献
[1]李德毅.人工智能导论[M].安徽:中国科学技术大学出版社,2018.83-88
**

最后复完整代码，供大家学习，大家多多支持

**

function [X,Y,Fxy]=GA( )
    clear
    clc
    global N L xRange yRange num
    xRange=[-128, 128]; % x的范围
    yRange=[-128 ,128];% y的范围 

    num=2 ;%变量个数
    N = 100;                          % 种群上限
    ger = 5;                  % 迭代次数
    L = 13;                           % 单个变量的编码基因长度
    pc = 0.8;                        % 交叉概率
    pm = 0.1;                        % 变异概率
    x=linspace(xRange(1), xRange(2), 100); %一列数100个从-128-128
    y=linspace(yRange(1), yRange(2), 100); %一列数100个从-128-128
    [xx, yy]=meshgrid(x, y); %2-D and 3-D grids
    surface( xx , yy ,Fx(xx, yy ) )  %Create surface object

    dna1 = randi([0, 1], [N, L]);   % 基因100hang 13lie的均匀矩阵
    dna2=  randi([0, 1], [N, L]) ;    % 基因100hang 13lie的均匀矩阵
    hold on
    x = decode(dna1,100)';  % 对初始种群解码
    y = decode(dna2,100)';
    %plot3(x, y, Fx(x, y), '*', 'linewidth',3) % 画出初始解的位置
    
    dna=[dna1,dna2]; %合并两个变量的基因

    for epoch = 1: ger               % 进化代数为50
        %交叉
        [x1,x2]=genetic(dna, pc);  % pc = 0.8
        % 变异操作
        x3=variation(dna, pm); 
        
        dna = [dna; x1; x2; x3];          % 合并新旧基因
        dna1=dna(:, 1:L);
        dna2=dna(:, L+1 : 2*L);
        
        fi = zeros(400, 1); % 初始化适应度，提速 
        
        xx = decode(dna1, 400)';  % 对种群解码
        yy = decode(dna2, 400)';
       
        fi=-Fx(xx , yy); %求最小值
       %求最小值，或者最大值，只需将适应度取逆即可！y
      %{
      size_x = size(x)
      size_y = size(y)
      size_dna = size(dna)
      size_fi = size(fi)
      %}
        dna = [dna, fi]; %适应度放在最后一列  
        %自然选择
        dna=choice(dna, fi);
        dna = dna(:, 1:num*L);  %去掉矩阵末尾的适应度，进入下一次迭代
    end
    
%     gene=dna(1,:); %最优基因
    x = decode( dna(:, 1:L) ,100)';         % 对最终种群解码 (dco(1)-1)
    y = decode( dna(:, L+1:num*L ) ,100)';
    plot3(x, y, Fx(x, y) , 'ko' ,'linewidth', 3)           % 画出最终解的位置
%     fprintf('最优解为x=%.10f,y=%.10f\n', x(1),y(1) );
%     fprintf('最优值为z=%.10f\n',fi(1));
    X_m = x(1)
    Y_m = y(1)
    Fxy_m=Fx(X_m,Y_m)
    
end

%解码
function numbers = decode(t, row)
randshu = [];
for i = 1: row
% 整数部分
judgement = t(i, 1);
bit1 = 7;
bit_integer = t(i,2:bit1+1);
% 小数部分
bit2 = 5;
bit_decimal = t(i,bit1+2:13);
integer = 0;
decimal = 0;
% 计算整数部分
for i = 1 : bit1
    integer = integer + bit_integer(i) * (2^(bit1 -i));
end
% 计算小数部分
for j = 1 : bit2
    decimal = decimal + bit_decimal(j) * (2^(-j));
end
% 整合
num = integer + decimal;

num = sum(num);
if judgement == 1
    num = -num;
end
randshu = [randshu num];

end

numbers = randshu;
end

%适应度函数
function result=Fx(x , y)

    result  = x.^2 + y.^2; %最小值
     
end


%选择
function result=choice(gene, fi)
    % 自然选择--排名法
    global L N num
    dna = flipud(sortrows(gene, L*num + 1));   % 对适应度进行降序排名，
    %sortrows(A,column) 基于向量 column 中指定的列对 A 进行排序。
    % flipud()上下翻转矩阵
    
    while size(dna, 1) > N                 % 自然选择,只保留50条基因
        d = randi(size(dna, 1));           % 排名法 size（dna，1）表示行数，介于0-行数之间的一个随机数
        %当d非常大时，则容易被淘汰
        if rand < (d - 1) / size(dna, 1)  %根据概率大小，淘汰排名靠后的小概率基因 
            dna(d,:) = [];
            fi(d, :) = [];
        end
    end
    result=dna ; 
end


%变异
function result=variation(gene, pm)
    global N L num
    for i = 1: N                           % 变异操作
        if rand < pm
            gene(i,randi(num*L)) = randi([0, 1]);
            %产生一个i行num*L列的变异基因,只替换其中的一个基因
        end
    end
    result=gene;

end

%交叉函数，一条dna，产生两条交叉后的dna
function [x1,x2]=genetic(gene, pc)
    global N L num % n = 100, l  = 13, num = 2
    x1 = zeros(N, num*L);                % 初始化子代基因，提速用
    x2 = x1;                         % 同上
    for i = 1: N                 % 交叉操作
        if rand < pc % rand返回一个0-1的随机数
           d = randi(N);    % rand返回一个0-N的随机数
           m = gene(d,:);            % 随机抽选一个交叉的个体，取一行dna
           d = randi(num*L-1);          % 确定交叉断点 1-26随机去一个数字
           x1(i,:) = [gene(i,1:d), m(d+1: num*L)];  % 新个体 1   前半部分是源基因的，后半部分是随机抽取的基因后半部分     
           x2(i,:) = [m(1:d), gene(i, d+1: num*L)];  % 新个体 2  前半部分是随机抽取基因的，后半部分源基因后半部分 
        end
    end
end