本文用经典遗传算法框架模板,对matlab新手友好,快速上手看懂matlab代码,快速应用实践,源代码在文末给出。

基本原理:

遗传算法(Genetic Algorithm,GA)是一种受生物学启发的优化算法,它模拟了生物进化中的基因遗传和自然选择过程,通过模拟“进化”过程来搜索问题的最优解或近似最优解。

遗传算法的基本思想是通过对候选解(称为个体或染色体)的不断进化和变异来搜索最优解。在遗传算法中,每个个体都代表了问题的一个潜在解,并且具有一个与其适应度相关联的数值,用于评价个体的优劣。

遗传算法的主要步骤包括:

  • 初始化:随机生成一组初始个体(种群),并为每个个体计算适应度。
  • 选择:根据个体的适应度值,使用选择算子(如轮盘赌算法、竞争选择等)选择一部分个体作为父代。
  • 交叉:从父代中选取一对个体,通过交叉操作产生新的个体(子代),以增加种群的多样性。
  • 变异:对新生成的个体进行变异操作,以引入新的基因组合,增加种群的多样性。
  • 替换:使用替换策略(如保留最优个体、精英保留等)更新种群,产生下一代。
  • 重复迭代:重复进行选择、交叉、变异和替换步骤,直到满足停止条件(如达到最大迭代次数、找到满意解等)。

        通过不断重复这些步骤,遗传算法能够在解空间中搜索出较好的解,逐步逼近最优解。由于其并行性和全局搜索能力,遗传算法被广泛应用于各种领域的优化问题,如工程优化、机器学习、数据挖掘等。

参数与子函数定义:

此代码的待解决问题为function函数,也就是多元函数平方和求极小值问题,看懂代码后,可将function应用为自己的待解决问题。

function MySGA
NP = 50;                     %种群规模
Max_N = 1000;                %迭代次数
Pc = ones(1,NP)*0.8;         %ones表示创建一个全为1的数组,参数为数组大小
Pm = ones(1,NP)*0.1;         %变异概率
D = 500;
MinX = -100;
MaxX = 100;
Error = 1.0e-10;
X = MinX + (MaxX - MinX)*rand(NP,D);
F=fun(X);
[BestF,BestLow] = min(F);

function F = fun(X)
for i = 1:1:length(X(:,1))
    F(i) = sum(X(i,:).^2)
end

代码讲解:

ones函数:

这里出现了ones函数

  • 用法为:A = ones(sz)
  • sz为参数,可以是一个表示数组大小的向量,可以是一个数字标量、一个向量或一个矩阵
  • 这个函数表示将一个空间内的元素都置为1
  • ones(3)创建一个 3x3 的矩阵,所有元素都为 1。只有一个参数默认还是2维x2维的矩阵。尽量使用ones(3,3)来创建,增加代码可读性
  • ones(2,4) 创建一个 2x4 的矩阵,所有元素都为 1。
  • ones([2,3,4])创建一个 2x3x4 的数组,所有元素都为 1。

标量,向量,矩阵

这里提到了标量,向量,矩阵,这里继续学习一下:

  • 数字标量:就是一个普通的数
  • 向量:向量是一种特殊的数组,它只有一个维度,有水平和垂直两种向量,也可以说行向量或者列向量,例如:行向量:v = [1, 2, 3, 4];列向量:v = [1; 2; 3; 4];
  • 矩阵:矩阵是一个由数字按行和列排列成的二维数组。
  • 创建一个矩阵:A = [1, 2, 3; 4, 5, 6; 7, 8, 9];

这里的Pc = ones(1,NP)*0.8就表示创建一个包含1行,NP列的数组,全为0.8,表示变异概率。

这里出现了Error = 1.0e-10;,这种表示方法为科学计数法表示1 \times 10^{-10} 

接下来是X = MinX + (MaxX - MinX)*rand(NP,D);这里先看rand函数

rand函数:

  • 定义:rand函数用于生成指定大小的服从均匀分布的随机数矩阵或随机数向量
  • 参数:两个参数表示生成mxn维的矩阵,矩阵元素为0-1之间的随机数,若有一个参数为n则表示生成一个包含n个0-1之间随机数的向量,若无参数则生成一个0-1之间的随机数
  • 返回:例如r = rand(3, 3),这会生成一个3x3的矩阵返回给r,且矩阵的元素为(0,1)之间的随机值

        再看这个代码, (MaxX - MinX)*rand(NP,D)表示生成为(0,MaxX - MinX)直接的随机值的NP行,D列的矩阵,前面再加上MinX就表示生成范围为(MinX,MaxX)之间的NP行,D列的矩阵。

