基于遗传算法和非线性规划的函数寻优问题

发表于
| 文章字数:2.7k | 阅读时长:11min
matlab
这是一篇更新于 568 天前的文章，其中的信息可能已经有所发展或是发生改变。

本篇文章介绍基于遗传算法和非线性规划的函数寻优算法，目前对于此算法还没有完全的理解掌握，或许有的地方讲的并不是十分正确十分清楚，欢迎指正(●’◡’●)

非线性规划

起源非线性规划是20世纪50年代形成的一门新兴学科。1951年库恩和塔克发表的关于最优线性条件（库恩.塔克条件）的论文是非线性规划诞生的标志。

非线性规划函数

函数fmincon是MATLAB最优化工具箱中求解非线性规划问题的函数。从一个预估值出发，搜素约束条件下非线性多元函数的最小值。

函数fmincon的约束条件为

其中，x、b、beq、lb、和ub是矢量；A和Aeq为矩阵；c(x)和ceq(x)返回矢量的函数；f(x)、e(x)和ceq(x)是非线性函数

函数fmincon的基本用法为

x=fmincon(fun,x0,A,b)

x=fmincon(fun,x0,A,b,Aeq,beq)

x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub)

x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)

x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)

其中，nonlcon为非线性约束条件；lb和ub分别为x的下界和上界。当函数输入参数不包括A,b,Aeq,beq时，默认A=0，b=0,Aeq=[],beq=[].x0为x的初值。

案例背景

采用遗传算法和非线性规划的方法求解如下函数的最小值：
$$
f(x)=-5sinx1.sinx2.sinx3.sinx4.sinx5-sin5x1.5x2.5x3.5x4.5x5+8
$$
其中，x1、x2、x3、x4、x5是0~9π之间的实数。其函数最小值为-2，最小值的位置为（π/2，π/2，π/2，π/2）

非线性规划遗传算法的算法流程图如下所示

遗传算法实现

种群初始化

由于遗传算法不能直接除了问题空间的参数，因此必须通过编码把可行解表示成遗传空间或者染色体个体。常用的编码方法有位串编码、Grey编码、实数编码(浮点法编码)、多级参数编码、有序串编码、结构式编码等。

适应度函数

适应度函数是用来区分群体中个体好坏的标准，是进行自然选择的唯一依据，一般是由目标函数加以变换得到的。本案例是求函数的最小值，把函数值的倒数作为个体的适应度值。函数值越小的个体，适应度值越大，个体越优。适应度计算函数为
$$
F=f(x)
$$

选择操作

选择操作从就旧的群体中以一定的概率选择优良的个体组成新的种群，以繁殖下一代个体。个体被选中的概率跟适应度值有关，个体适应度值越高，被选中的概率越大。遗传算法选择操作有轮盘赌法、锦标赛法等多种方法，此处选择轮盘赌法，即基于适应度比例的选择策略，个体i被选中的概率为
$$
p_i= \frac{F_i}{\sum_{j=1}^N F_j }
$$

$$
其中F_i为个体i的适应度值；N为种群个体数目。
$$

交叉操作

交叉操作是指从种群中随机选择两个个体，通过染色体的交换组合，把父串的优秀特征遗传给子串，从而产生新的优秀个体。由于个体采用实数编码，所以交叉操作采用实数交叉操作法，第k个染色体ak和第l个染色体al在j位的交叉操作方法为：

akj=aij(1-b)+aljb

alj=alj(1-b)+akjb

其中b是[0,1]之间的随机数。

变异操作

变异操作的主要目的是维持种群的多样性。变异操作从种群中随机选取一个个体，选择个体中的一点进行变异以产生更优秀的个体。第i个个体的第j个基因aij进行的变异操作方法为：
$$
\begin{aligned}
f_Y(y) & = f_X[h(y)]|h’(y)| \[2ex]
& = f_X[h(y)]h’(y) \[2ex]
& = \frac{1}{\theta}e^{-\frac{x}{\theta}}[\frac{dx}{dy}(-\frac{\theta}{ln(1-y)})] \[2ex]
& = \frac{1}{\theta}e^{-\frac{-\frac{\theta}{ln(1-y)}}{\theta}}\frac{\theta}{1-y} \[2ex]
& = \frac{1}{\theta}e^{ln(1-y)}\frac{\theta}{1-y} \[2ex]
& = \frac{1-y}{\theta}\frac{\theta}{1-y} \[2ex]
& = 1
\end{aligned}
$$

程序实现

适应度函数：

1
2
3
function y = fun(x)

y=-5*sin(x(1))*sin(x(2))*sin(x(3))*sin(x(4))*sin(x(5))-sin(5*x(1))*sin(5*x(2))*sin(5*x(3))*sin(5*x(4))*sin(5*x(5))+8;

