遗传算法(GA)学习笔记.md

本文受到樊哲勇老师博客中的一个用matlab实现的50行的遗传算法程序 启发，并在此基础上增加了注释，并将遗传算法中的一些关键步骤重写了函数，使得整个代码结构、逻辑清晰。代码主要有六个文件：

main.m : 主函数，在这里调用遗传算法

my_ga.m : 完整的遗传算法函数

my_gene.m : 基因编码，产生初代基因

my_fitness.m : 适应度函数

my_cross.m : 用于染色体交叉

my_mutations.m : 基于基因突变突变方式

下面我将解释每个文件：

main.m 主函数

clear; close all;



% 调用 my_ga 进行计算

[best_fitness, elite, generation] = my_ga(15, 'my_fitness', 100, 50, 0.3, 1000, 1.0e-6);

%    number_of_variables, ...    % 求解问题的参数个数
%    fitness_function, ...       % 自定义适应度函数名
%    population_size, ...        % 种群规模（每一代个体数目）
%    parent_number, ...          % 每一代中保持不变的数目（除了变异）
%    mutation_rate, ...          % 变异概率
%    maximal_generation, ...     % 最大演化代数
%    minimal_cost ...            % 最小目标值（函数值越小，则适应度越高）

min(best_fitness)      % 最短路径
elite(generation,:)	   % 最优路径

% 最佳适应度的演化情况
figure
plot(best_fitness);
xlabel('Generation','fontsize',15);
ylabel('Best Fitness','fontsize',15);

在主函数中主要调用遗传算法函数，传入的参数主要有参数个数、自定义适应度函数名称、种群规模、种群中父代的个数、变异概率、最大演化代数和最小目标值。

---

my_ga.m 遗传算法函数

function [best_fitness, elite, generation, last_generation] = my_ga( ...
    number_of_variables, ...    % 求解问题的参数个数
    fitness_function, ...       % 自定义适应度函数名
    population_size, ...        % 种群规模（每一代个体数目）
    parent_number, ...          % 每一代中保持不变的数目（除了变异）
    mutation_rate, ...          % 变异概率
    maximal_generation, ...     % 最大演化代数
    minimal_cost ...            % 最小目标值（函数值越小，则适应度越高）
)

% 累加概率
% 假设 parent_number = 10
% 分子 parent_number:-1:1 用于生成一个数列
% 分母 sum(parent_number:-1:1) 是一个求和结果（一个数）
%
% 分子 10     9     8     7     6     5     4     3     2     1
% 分母 55
% 相除 0.1818    0.1636    0.1455    0.1273    0.1091    0.0909    0.0727    0.0545    0.0364    0.0182
% 累加 0.1818    0.3455    0.4909    0.6182    0.7273    0.8182    0.8909    0.9455    0.9818    1.0000
%
% 运算结果可以看出, 累加概率函数是一个从0到1增长得越来越慢的函数
% 因为后面加的概率越来越小（数列是降虚排列的）
cumulative_probabilities = cumsum((parent_number:-1:1) / sum(parent_number:-1:1)); % 1个长度为parent_number的数列,其值为[0,1]，越往后越慢

% 最佳适应度
% 每一代的最佳适应度都先初始化为1
best_fitness = ones(maximal_generation, 1);

% 精英
% 每一代的精英的参数值都先初始化为0
elite = zeros(maximal_generation, number_of_variables);

% 子女数量
% 种群数量 - 父母数量（父母即每一代中不发生改变的个体）
child_number = population_size - parent_number; % 每一代子女的数目

% 初始化种群 --- 产生基因
% population_size 对应矩阵的行，每一行表示1个个体，行数=个体数（种群数量）
% number_of_variables 对应矩阵的列，列数=参数个数（个体特征由这些参数表示） 
% population = rand(population_size, number_of_variables);
population = my_gene(population_size, number_of_variables);

last_generation = 0; % 记录跳出循环时的代数

