一、引言
在上一篇中我們詳細介紹了什么是遺傳算法，但是光說不練是不行的，因此，在這一篇中，我們將舉一個例子，並且利用遺傳算法來解決我們的例子。
二、問題
已知：$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次的時候，整個種群的平均函數最大值就已經接近真正最大值了。