选择操作

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
function ret=Select(individuals,sizepop)
% 本函数对每一代种群中的染色体进行选择，以进行后面的交叉和变异
% individuals input : 种群结构体信息
% sizepop input : 种群规模100个种群
% opts input : 选择方法的选择
% ret output : 经过选择后的种群
individuals.fitness= 1./(individuals.fitness);%%根据适应度值归一
sumfitness=sum(individuals.fitness);
sumf=individuals.fitness./sumfitness;
index=[];
for i=1:sizepop %转sizepop次轮盘
    pick=rand;%%防止初始化为0    
    while pick==0
        pick=rand;
    end
    for j=1:sizepop
        pick=pick-sumf(j);
        if pick<0
            index=[index j];
            break; %寻找落入的区间，此次转轮盘选中了染色体i，注意：在转sizepop次轮盘的过程中，有可能会重复选择某些染色体
        end
    end
end
individuals.chrom=individuals.chrom(index,:);
individuals.fitness=individuals.fitness(index);
ret=individuals;

交叉操作

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
function ret=Cross(pcross,lenchrom,chrom,sizepop,bound)
%本函数完成交叉操作
% pcorss                input  : 交叉概率
% lenchrom              input  : 染色体的长度
% chrom                 input  : 染色体群
% sizepop               input  : 种群规模
% ret                   output : 交叉后的染色体

for i=1:sizepop 
    
    % 随机选择两个染色体进行交叉
    pick=rand(1,2);
    while prod(pick)==0
        pick=rand(1,2);
    end
    index=ceil(pick.*sizepop);
    % 交叉概率决定是否进行交叉
    pick=rand;
    while pick==0
        pick=rand;
    end
    if pick>pcross
        continue;
    end
    flag=0;
    while flag==0
        % 随机选择交叉位置
        pick=rand;
        while pick==0
            pick=rand;
        end
        pos=ceil(pick.*sum(lenchrom)); %随机选择进行交叉的位置，即选择第几个变量进行交叉，注意：两个染色体交叉的位置相同
        pick=rand; %交叉开始
        v1=chrom(index(1),pos);
        v2=chrom(index(2),pos);
        chrom(index(1),pos)=pick*v2+(1-pick)*v1;
        chrom(index(2),pos)=pick*v1+(1-pick)*v2; %交叉结束
        flag1=test(lenchrom,bound,chrom(index(1),:));  %检验染色体1的可行性
        flag2=test(lenchrom,bound,chrom(index(2),:));  %检验染色体2的可行性
        if   flag1*flag2==0
            flag=0;
        else flag=1;
        end    %如果两个染色体不是都可行，则重新交叉
    end
end
ret=chrom;

变异操作

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
function ret=Mutation(pmutation,lenchrom,chrom,sizepop,pop,bound)
% 本函数完成变异操作
% pcorss                input  : 变异概率
% lenchrom              input  : 染色体长度
% chrom                 input  : 染色体群
% sizepop               input  : 种群规模
% pop                   input  : 当前种群的进化代数和最大的进化代数信息
% ret                   output : 变异后的染色体

for i=1:sizepop  
    % 随机选择一个染色体进行变异
    pick=rand;
    while pick==0
        pick=rand;
    end
    index=ceil(pick*sizepop);
    % 变异概率决定该轮循环是否进行变异
    pick=rand;
    if pick>pmutation
        continue;
    end
    flag=0;
    while flag==0
        % 变异位置
        pick=rand;
        while pick==0
            pick=rand;
        end
        pos=ceil(pick*sum(lenchrom));  %随机选择了染色体变异的位置，即选择了第pos个变量进行变异
        v=chrom(i,pos);
        v1=v-bound(pos,1);
        v2=bound(pos,2)-v;
        pick=rand; %变异开始
        if pick>0.5
            delta=v2*(1-pick^((1-pop(1)/pop(2))^2));
            chrom(i,pos)=v+delta;
        else
            delta=v1*(1-pick^((1-pop(1)/pop(2))^2));
            chrom(i,pos)=v-delta;
        end   %变异结束
        flag=test(lenchrom,bound,chrom(i,:));     %检验染色体的可行性
    end
end
ret=chrom;

Code函数

1
2
3
4
5
6
7
8
9
10
11
12
function ret=Code(lenchrom,bound)
%本函数将变量编码成染色体，用于随机初始化一个种群
% lenchrom   input : 染色体长度
% bound      input : 变量的取值范围
% ret        output: 染色体的编码值

flag=0;
while flag==0
    pick=rand(1,length(lenchrom));
    ret=bound(:,1)'+(bound(:,2)-bound(:,1))'.*pick; %线性插值
    flag=test(lenchrom,bound,ret);             %检验染色体的可行性