% 后面的代码都在for循环中
for generation = 1 : maximal_generation % 演化循环开始
    
    % feval把数据带入到一个定义好的函数句柄中计算
    % 把population矩阵带入fitness_function函数计算
    cost = feval(fitness_function, population); % 计算所有个体的适应度（population_size*1的矩阵）

    % index记录排序后每个值原来的行数
    [cost, index] = sort(cost); % 将适应度函数值从小到大排序

    % index(1:parent_number) 
    % 前parent_number个cost较小的个体在种群population中的行数
    % 选出这部分(parent_number个)个体作为父母，其实parent_number对应交叉概率
    population = population(index(1:parent_number), :); % 先保留一部分较优的个体
    % 可以看出population矩阵是不断变化的

    % cost在经过前面的sort排序后，矩阵已经改变为升序的
    % cost(1)即为本代的最佳适应度
    best_fitness(generation) = cost(1); % 记录本代的最佳适应度

    % population矩阵第一行为本代的精英个体
    elite(generation, :) = population(1, :); % 记录本代的最优解（精英）

    % 若本代的最优解已足够好，则停止演化
    if best_fitness(generation) < minimal_cost; 
        last_generation = generation;
        break; 
    end
    % 交叉变异产生新的种群
        
    % 染色体交叉开始
    population = my_cross(population, cumulative_probabilities, child_number, number_of_variables, parent_number);
    % 染色体交叉结束
     
    % 染色体变异开始
    mutation_population = population(2:population_size, :); % 精英不参与变异，所以从2开始
    population(2:population_size, :) = my_mutations(mutation_population, mutation_rate); % 发生变异之后的种群
    % 染色体变异结束
     
    if(mod(generation, 100) == 0)
        disp(['程序已迭代：', num2str(generation) ,'次']);
    end
   
end % 演化循环结束

我直接跳过前面的部分参数的初始化的过程，直接说一些关键的点。

在进入循环之前，先调用my_gene函数生成初代种群，然后进入循环中。在每次循环中，先进行适应度计算，这个时候调用的是函数my_fitness，将计算好的适应度函数值从小到大进行排序，保留前parent_number个体作为父代，后面的个体就淘汰了，并记录本代中最优的适应度和路径(即为第一个数据)。因为淘汰了部分个体，所以要从保留下来的父代中生成子代，并组成下一次迭代的种群，这个时候就调用了函数my_cross，主要通过染色体的交叉来实现，这样一个两个父代就可以产生两个子代。为了避免算法陷入局部最优解，要使部分染色体进行变异，变异调用的函数就my_mutations。这就是遗传算法的主要流程，也就是我们经常听到的优胜劣汰->染色体交叉->染色体变异。

---

下面我将按照上一段所提到的函数顺序，依次来介绍这些函数：

my_gene.m 产生初代种群

function population = my_gene(population_size, number_of_variables)
  % 产生初代基因（基因编码）
  % population_size ： 种群规模
  % number_of_variables ： 变量的个数/基因的长度
  % population = rand(population_size, number_of_variables);
  
  population = rand(population_size, number_of_variables);
  
  for i = 1 : population_size
    population(i,:) = randperm(number_of_variables);
  end

这里以一个简单的TSP问题作为例子，使用的编码就是整数编码，每个整数都代表一个城市。因为整数编码要比01编码复杂的多，所以这里写出整数编码的遗传算法。

可以看到，这里只是使用了一个循环和一个函数rendperm()，这个函数可以生成1~number_of_variables之间所有数字的任意排列，通过这种方式，我们就可以迭代生成初始种群。

---

my_fitness.m 适应度函数计算

function y = my_fitness(population)

% 简答的TSP问题
costM = [0 159 52 186 194 62 149 89 42 101 42 159 40 39 172 ;
159 0 189 162 198 78 113 178 188 38 176 6 119 184 81 ;
52 189 0 14 180 73 30 92 184 166 115 11 171 192 11 ;
186 162 14 0 89 62 74 191 81 62 164 31 117 196 83 ;
194 198 180 89 0 39 199 87 109 143 78 151 3 101 41 ;
62 78 73 62 39 0 123 22 159 4 89 11 11 53 129 ;
149 113 30 74 199 123 0 116 197 37 167 11 197 138 106 ;
89 178 92 191 87 22 116 0 194 144 33 64 19 18 172 ;
42 188 184 81 109 159 197 194 0 184 46 53 185 27 166 ;
101 38 166 62 143 4 37 144 184 0 11 143 33 70 56 ;
42 176 115 164 78 89 167 33 46 11 0 79 47 139 79 ;
159 6 11 31 151 11 11 64 53 143 79 0 114 3 116 ;
40 119 171 117 3 11 197 19 185 33 47 114 0 14 169 ;
39 184 192 196 101 53 138 18 27 70 139 3 14 0 159 ;
172 81 11 83 41 129 106 172 166 56 79 116 169 159 0 ];

y = zeros(size(population,1),1);

for i = 1 : size(population,1)
    for j = 1 : 12
      y(i) = y(i) + costM(population(i,j),population(i,j+1));
    end
    y(i) = y(i) + costM(population(i,13),population(i,1));
end

通常我们把TSP中的OD矩阵放在适应度函数中，这样我们在计算的时候可以直接调用，在TSP问题中，适应度函数就这条路径的总路程，我们进行循环计算就可以求出每条染色体指代的路径的总和。针对于不同的问题，可以设置不同的适应度函数。

---

my_cross.m 染色体交叉

