1、引言
在上一篇中咱們詳細介紹了什麼是遺傳算法，可是光說不練是不行的，所以，在這一篇中，咱們將舉一個例子，而且利用遺傳算法來解決咱們的例子。算法
2、問題
已知：$f(x) = x + 10sin5x + 7cos4x, x \in [0, 9]$函數
求：函數$f(x)$的最大值spa
3、通常求解
在MATLAB中輸入以下代碼：ssr
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)。code
好，接下來，咱們利用遺傳算法來設計代碼，對咱們這個問題進行求解，看看會是怎樣。blog
3、遺傳算法求解
咱們來回顧下上個篇章所講的，遺傳算法的步驟，以下：
內存
1. 種羣初始化
2. 計算每一個種羣的適應度值
3. 選擇(Selection)
4. 交叉(Crossover)
5. 變異(Mutation)
6. 重複2-5步直至達到進化次數ci
咱們從第一步開始編寫咱們的代碼，爲了讓咱們的遺傳算法的代碼具備更好的包容性，即針對不一樣的題目，咱們不想每一次都要大幅度的重寫咱們的代碼，所以，咱們但願可以把一些步驟的功能編寫成函數，這樣針對不一樣的題目，咱們只須要調用咱們事先編寫好的函數，輸入不一樣的參數和數據便可。好了，廢話很少說，咱們開始吧。get
（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次的時候，整個種羣的平均函數最大值就已經接近真正最大值了。