end

test函数

1
2
3
4
5
6
7
8
9
10
11
function flag=test(lenchrom,bound,code)
% lenchrom input : 染色体长度
% bound input : 变量的取值范围
% code output: 染色体的编码值
flag=1;
[n,m]=size(code);%%n=1 m=5
for i=1:n
    if code(i)<bound(i,1) || code(i)>bound(i,2)%%如果编码后的值在取值范围之内，则flag标记为0
    flag=0;
    end
end

非线性寻优

1
2
3
4
5
6
function ret = nonlinear(chrom,sizepop)

for i=1:sizepop
    x=fmincon(inline('-5*sin(x(1))*sin(x(2))*sin(x(3))*sin(x(4))*sin(x(5))-sin(5*x(1))*sin(5*x(2))*sin(5*x(3))*sin(5*x(4))*sin(5*x(5))'),chrom(i,:)',[],[],[],[],[0 0 0 0 0],[2.8274 2.8274 2.8274 2.8274 2.8274]);
    ret(i,:)=x';
end

算法主函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
%% 清空环境
clc
clear
warning off

%% 遗传算法参数
maxgen=30;                         %进化代数
sizepop=100;                       %种群规模
pcross=[0.6];                      %交叉概率
pmutation=[0.01];                  %变异概率
lenchrom=[1 1 1 1 1];              %变量字串长度
bound=[0 0.9*pi;0 0.9*pi;0 0.9*pi;0 0.9*pi;0 0.9*pi];  %变量范围

%% 个体初始化
individuals=struct('fitness',zeros(1,sizepop), 'chrom',[]);  %种群结构体
avgfitness=[];                                               %种群平均适应度
bestfitness=[];                                              %种群最佳适应度
bestchrom=[];                                                %适应度最好染色体
% 初始化种群
for i=1:sizepop
    individuals.chrom(i,:)=Code(lenchrom,bound);       %随机产生个体
    x=individuals.chrom(i,:);
    individuals.fitness(i)=fun(x);                     %个体适应度
end

%找最好的染色体
[bestfitness bestindex]=min(individuals.fitness);
bestchrom=individuals.chrom(bestindex,:);  %最好的染色体
avgfitness=sum(individuals.fitness)/sizepop; %染色体的平均适应度
% 记录每一代进化中最好的适应度和平均适应度
trace=[];

%% 进化开始
for i=1:maxgen
    
    % 选择操作
    individuals=Select(individuals,sizepop);
    avgfitness=sum(individuals.fitness)/sizepop;
    % 交叉操作
    individuals.chrom=Cross(pcross,lenchrom,individuals.chrom,sizepop,bound);
    % 变异操作
    individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,sizepop,[i maxgen],bound);
    
    if mod(i,10)==0
        individuals.chrom=nonlinear(individuals.chrom,sizepop);
    end
    
    % 计算适应度
    for j=1:sizepop
        x=individuals.chrom(j,:);
        individuals.fitness(j)=fun(x);
    end
    
    %找到最小和最大适应度的染色体及它们在种群中的位置
    [newbestfitness,newbestindex]=min(individuals.fitness);
    [worestfitness,worestindex]=max(individuals.fitness);
    % 代替上一次进化中最好的染色体
    if bestfitness>newbestfitness
        bestfitness=newbestfitness;
        bestchrom=individuals.chrom(newbestindex,:);
    end
    individuals.chrom(worestindex,:)=bestchrom;
    individuals.fitness(worestindex)=bestfitness;
    
    avgfitness=sum(individuals.fitness)/sizepop;
    
    trace=[trace;avgfitness bestfitness]; %记录每一代进化中最好的适应度和平均适应度
end
%进化结束

结果显示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
%% 结果显示

[r c]=size(trace);

figure

plot([1:r]',trace(:,1),'r-',[1:r]',trace(:,2),'b--');

title(['函数值曲线 ' '终止代数＝' num2str(maxgen)],'fontsize',12);

xlabel('进化代数','fontsize',12);ylabel('函数值','fontsize',12);

legend('各代平均值','各代最佳值');

disp('函数值 变量');

ylim([1.5 8])

%xlim([1,size(trace,1)])

grid on

% 窗口显示
disp([bestfitness x]);

函数值 变量
2.0000 1.5708 1.5708 1.5708 1.5708 1.5708

可知当种群进化到30代时，函数值收敛到2.0000，在x1、x2、x2、x3、x4、x5分别取1.5708、1.5708、1.5708、1.5708、1.5708 时达到该值。

本文标题：
基于遗传算法和非线性规划的函数寻优问题
文章作者：
嘲弄i
原始链接：
许可协议：
署名–非商业性使用–禁止演绎 4.0 国际，转载请保留原文链接及作者信息。
扫码加我微信