F=fun(X)

子函数代码:

function F = fun(X)
for i = 1:1:length(X(:,1))
    F(i) = sum(X(i,:).^2)
end
 for函数:
for index = start:increment:end 
% 在这里执行循环体代码 
end

 start是开始值,increment是步长,end是结束值,步长为1的话可以省略increment

length函数:

length函数用于返回数组的长度或矩阵的最大维度

例如:

vec = [1, 2, 3, 4, 5];
len = length(vec); % 返回 5,因为 vec 含有 5 个元素

 len为5,因为有5个元素

mat = [1, 2, 3; 4, 5, 6];
len = length(mat); % 返回 3,因为矩阵的列数为 3

len为3,因为矩阵有3列 

X(:,1)用法:

 这表示选中X矩阵的所有行的第一个元素,:表示所有,也就是提取出X的第一列的内容

分析此函数功能

for i = 1:1:length(X(:,1))表示根据X的第一列的长度进行遍历,也就是循环次数为X的行数。

F(i) = sum(X(i,:).^2)这个代码中,X(i,:)表示选中第i行的所有列,也就是提取出第i行元素,.^2表示进行乘方运算,那么这个代码的运行结果就是将第i行的所有元素平方和相加赋值给F的第i个元素。

min函数:

  • 定义:min函数用于找到数组或矩阵中的最小值及其对应的索引(若有多个列则按列索引)
  • 格式:[M, I] = min(A)
  • 参数:参数A为待索引矩阵或向量
  • 返回值:返回值为一个一维矩阵,一行或一列矩阵返回该行或该列的最小值和所在行/列位置,多行矩阵的话,第一个元素为A第一列的最小值,第二个元素为A第一列最小值的索引位置

 主函数实现:

选择:

