Matlab与遗传算法

文章目录

• 引入
• 1 步骤描述
• 1.1 编码
• • 1.2 解码
  • 1.3 交配
  • 1.4 突变
  • 1.5 倒位
  • 1.6 个体适应度评估
  • 1.7 复制
• 2 示例
• • 2.1 问题1
  • • 2.1.1 代码
    • 2.1.2 求解步骤

引入

  J. Holland教授于1975年提出Simple Genetic Algorithm (SGA) ¹算法，其由编解码、个体适应度评估和遗传运算三大模块构成；遗传运算包括染色体复制、交叉、变异及倒位等。
  参考书目：卓金武等，MATLAB在数学建模中的应用。

1 步骤描述

1.1 编码

  遗传算法的编码 (基因型)有浮点编码和二进制编码两种，这里只对后者进行介绍。设某一参数的取值范围为$(L,U)$，使用长度为$k$的二进制编码表示该参数，则它共有$2^k$种不一样的编码。该参数编码时的对应关系为：

$$\begin{array}{l}000000000000000000=0\to L\\ 000000000000000001=1\to L+\delta \\ 000000000000000010=2\to L+2\delta \\ 000000000000000011=3\to L+3\delta \\ \cdots \cdots \\ 1111111111111111111=2^k-1\to U\end{array}$$其中

$$\delta =\frac{U-L}{2^k-1}.$$

1.2 解码

  解码 (表现型)的目的是为了将不直观的二进制数据还原成十进制。设某一个体的二进制编码为$b_k{b}_{k-1}\cdots b_2{b}_1$，则相应解码公式为：

$$x=L+(\sum _{i=1}^k{b}_i{2}^{i-1})\frac{U-L}{2^k-1}.$$  例如，设有参数$x\in [2,4]$，现用$5$位二进制数对$x$进行编码，可得$2^5=32$条染色体：

$$00000,00001,\cdots ,11111$$  对于任意二进制数据，例如$x_{22}=10101$，解码后对应的$x$值为：

$$x=2+21\times \frac{4-2}{2^5-1}=3.3548.$$

1.3 交配

  交配运算是使用单点或多点进行交叉的算子。首先用随机数产生一个或多个交配点位置，而后两个个体在交配点位置互换部分基因，造成两个子个体。例以下图中，父染色体$S_1=01001011$、$S_2=10010101$交换其后四位基因，获得$S_1^{\prime }=01000101$、$S_2^{\prime }=10011011$两条子染色体。

1.4 突变

  突变运算是使用基本位进行基因突变。为了不在算法迭代后期出现种群过早收敛，对于二进制基因码组成的个体种群，实行基因码的小概率翻转。例如，$S=11001101$第$3$位翻转，即$S=11001101\to 11101101=S^{\prime }$。$S^{\prime }$可看做是$S$的子染色体。

1.5 倒位

  倒位运算是指一个染色体某区段正常排列顺序发生$180^{\circ }$的颠倒，形成染色体内的DNA从新排列。例如对$S$下划线的部分进行倒位：$S=11\underline{0011}01\Rightarrow S^{\prime }=11\underline{1100}01$。

1.6 个体适应度评估

  **遗传算法依照与个体适应度成正比的概率决定种群中各个个体遗传到下一代的机会。**个体适应度大的个体更容易被遗传到下一代。一般，求目标函数最大值的问题能够直接把目标函数做为检测个体适应度大小的函数。

1.7 复制

  复制运算是根据个体适应度大小决定其下代遗传的可能性。设种群中个体总数为$N$，个体$i$的适应度为$f_i$则个体$i$被选取的概率为：

$$P_i=\frac{f_i}{{\sum }_{k=1}^N{f}_k}.$$  当个体复制概率决定后，再产生$[0,1]$区间的均匀随机数来决定哪一个个体参加复制。

2 示例

2.1 问题1

