数学建模方法-遗传算法(实战篇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次的时候，整个种群的平均函数最大值就已经接近真正最大值了。