Tianyou Zhao 的个人博客
0%

遗传算法学习笔记

发表于 更新于 分类于 Valine：
本文字数： 8k 阅读时长 ≈ 7 分钟
源自一次大一下学期的数模比赛

遗传算法

演化思想

遗传算法的本质

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

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

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

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

Ras函数

Ras函数（Rastrigin’s Function）是一个典型的非线性==多峰==函数，具有多个局部的极大值点和极小值点，常规算法很容易陷入局部最优解中。
在n维空间中，Ras函数表达式为：
$$ f(x)=An+\sum_{i=1}^{n}\ [x^{2}{i}-Acos(2\pi x{i})] $$
其中，$A$是可以人为定义的常量，$n$是Ras函数的维数
当$n=2$时，Ras函数为二维形式，这个函数在$(0,0)$有全局最小值。

下图是$A=10,n=2$时，Ras函数在$x,y\in[-5,5]$上的三维图象

hello

遗传算法

==目标函数== objective function

1
2
3
4
5
6
7
8
9
%%
% 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$,$x$,$y$。其中$A$是我们自行设定的常数，而$x，y$则对应二维平面上的点
补充：
function------->matlab官方文档-function
fsurf------->matlab官方文档-fsurf

==初始化种群== initialization population

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

==解码== decoding：from binary to decimal

1
2
3
4
5
6
7
8
9
10
11
12
13
14
% 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的时候注意元胞数组的索引方式
    每个个体与向量$vector$相乘，可以将16位二进制数转化成十进制数$$vector=[2^{15} ，2^{14},…,2,1]$$
    得到的乘积是一个范围落在在$(0,65535)$内的整数。为了将这个范围缩小到10附近，再将结果除以6553。
    函数最后输出一个50*2的矩阵，储存50个个体的二维坐标值
    补充：
    mat2cell函数------->matlab官方文档-mat2cell
    ones函数------->matlab官方文档-ones

==计算目标函数== calculate the real value

1
2
3
4
5
6
7
8
9
10
11
% 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

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
% 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 = \frac{x-x_{min}}{x_{max}-x_{min}}$$
如果求解最大值，则$X$越大，个体适应度越高；
如果求解最小值时，则$(1-X)$的结果越大，个体适应度越高。

==选择最优个体== choose the best individual

1
2
3
4
5
6
7
8
% best
function best = bestIndividual(pop)
    popDe = decoding(pop);
    obj = calobjvalue(popDe);
    fit = calfitvalue(obj);
    [r,c] = max(fit);
    best = popDe(r,:);
end

==选择复制== selection copy

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
% 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)$。

==交叉互换== crossover

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
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
)$

==随机变异== mutation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
% 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)$，更加逼近全局最优解。

==总代码==

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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
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

==未回答的问题==

算法虽好，仍有遗憾

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

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