$$\begin{gathered} maxf(x_1,x_2)=21.5+x_1\sin (4\pi x_1)+x_2\sin (20\pi x_2) \\ \text{s.t.}\{\begin{array}{l}-3.0\le x_1\le 12.1;\\ 4.1\le x_2\le 5.8.\end{array} \end{gathered}$$其三维图以下：

2.1.1 代码

  这里给出完整求解代码，具体的每一步将在代码后进行介绍。

function f = test()
    %% 参数设置
    para.num = 10; % 初始个体数量，建议<100
    para.num_loop = 1000; % 迭代次数，建议100~500
    para.len_x1 = 18; % x1编码长度
    para.len_x2 = 15; % x2编码长度
    para.l_x1 = -3.0; % x1下界
    para.l_x2 = 4.1; % x2下界
    para.u_x1 = 12.1; % x1上界
    para.u_x2 = 5.8; % x2上界
    para.pc = 0.25; % 交配几率，建议0.4~0.99
    para.pm = 0.01; % 基因突变几率，建议0.0001~0.2
    para.num_pc = floor(0.25 * para.num); % 染色体交配数量
    
    %% 参数计算
    para.len = para.len_x1 + para.len_x2;
    para.num_pc = para.num_pc - mod(para.num_pc, 2); % 染色体交配数量默认设置为偶数
    f.obj = 0;
    f.para = [0, 0];

    %% 步骤1.1：生成随机个体
    ini.geni = get_init_geni(para.num, para.len);
    
    for loop = 1: para.num_loop
        fprintf('%d-th looping...\n', loop);
       %% 步骤1.2：解码
        ini.decode = get_decode(ini.geni, para);

        %% 步骤2.1：计算适应度值
        ini.eval = zeros(1, para.num);
        for i = 1: para.num
           ini.eval(i) = get_eval(ini.decode(i, 1), ini.decode(i, 2));
        end
        ini.sum_eval = sum(ini.eval);
        temp_max = max(ini.eval);
		
		%% 保存当前最优解
        f.obj = max(temp_max, f.obj);
        f.para = ini.decode(ini.eval == temp_max, :);

        %% 步骤2.2：计算选取几率
        ini.prob = ini.eval / ini.sum_eval;
        ini.sum_prob = cumsum(ini.prob);
        
        %% 步骤3：新种群复制
        temp_rand = rand(1, para.num);
        temp_choose = zeros(1, para.num);
        for i = 1: para.num
            temp_idx = find(ini.sum_prob > temp_rand(i));
            temp_choose(i) = temp_idx(1);
        end
        ini.geni = ini.geni(temp_choose, :);
        
        %% 步骤4.1：肯定交配个体
        temp_rand = rand(1, para.num);
        [~, temp_min_idx] = sort(temp_rand);
        temp_min_idx = temp_min_idx(1 : para.num_pc);
        
        %% 步骤4.2：开始交配
        for i = 1: floor(para.num_pc / 2)
           % 随机产生交配点
           temp_idx = round(rand(1,1) * (para.len - 1)) + 1;
           temp_geni = ini.geni(temp_min_idx(2 * i - 1), temp_idx: para.len);
           ini.geni(temp_min_idx(2 * i - 1), temp_idx: para.len) = ini.geni(temp_min_idx(2 * i), temp_idx: para.len);
           ini.geni(temp_min_idx(2 * i), temp_idx: para.len) = temp_geni;
        end
        
        %% 步骤5：基因突变
        temp_rand = rand(1, para.num * para.len);
        temp_rand = find(temp_rand <= para.pm);
        temp_idx = floor(temp_rand / para.len) + 1;
        temp_idx(temp_idx == para.num + 1) = para.num;
        temp_geni_idx = mod(temp_rand, para.len);
        temp_geni_idx(temp_geni_idx == 0) = para.len;
        
        if size(temp_idx, 2) == 0
            disp(loop);
            continue;
        end
        
        for i = 1: size(temp_idx, 2)
           if ini.geni(temp_idx(i), temp_geni_idx(i)) == 0
              ini.geni(temp_idx(i), temp_geni_idx(i)) = 1;
           else
               ini.geni(temp_idx(i), temp_geni_idx(i)) = 0;
           end
        end
        
    end
end

function f = get_init_geni(num, len)
    %% 生成初始随机个体
    %   num：初始个体数
    %   len：基因长度
    f = zeros(num, len);
    for i = 1: num
       f(i, :) = rand(1, len) < 0.5; 
    end
end

function f = get_decode(geni, para)
    %% 获取解码值
    f = zeros(para.num, 2);
    for i = 1: para.num
        f(i, 1) = decode(geni(i, 1: para.len_x1), para.l_x1, para.u_x1);
        f(i, 2) = decode(geni(i, para.len_x1 + 1: para.len), para.l_x2, para.u_x2);
    end
end

function f = decode(geni, l, u)
    %% 解码每个个体
    m = size(geni, 2);
    sum_bit = 0;
    for i = 1: m
       sum_bit = sum_bit + geni(i) * 2^(i - 1); 
    end
    f = l + sum_bit * (u - l) / (2^m - 1);
end

function f = get_eval(x1, x2)
    %% 适应度函数，对于极大值问题，适应度函数即优化目标函数
    f = 21.5 + x1 * sin(4 * pi * x1) + x2 * sin(20 * pi * x2);
end

  这里给出迭代一千次的结果：

obj: 38.7462
    para: [10.6994 4.8154]

  适应度的变化以下图：

2.1.2 求解步骤

  这里对详细的步骤进行介绍，整体的步骤与第一节中步骤、代码步骤一致 (假设初始个体数为$10$)。

  1）编 码 解 码
  编码的目的是将变量转换为二进制数串。数串的长度$m_j$取决于所要求的精度$p$，即小数点保留位数：

$$2^{m_j-1}<(U-L)\times 10^p\le 2^{m_j}-1$$所以，对于上题的约束条件，取精度为$4$，有：

$$\{\begin{array}{l}-3.0\le x_1\le 12.1\\ [12.1-(-3.0)]\times 10^4=151000\\ 2^{m_1-1}<151000\le 2^{m_1}-1\\ m_1=18\\ 4.1\le x_1\le 5.8\\ [5.8-4.1]\times 10^4=17000\\ 2^{m_2-1}<17000\le 2^{m_2}-1\\ m_1=15\\ m=m_1+m_2=33\end{array}$$即本例中任意染色体数串的长度均为$33$，其中前$18$位表示$x_1$，后$15$位表示$x_2$。解码时，只需明确当前变量的上下界、数串、数串及起始位置。

  2）评 价 个 体 适 应 度
  对于极大值函数，适应度函数为该函数便可；极小值时，只需目标函数总体取负值。所以，每一个个体 (染色体)的适应度，即将解码后的变量带入适应度函数便可。
  2.1）计算每一个个体的适应度$f_1$，须要保留当前最佳适应度即相应参数；
  2.2）计算适应度总和$F$；
  2.3）计算每一个个体被复制的几率：

$$P_k=\frac{f_k}{F}.$$  2.4）计算每一个个体被复制的累加几率：

$$Q_k=\sum _{j=1}^k{P}_k.$$该函数对应matlab中的cumsum函数。

  3）新 种 群 复 制
  3.1）为每一个个体生成属于$0\sim 1$的随机数，假设十个个体的相关系数以下 (图片源自引入中所提书籍)：

  每次选择时，根据当前个体的随机数进行选择，例如随机数$0.300000$，其大于$Q_3$小于$Q_4$，所以$U_4$被选择。显然，复制的样本源自原样本，且新样本中可能有重复个体。这样选择的依据在于适应度大的个体，其几率区间$[Q_k,Q_{k+1}]$也更大，即更易被选中。

  4）新 种 群 交 配
  4.1）肯定交配染色体数量，例如交配几率$P_c=0.25$，则交配数量为2，代码中为了简便，最终的交配数量将强制为偶数；
  4.2）肯定交配对象，首先为每一个个体生成[0, 1]区间随机数，选取符合交配数量的最小数的个体做为交配对象；
  4.3）交配，首先生成一个介于$0\sim (m-1)$的整数，例如产生的数为$1$，则将选中个体的$[1:]$区间的数值进行交互，示例以下 (图片源自引入中所提书籍)：

  5）基 因 突 变
  本例中，基因的总长度为$10\times 33=330$，假设突变几率为$P_m=0.01$：首先为生成330个介于$0\sim 1$的随机数，并进行编号；获取小于$P_m$的随机数的编号，并经过如下方式得到突变基因所在的位置 (图片源自引入中所提书籍)：

  至此，全新的子代产生，并回退至步骤2，直至达到给定迭代次数。

---

1. J. Holland, Adaptation in Natural and Artificial Systems. ↩︎

本文同步分享在 博客“因吉”（CSDN）。
若有侵权，请联系 support@oschina.cn 删除。
本文参与“OSC源创计划”，欢迎正在阅读的你也加入，一块儿分享。