for gen = 1:1:Max_N %循环迭代次数
    time(gen) = gen;
    tmpF = cumsum(F/sum(F);     %归一化
    disp(tmpF);
    CopyX = X;
    for i = 1:1:NP %通过适应度选择个体
       rnd = rand;%生成0-1之间的数
       for flag = 1:1:NP
          if(rnd  < tmpF(flag))
              break;
          end
          X(i,:) = CopyX(flag,:);
       end    
    end
    
end

        如上述代码所示,首先第一个for循环是循环迭代次数,下面的代码为每一次迭代所需要的操作,其中,选择操作采用了轮盘赌选择法。

        time数组记录当前迭代次数,tmpF是制作好的轮盘。

分析tmpF的实现:

        首先轮盘赌选择法是基于适应度进行选择的,适应度与选中概率成正比,适应度越高那么被选中的概率越大。轮盘的制作分为两步:

1.根据适应度计算概率

设种群规模为NF_{i}表示第i个个体的适应值,那么这个个体被选中的概率为:P_{i}=\frac{F_{i}}{\sum_{j=1}^{N}F_{j}}

代码中F/sum(F)实现此功能。

2.归一化制作轮盘

        要实现概率越高选中越大就要将这个所有个体的概率归一化,那么使用cumsum函数计算每个个体与它前所有元素的累加和,这样就制作好了一个轮盘,实现所有个体被选中的累计概率为1。如下图所示:因此第一个for循环代码表示遍历整个种群对每个个体进行更新,rnd=rand表示模拟轮盘转动。第二个for循环进行个体选择。

        上述提到的是经典遗传算法问题的轮盘制作过程,本模板需要计算的是多元函数求极值问题,具体为多个未知数的平方和最小问题,因此这里用F表示适应度的话,F越小、适应度越小被选中概率越高,因此轮盘制作代码需要改为如下:

tmpF = cumsum((1./F)/sum(1./F));

交叉:

    Pc_rand=rand(1,NP);
    for i=1:2:(NP-1)
        if Pc(i) > Pc_rand(i)
            alfa=rand;
            CpX = X(i,:);
            X(i,:)=alfa*X(i+1,:)+(1-alfa)*X(i,:);
            X(i+1,:)=alfa*CpX+(1-alfa)*X(i+1,:);
        end
    end

 Pc_rand为1行NP列,也就是一行NP个元素,每个元素为0-1随机值的矩阵。使用它判断每两个个体是否进行交叉操作,Pc为1行NP列,元素全为0.8的矩阵,代表进行交叉操作的概率。

开始循环,循环次数为种群的一半,代表一次是否选中两个个体进行交叉操作,循环里进行判断是否进行交叉操作,实现了交叉操作概率为0.8。

交叉操作:生成一个0-1之间的随机数alfa = rand;将第i个个体乘以(1-alfa)加上第i+1个个体乘以alfa更新为第i个个体,第i+1个个体同理。

变异:

    Pm_rand=rand(NP,D);
    for i=1:1:NP
        for j=1:1:D
            if Pm(i) > Pm_rand(i,j)
                X(i,j)=MinX+(MaxX-MinX)*rand;
            end
        end
    end

这个代码就很简洁了,遍历每个个体的每个参数,当0.01比这个值大时(0.01就是变异概率),就随机改变这个参数的值,也就是变异了。 

后续处理:

    X(NP,:)=bestX;%保留最佳个体防止淘汰
    F=fun(X); %评估当前种群
    [BestF,BestLow] = min(F);
    bestX=X(BestLow,:);
    result(gen)=BestF;  %记录结果 

        X(NP,:)=bestX就是将上一次循环后的最优个体保留,防止种群的最后一个个体位置,确保每次迭代的最佳个体被保留在种群中,以防止其在后续的进化过程中被淘汰掉。

整体代码:

function MySGA
NP = 100;                     %种群规模
Max_N = 10000;                %迭代次数
Pc = ones(1,NP)*0.8;         %ones表示创建一个全为1的数组,参数为数组大小
Pm = ones(1,NP)*0.01;         %变异概率
flagc=[0,Max_N]; 
D = 50;
MinX = -100;
MaxX = 100;
Error = 1.0e-10;
X = MinX + (MaxX - MinX)*rand(NP,D);
F=fun(X);
[BestF,BestLow] = min(F);
bestX=X(BestLow,:);
%disp(F);
for gen = 1:1:Max_N %循环迭代次数
    time(gen) = gen;
    tmpF = cumsum((1./F)/sum(1./F));     %归一化
    %disp(tmpF);
    CopyX = X;
    for i = 1:1:NP %通过适应度选择个体
       rnd = rand;%生成0-1之间的数
       for flag = 1:1:NP
          if(rnd  < tmpF(flag))
              break;
          end
          X(i,:) = CopyX(flag,:);
       end    
    end
    %种群交叉
    Pc_rand=rand(1,NP);
    for i=1:2:(NP-1)
        if Pc(i) > Pc_rand(i)
            alfa=rand;
            CpX = X(i,:);
            X(i,:)=alfa*X(i+1,:)+(1-alfa)*X(i,:);
            X(i+1,:)=alfa*CpX+(1-alfa)*X(i+1,:);
        end
    end
    %种群变异
    Pm_rand=rand(NP,D);
    for i=1:1:NP
        for j=1:1:D
            if Pm(i) > Pm_rand(i,j)
                X(i,j)=MinX+(MaxX-MinX)*rand;
            end
        end
    end
    X(NP,:)=bestX;%保留最佳个体防止淘汰
    F=fun(X);
    %------------------------------ 种群评估 -------------------------------
    [BestF,BestLow] = min(F);
    bestX=X(BestLow,:);
    %--------------------------- 记录结果 ----------------------------------
    result(gen)=BestF;
    if (BestF<Error) & (flagc(1)==0)
        flagc(1)=1;flagc(2)=gen;
    end
    if mod(gen,1000)==0
        disp(['代数:',num2str(gen),'----最优:',num2str(BestF)]);
        %disp(X);
    end
end
plot(time,result)
function F = fun(X)
for i = 1:1:length(X(:,1))
    F(i) = sum(X(i,:).^2);
end

此代码在交叉操作中还可以采用别的方式,可以读懂代码后,继续尝试增加多种方式。

Logo

旨在为数千万中国开发者提供一个无缝且高效的云端环境,以支持学习、使用和贡献开源项目。

更多推荐