function y = my_cross(population, cumulative_probabilities, child_number, number_of_variables, parent_number)
    % population : 族群
    % cumulative_probabilities : 累计概率
    % child_number : 后代数量
    % number_of_variables : 变量数量
    % parent_number : 父代数量
    
    for child = 1:2:child_number % 步长为2是因为每一次交叉会产生2个孩子
      % cumulative_probabilities长度为parent_number
      % 从中随机选择2个父母出来  (child+parent_number)%parent_number
      mother_index = find(cumulative_probabilities > rand, 1); % 选择一个较优秀的母亲
      father_index = find(cumulative_probabilities > rand, 1); % 选择一个较优秀的父亲
      
      mother = population(mother_index, :);
      father = population(father_index, :);
      
      % 生成两个交换点
      while 1
        crossover_point1 = ceil(rand*number_of_variables); % 随机地确定一个染色体交叉点
        crossover_point2 = ceil(rand*number_of_variables); % 随机地确定一个染色体交叉点
        if(crossover_point1 ~= crossover_point2)
          break;
        end
      end
      
      temp = 0; %临时变量
      
      % 开始交换
      % 按位交换，在本染色体中，将被交换的基因与本染色体内的基因交换即可
      for i = min(crossover_point1, crossover_point2) : max(crossover_point1, crossover_point2)
        % mother交换
        temp = mother(i);
        mother(find(mother == father(i))) = mother(i);
        mother(i) = father(i);
        %father交换
        father(find(father == temp)) = father(i);
        father(i) = temp;
      end
     
      % 得到下一代
      population(parent_number + child, :) = mother; % 一个孩子
      population(parent_number+child+1, :) = father; % 另一个孩子
  end 
    
  y = population;

这里首先利用了一个列加概率函数cumulative_probabilities，这是一个缓慢增长的函数，最大值为1。然后利用find函数取随机值，选出cumulative_probabilities大于这个随机值中的第一个值。这种操作相当于随机选取一个母本或者父本，**但是为什么不直接使用随机函数，然后再取整如ceil(rand * parent_number)，而要使用一个累加概率函数？**留待之后思考。

选取出父本和母本之后就要进行“繁殖”了，“繁殖”的方式就是再随机产生两个交叉点，在交叉点之间的染色体进行交换。**这个时候需要特别注意：**01的染色体可以直接交换，但是在整数染色体中，直接交换会导致染色体中有重复的数字。所以这个时候的交换就采取了一种巧妙的方式，比如父本中基因‘3’要与母本中的基因‘9’进行交换，那么操作就是父本中的‘3’与父本自己的‘9’进行交换，而母本的‘9’也与自己的‘3’进行交换，这样，既达到了交叉的目的，还避免了染色体中有重复的基因。

---

my_mutations.m 染色体变异

function mutation = my_mutations(mutation_population, mutation_rate)
  % 基因发生变异
  % mutation_population : 变异的种群
  % mutation_rate ： 变异率
  
  % 简单的TSP问题
  number_of_mutations = ceil(size(mutation_population,1) * mutation_rate); % 变异的基因数目（基因总数*变异率）
    
  for i = 1 : number_of_mutations
    row = ceil(rand * size(mutation_population,1));
    r1 = ceil(rand * size(mutation_population,2));
    r2 = ceil(rand * size(mutation_population,2));
    temp = mutation_population(row, r1);
    mutation_population(row, r1) = mutation_population(row, r2);
    mutation_population(row, r2) = temp;
  end
  
  %进化逆转
  for i = 1 : number_of_mutations
    row = ceil(rand * size(mutation_population,1));
    r1 = ceil(rand * size(mutation_population,2));
    r2 = ceil(rand * size(mutation_population,2));
    
    mutation_population(row, min(r1, r2) : max(r1, r2)) = mutation_population(row, max(r1, r2) : -1 : min(r1, r2));
    
  end
  
  mutation = mutation_population;

在变异中主要进行两个操作，一个变异、一个是进化逆转。首先我们根据变异率计算出需要变异的染色体个数，然后进入循环，随机得到一个需要变异的染色体，再随机得到两个变异的基因，交换这两个基因就可以实现变异操作。所谓进化逆转指的是对染色体中随机两个基因之间的序列进行倒转，比如之前某两个基因之间是1234，进化逆转之后就变化为4321。其主要目的还是为了避免算法陷入局部最优之中。

---

最后运行的结果示意如下：

main
程序已迭代：100次
程序已迭代：200次
程序已迭代：300次
程序已迭代：400次
程序已迭代：500次
程序已迭代：600次
程序已迭代：700次
程序已迭代：800次
程序已迭代：900次
程序已迭代：1000次

ans =

266

ans =

11    10     6     8    13     5    15     3     7    12    14     9     1     2     4