數學建模方法-遺傳算法(實戰篇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次的時候，整個種群的平均函數最大值就已經接近真正最大值了。