数学建模方法-遗传算法(实战篇part 1)
一、引言
在上一篇中我们详细介绍了什么是遗传算法，但是光说不练是不行的，因此，在这一篇中，我们将举一个例子，并且利用遗传算法来解决我们的例子。
二、问题
已知：$f(x) = x + 10sin5x + 7cos4x, x \in [0, 9]$
求：函数$f(x)$的最大值
三、一般求解
在MATLAB中输入如下代码：
x = 0: 0.0001: 9; y = x + 10*sin(5*x) + 7*cos(4*x); [maxY, index] = max(y) maxX = x(index)
则会输出以下结果：
maxY =
24.8554
index =
78568
maxX =
7.8567
对此，我们得到$f(x)$在其定义域内的最大坐标值为(7.8567, 24.8554)。
好，接下来，我们利用遗传算法来设计代码，对我们这个问题进行求解，看看会是怎样。
三、遗传算法求解
我们来回顾下上个篇章所讲的，遗传算法的步骤，如下：
1. 种群初始化
2. 计算每个种群的适应度值
3. 选择(Selection)
4. 交叉(Crossover)
5. 变异(Mutation)
6. 重复2-5步直至达到进化次数
我们从第一步开始编写我们的代码，为了让我们的遗传算法的代码具有更好的包容性，即针对不同的题目，我们不想每一次都要大幅度的重写我们的代码，因此，我们希望能够把一些步骤的功能编写成函数，这样针对不同的题目，我们只需要调用我们事先编写好的函数，输入不同的参数和数据即可。好了，废话不多说，我们开始吧。
（1）种群初始化函数popInit()，能根据提供的种群大小和染色体长度产生初始的种群。代码如下：
% [name] -- popInit(种群初始化函数) % [function]-- 构建种群矩阵，其中行数为种群个数，列数为染色体长度（即基因的个数），并随机分配染色体的样式 % [input] -- 1. popSize（种群大小/个数) % -- 2. cLength（染色体长度） % [output] -- popMat（种群矩阵） function popMat = popInit(popSize, cLength) popMat = zeros(popSize, cLength); % 预分配内存 for i = 1: popSize for j = 1: cLength popMat(i, j) = round(rand); % rand产生(0,1)之间的随机数，round()是四舍五入函数 end end clear i; clear j;
（2）计算每个种群的适应度值函数getfitnessValue()，这个函数，针对不同的题目有不同的适应度值，因此，对于不同的题目，这个函数需要更改。在基于本章的题目中，该函数可以实现对种群的适应度值的计算。代码如下：
% [name] -- getfitnessValue(计算种群个体适应度值) % [function]-- 根据不同的题目要求，设计适应度方程计算的算法，本例中，适应度函数为f(x) = x + 10*sin(5*x) + 7*cos(4*x)，x ∈[0, 9],其解码规则为： % x = lower_bound + decimal(chromosome) * (upper_bound - lower_bound) / (2 ^ chromosome_length - 1) % [input] -- popMat (种群矩阵) % [output] -- fitValMat（每个个体的适应度值） function fitnessValueMatrix = getfitnessValue(popMat) [popSize, cLength] = size(popMat); % 获取种群的数目（行）和染色体长度（列） fitnessValueMatrix = zeros(popSize, 1); % 初始化适应度矩阵 X = zeros(popSize, 1); % 预分配内存 lower_bound = 0; % 自变量区间下限 upper_bound = 9; % 自变量区间上限 % 1. 首先先将种群中个体的染色体所对应的十进制求出来 for i = 1: popSize for j = 1: cLength if popMat(i, j) == 1 X(i) = X(i) + 2 ^ (j - 1); end end % 2. 根据十进制值进行解码 X(i) = lower_bound + X(i) * (upper_bound - lower_bound) / (2 ^ cLength - 1); % 3. 获取适应度值 fitnessValueMatrix(i) = X(i) + 10 * sin(5 * X(i)) + 7 * cos(4 * X(i)); end clear i; clear j;
（3）选择函数selection()，可以对种群进行选择，具体算法和代码如下：
% [name] -- selection(选择操作)
% [function]-- 采用轮盘赌的一个形式来进行选择，增加不确定性，这样种群就不容易趋向局部最优
% [input] -- 1. fitnessValueMatrix (适应度值矩阵)
% -- 2. popMat(未选择的种群矩阵)
% -- 3. type -- 1: 保留亲代最优个体
% 0: 不保留亲代最优个体
% [output] -- updatedPopMat（经选择后的种群矩阵）
function updatedPopMat = selection(fitnessValueMatrix, popMat, type)
[popSize, cLength] = size(popMat);
updatedPopMat = zeros(popSize, cLength);
% 剔除适应值为负的种群，使其适应值为0
for i = 1: popSize
if fitnessValueMatrix(i, 1) < 0
fitnessValueMatrix(i, 1) = 0;
end
end
% 轮盘赌算法
P = fitnessValueMatrix / sum(fitnessValueMatrix);
Pc = cumsum(P);
for i = 1: popSize
index = find(Pc >= rand);
updatedPopMat(i, :) = popMat(index(1), :);
end
% 是否要保留亲本适应度值最高的，若是，则type = 1,否则为0
if type
[~, bestIndex] = max(fitnessValueMatrix);
updatedPopMat(popSize, :) = popMat(bestIndex, :); % 将亲代最优染色体放到子代的最后一个个体中
end
clear i;
clear j;
（4）交叉函数crossover()，可以对种群进行交叉，交叉的方式又分为单点交叉和多点交叉，根据输入不同的参数来选择不同的实现方式，具体算法和代码如下：
% [name] -- crossover(交叉操作)
% [function]-- 选定交叉点并进行互换
% [input] -- 1. popMat (未交叉的种群矩阵)
% -- 2. type -- 1: 单点交叉
% -- 2: 多点交叉
% -- 3. crossrate (交叉率) -- 建议值为0.6
% [output] -- updatedPopMat（经交叉后的种群矩阵）
function updatedPopMat = crossover(popMat, type, crossrate)
[popSize, cLength] = size(popMat);
if type == 1
% 单点交叉
for i = 1: 2: popSize
if crossrate >= rand
crossPosition = round(rand() * (cLength - 2) + 2); % 随机获取交叉点，去除0和1
% 对 crossPosition及之后的二进制串进行交换
for j = crossPosition: cLength
temp = popMat(i, j);
popMat(i, j) = popMat(i + 1, j);
popMat(i + 1, j) = temp;
end
end
end
updatedPopMat = popMat;
elseif type == 2
% 多点交叉
for i = 1: 2: popSize
if crossrate >= rand
crossPosition1 = round(rand() * (cLength - 2) + 2); % 第一个交叉点
crossPosition2 = round(rand() * (cLength - 2) + 2); % 第二个交叉点
first = min(crossPosition1, crossPosition2);
last = max(crossPosition1, crossPosition2);
for j = first: last
temp = popMat(i, j);
popMat(i, j) = popMat(i + 1, j);
popMat(i + 1, j) = temp;
end
end
end
updatedPopMat = popMat;
else
h = errordlg('type的类型只能为1(单点交叉)或者2(多点交叉)', '进行交叉时发生错误');
end
clear i;
clear j;
（5）变异函数mutation()，可以对种群进行变异，具体算法和代码如下：
% [name] -- mutation(变异操作) % [function]-- 单点变异：随机选择变异点，将其变为0或1 % [input] -- 1. popMat (未交叉的种群矩阵) % -- 2. mutateRate (交叉率) -- 建议值为0.1 % [output] -- updatedPopMat（经交叉后的种群矩阵） function updatedPopMat = mutation(popMat, mutateRate) [popSize, cLength] = size(popMat); for i = 1: popSize if mutateRate >= rand mutatePosition = round(rand() * (cLength - 1) + 1); % 随机获取交叉点，去除0 % 对mutatePosition点进行变异 popMat(i, mutatePosition) = 1 - popMat(i, mutatePosition); end end updatedPopMat = popMat; clear i; clear j;
（6）另外，针对本章的题目，我们需要将二进制的染色体与题目十进制的自变量x值关联起来，因此，我们需要编写一个函数getDecimalValue()来实现这样的工作。代码如下：
% [name] -- getDecimalValue(根据染色体个体（二进制）算出所对应的x值（十进制）)
% [function]-- 根据不同的题目要求，设计适应度方程计算的算法，本例中，适应度函数为f(x) = x + 10*sin(5*x) + 7*cos(4*x)，x ∈[0, 9],其解码规则为：
% x = lower_bound + decimal(chromosome) * (upper_bound - lower_bound) / (2 ^ chromosome_length - 1)
% [input] -- chromosome (染色体)
% [output] -- decimalValue（每个个体的物理意义值）
function decimalValue = getdecimalValue(chromosome)
decimalValue = 0.0; % 初始化
cLength = size(chromosome, 2); % 获取种群的数目（行）和染色体长度（列）
lower_bound = 0; % 自变量区间下限
upper_bound = 9; % 自变量区间上限
% 1. 首先先将种群中个体的染色体所对应的十进制求出来
for i = 1: cLength
if chromosome(1, i) == 1
decimalValue = decimalValue + 2 ^ (i - 1);
end
end
% 2. 解码
decimalValue = lower_bound + decimalValue * (upper_bound - lower_bound) / (2 ^ cLength - 1);
clear i;
（7）另外，我们还编写了一个画图的函数，用于直观的显示在进化的过程中，种群平均适应度值的变化。代码如下：
% [name] -- plotGraph(画图)
% [function]-- 直观的展示在进化过程中，平均适应度值的趋势变化
% [input] -- 1. generationTime(进化次数)
% 2. fitnessAverageValueMatrix(平均适应度值矩阵)
% [output] -- none
function plotGraph(fitnessAverageValueMatrix, generationTime)
x = 1: 1: generationTime;
y = fitnessAverageValueMatrix;
plot(x, y);
xlabel('迭代次数');
ylabel('平均函数最大值');
title('种群平均函数最大值的进化曲线');
好了，至此，我们就完成了解决本题目需要的函数块。接下来，我们只需要编写主函数main.m，针对本章的题目，对其进行求解即可。代码如下：
%% I. 清空变量
clc
clear all
%% II. 参数初始化
popSize = 100; % 种群大小
cLength = 17; % 染色体的长度
generationTime = 200; % 进化次数
crossRate = 0.6; % 交叉概率
mutateRate = 0.01; % 变异概率
fitnessAverageValueMatrix = zeros(generationTime, 1); % 每代平均适应度值
fitnessBestValueMatrix = zeros(generationTime, 1); % 每代最优适应度值
bestIndividual = zeros(generationTime, 1); % 每代最佳自变量（十进制）
%% III. 初始化种群
popMat = popInit(popSize, cLength);
%% IV. 迭代繁衍获取更好的个体(选择、交叉、变异)
for currentTime = 1: generationTime
% 求适应度值，返回适应度值矩阵
fitnessValueMatrix = getfitnessValue(popMat);
% 记录当前最好的数据
if currentTime == 1
[bestValue, bestIndex] = max(fitnessValueMatrix);
fitnessBestValueMatrix(currentTime) = bestValue;
fitnessAverageValueMatrix(currentTime) = mean(fitnessValueMatrix);
else
[bestValue, bestIndex] = max(fitnessValueMatrix);
fitnessBestValueMatrix(currentTime) = max(fitnessBestValueMatrix(currentTime - 1), bestValue);
fitnessAverageValueMatrix(currentTime) = mean(fitnessValueMatrix);
end
bestChromosome = popMat(bestIndex, :); % 最佳适应度值所对应的个体（二进制）
bestIndividual(currentTime) = getdecimalValue(bestChromosome); % 解码，将二进制转为十进制，得到最佳适应度值所对应的x值
if currentTime ~= generationTime
% 选择
popMat = selection(fitnessValueMatrix, popMat, 1); % 保留亲代最佳个体
% 交叉
popMat = crossover(popMat, 1, crossRate); % 单点交叉
% 变异
popMat = mutation(popMat, mutateRate);
end
end
%% V. 画图并展示结果
% 画图
plotGraph(fitnessAverageValueMatrix, generationTime);
% 展示数据
[fitnessBestValue, index] = max(fitnessBestValueMatrix);
disp 最优适应度
fitnessBestValue
disp 最优个体对应的自变量值
bestIndividual(index)
运行上述程序，可以得到：
最优适应度
fitnessBestValue =
24.8554
最优个体对应的自变量值
ans =
7.8567
将上述结果跟我们第二步用一般求解的对比发现，答案一致。因此，我们的遗传算法是可行的。最后po一张图，可以明显看到，差不多迭代到30次的时候，整个种群的平均函数最大值就已经接近真正最大值了。
