• matlab新手快速上手6（引力搜索算法）

            本文根据一个较为简单的matlab引力搜索算法框架详细分析蚁群算法的实现过程，对matlab新手友好，源码在文末给出。

    引力搜索算法简介：

            引力搜索算法是一种启发式优化算法，最初于2009年由伊朗的Esmat Rashedi、Hossein Nezamabadi-pour和Saeid Saryazdi提出。这种算法灵感来源于引力的物理现象，其中个体之间的相互吸引力和排斥力决定了它们的运动轨迹，进而影响到最终的优化结果。

            这个算法的核心思想是模拟物体之间的引力和排斥力，以在解空间中搜索最优解。具体来说，每个解（个体）都被视为具有质量的物体，它们之间的相互作用由引力和排斥力来描述。通过计算每个解受到的引力和排斥力，可以更新它们的位置，以期望获得更优的解。

    引力搜索算法的一般步骤如下：

    1. 初始化：随机生成初始解（个体）的位置。
    2. 计算适应度：计算每个解的适应度，也就是目标函数的值。
    3. 计算引力和排斥力：根据每个解之间的距离和适应度，计算相互之间的引力和排斥力。
    4. 更新位置：根据引力和排斥力的作用，更新每个解的位置。
    5. 重复迭代：重复执行步骤2到步骤4，直到达到终止条件（如达到最大迭代次数）。
    6. 输出结果：输出最优解或者最优解对应的适应度值。

            引力搜索算法的性能取决于参数的选择、种群大小和迭代次数等因素。这种算法适用于解决各种优化问题，包括连续型和离散型优化问题。

    开始编程：

    参数与子函数定义：

    1. %============================== 引力搜索算法 ==============================
    2. function GSA
    3. %--------------------------------- 共性参数 -------------------------------
    4. NP=30;                             %种群规模
    5. D=10;                              %变量个数
    6. Max_N=1000;                        %限定代数
    7. G0=100;                            %引力常数
    8. alpha=20;                          %引力常数
    9. K0=NP;                             %更新常数
    10. K1=1;                              %更新常数
    11. %--------------------------------- 个性参数 -------------------------------
    12. MinX=-30; MaxX=30;
    13. %-------------------------------- 设置随机数 -------------------------------
    14. rand('state',round(sum(100*clock)));
    15. %---------------------------------- 初始化 ---------------------------------
    16. X=MinX+(MaxX-MinX)*rand(NP,D);
    17. V=zeros(NP,D);
    18. %子函数（目标函数）
    19. function fun=ackley(X)
    20. fun=20+exp(1)-20*exp(-0.2*(sum(X.^2)/length(X))^0.5)...
    21.     -exp(sum(cos(2*pi*X))/length(X));

    参数定义：

            与前几章相同，NP代表天体个数，D代表解的维度。这里rand('state',round(sum(100*clock)));这行代码表示设置随机数种子，以确保每次运行程序时生成的随机数序列是不同的。

            初始化，生成X矩阵，矩阵维度为NP行D列，元素值为(MinX,MaxX)之间的随机值。V矩阵为NP行D列的元素全为0的矩阵。

    子函数（目标函数）：

    数学公式：

    f(x_1, x_2, \ldots, x_n) = 20 + e - 20 \cdot e^{-0.2 \sqrt{\frac{1}{n} \sum_{i=1}^{n} x_i^2}} - e^{\frac{1}{n} \sum_{i=1}^{n} \cos(2\pi x_i)}

    函数性质：

            由公式可知，此函数有多个极小值点，但最小值点在原点处，也就是当x(i)都为0时函数最小，为0。 此算法的目标是找到函数的全局最小值点，即找到使函数值最接近0的变量取值。

    主函数：

    1. %--------------------------------- 优化开始 -------------------------------
    2. for gen=1:1:Max_N
    3.     G=G0*exp(-alpha*gen/Max_N);
    4.     K=round(K0+(K1-K0)*gen/Max_N);
    5.     for i=1:1:NP
    6.         F(i)=ackley(X(i,:));
    7.     end
    8.     [bestF,bestNo]=sort(F);
    9.     best=min(F);
    10.     worst=max(F);
    11.     if best
    12.         m=(F-worst-eps)/(best-worst);
    13.     else
    14.         m=ones(1,NP);
    15.     end
    16.     M=m/sum(m);
    17. %-------开始引力搜索-------
    18.     for i=1:1:NP
    19.         for k=1:1:D
    20.             Ft(i,k)=0;
    21.             for j=1:1:K
    22.                 if bestNo(j)~=i
    23.                     tF(i,bestNo(j),k)=G*M(i)*M(bestNo(j))*...
    24.                       (X(bestNo(j),k)-X(i,k))/norm(X(i,:)-X(bestNo(j),:));
    25.                     Ft(i,k)=Ft(i,k)+rand*tF(i,bestNo(j),k);
    26.                 end
    27.             end
    28.         end
    29.     end
    30.     for i=1:1:NP
    31.         a(i,:)=Ft(i,:)/M(i);
    32.     end
    33.     V=rand(NP,D).*V+a;
    34.     X=X+V;
    35.     %----------------------------- 记录结果 -------------------------------
    36.     GlobalMin_itr(gen)=best;
    37.     if mod(gen,100)==0
    38.         disp(['代数：',num2str(gen),'----最优：',num2str(best),...
    39.             '----中值：',num2str(median(F)),'----均值：',...
    40.             num2str(mean(F)),'----方差：',num2str(var(F))]);
    41.     end
    42. end
    43. GlobalMin=best;
    44. GlobalParams=X(bestNo(1),:);
    45. plot([1:Max_N],GlobalMin_itr);
    46. title('收敛曲线');

             第一个for循环表示循环迭代次数。

            G表示引力常数，公式为：G = G_0 \cdot e^{-\frac{alpha \cdot gen}{MaxN}}可知这是一个递减函数，也就是随着gen的增加，G值越小，也就是迭代次数越多，引力越小。

            K表示更新常数，公式为：K = \text{round}\left(K0 + \frac{(K1 - K0) \cdot gen}{MaxN}\right)，round函数表示四舍五入取整，这个公式表示通过线性插值的方式，将K在迭代过程中逐步从K0变为K1。也就是K0是初始更新常数，K1是最终的更新常数。

            接下来的for循环计算每个个体的适应度的值，存储在F对应的元素中。best存储最好的适应度的值，worst存储最差的，当最小值小于最差值时，计算归一化因子m对适应度的值进行归一化，也就是最好的个体对应的m为1，最差的对应为0。

            M为整体归一化过程，M中所有元素加起来为1，这个框架中的M的计算过程有点繁琐，可以直接采用M = ((1./F)/sum(1./F))来计算，效果是一样的，下面开始引力搜索，两个for循环遍历所有个体的所有元素，先让此元素为0。再次根据K更新常数进行遍历其他天体对此天体的引力影响。

            为什么根据K更新常数进行遍历呢，由上面可知，K是随着遍历次数增多而减少，线性地从NP减少到1，也就意味着在循环开始时，遍历所有个体对当前个体的引力，随着循环次数增多，K就能舍去最小的天体引力，也就是适应度最差的个体，在循环快到最后时，将只计算前几个引力强的个体对当前天体的影响。这就是K的作用。

            接下来是代码的核心部分：

            首先先了解一下引力公式：F = G\frac{Mm}{r^{2}}由引力公式可得出引力与M和m质量成正比，与r呈反比，因此下面实现通过引力更新位置的代码：

            if控制自身天体不会受自身引力影响，接下来就是计算当前天体受到前K个最优天体的引力影响后的方向与位置。tF矩阵中存储三个元素，表示第i个天体受到第bestNo(j)个天体在第k个维度的变化。G表示引力常数，G*M(i)*M(bestNo(j))表示对应上面引力公式的GMm，通过适应度表示质量，因此这个代码就实现了两个引力相互影响下的引力，这部分就实现了公式中的G\frac{Mm}{}部分。

            接下来实现r^{2}部分，由于引力与距离也有关系，距离越大引力越小那么继续编写代码

    首先先理解norm函数，再matlab中，norm([3,4])将返回 5，做这个运算： \sqrt{3^2 + 4^2} = 5 。norm函数就是计算得出两个天体的欧几里得距离.

    *(X(bestNo(j),k)-X(i,k))/norm(X(i,:)-X(bestNo(j),:));这里是难点 原代码中使用的是直接除以距离，虽然这样也可以，但是对比引力公式，这样这不便于理解，但是似乎效果更好，我在这里将源代码更改为此形式：

    (X(bestNo(j),k)-X(i,k))/norm(X(i,:)-X(bestNo(j),:))^2;

    来看此公式，/norm(X(i,:)-X(bestNo(j),:))^2;实现了/r^{2}这一部分，难点在于X(bestNo(j),k)-X(i,k)如何理解，为什么要乘以这个值呢？答案是控制引力方向。如下图所示，当不乘上这个值时，这样我们只计算出了具体的引力值大小F，但是我们需要的是将引力F映射到对应维度的力上，为了使F方向不改变，那么对应的F1和F2就要等比例缩放，因此再乘以X(bestNo(j),k)-X(i,k)这个差值就实现了将力分解到对应的维度上。

    如图表示：二维状态下的力：

            继续通过Ft(i,k)更新第i个天体的第k个维度的受力，用当前维度的力加上tF，tF(i, bestNo(j), k)表示是第i个天体再第bestNo(j)个引力影响下第k维的力。再乘以随机值增加多样性，这样就得到了某个天体在前K个天体的引力影响下，在所有维度的的引力大小。

            继续看下面的代码，a(i,:)=Ft(i,:)/M(i)表示第i个天体的加速度，将时间设为单位时间，那么v = \frac{\Delta s}{\Delta t}a = \frac{\Delta v}{\Delta t}，就变为v = \Delta s，a = \Delta v,因此V=rand(NP,D).*V+a;就表示速度的变化，X=X+V;就表示经过距离的变化后的X。这样就实现了在引力作用下，一个单位时间的天体位置更新。后续就是结果处理，绘制图像等过程。

    norm函数：

    在 MATLAB 中，norm函数用于计算向量的范数。它可以计算向量的不同类型的范数，包括：

    1. 二范数（默认）：向量元素的平方和的平方根。
    2. 一范数：向量元素的绝对值之和。
    3. 无穷范数：向量元素的绝对值的最大值。

    语法通常是norm(X)其中X是一个向量。例如,norm([3,4])将返回 5，因为这个向量的二范数是 \sqrt{3^2 + 4^2} = 5 ，norm([3,4,5]) = \sqrt{3^2 + 4^2 + 5^2} = 7.0711

    源代码：

    1. %============================== 引力搜索算法 ==============================
    2.  
    3. %                  一个伊朗人2009年提出的一个非常漂亮的算法
    4.  
    5. %============================== 引力搜索算法 ==============================
    6. function GSA
    7. %--------------------------------- 共性参数 -------------------------------
    8. NP=30;                             %种群规模
    9. D=10;                              %变量个数
    10. Max_N=10000;                        %限定代数
    11. G0=100;                            %引力常数
    12. alpha=20;                          %引力常数
    13. K0=NP;                             %更新常数
    14. K1=1;                              %更新常数
    15. %--------------------------------- 个性参数 -------------------------------
    16. MinX=-30; MaxX=30;
    17. %-------------------------------- 设置随机数 -------------------------------
    18. rand('state',round(sum(100*clock)));
    19. %rng(round(sum(100*clock)));
    20. %---------------------------------- 初始化 ---------------------------------
    21. X=MinX+(MaxX-MinX)*rand(NP,D);
    22. V=zeros(NP,D);
    23. %--------------------------------- 优化开始 -------------------------------
    24. for gen=1:1:Max_N
    25.     G=G0*exp(-alpha*gen/Max_N);
    26.     K=round(K0+(K1-K0)*gen/Max_N);
    27.     for i=1:1:NP
    28.         F(i)=ackley(X(i,:));
    29.     end
    30.     [bestF,bestNo]=sort(F);
    31.     best=min(F);
    32.     worst=max(F);
    33.     if best
    34.         m=(F-worst-eps)/(best-worst);
    35.         %if gen == 10000
    36.         %    disp(X(bestNo,:));
    37.         %end
    38.         
    39.     else
    40.         m=ones(1,NP);
    41.     end
    42.     %M=m/sum(m);
    43.     M = ((1./F)/sum(1./F));
    44.    %-------开始引力搜索-------
    45.     for i=1:1:NP
    46.         for k=1:1:D
    47.             Ft(i,k)=0;
    48.             for j=1:1:K
    49.                 if bestNo(j)~=i
    50.                     tF(i,bestNo(j),k)=G*M(i)*M(bestNo(j))*...
    51.                     (X(bestNo(j),k)-X(i,k))/norm(X(i,:)-X(bestNo(j),:));
    52.                     %tF(i,bestNo(j),k)=G*M(i)*M(bestNo(j))*...
    53.                     %(X(bestNo(j),k)-X(i,k))/norm(X(i,:)-X(bestNo(j),:))^2;
    54.                     Ft(i,k)=Ft(i,k)+rand*tF(i,bestNo(j),k);
    55.                 end
    56.             end
    57.         end
    58.     end
    59.     for i=1:1:NP
    60.         a(i,:)=Ft(i,:)/M(i);
    61.     end
    62.     V=rand(NP,D).*V+a;
    63.     X=X+V;
    64.     %----------------------------- 记录结果 -------------------------------
    65.     GlobalMin_itr(gen)=best;
    66.     if mod(gen,100)==0
    67.         disp(['代数：',num2str(gen),'----最优：',num2str(best),...
    68.             '----中值：',num2str(median(F)),'----均值：',...
    69.             num2str(mean(F)),'----方差：',num2str(var(F))]);
    70.     end
    71. end
    72. GlobalMin=best;
    73. GlobalParams=X(bestNo(1),:);
    74. plot([1:Max_N],GlobalMin_itr);
    75. title('收敛曲线');
    76.  
    77. function fun=ackley(X)
    78. fun=20+exp(1)-20*exp(-0.2*(sum(X.^2)/length(X))^0.5)...
    79.     -exp(sum(cos(2*pi*X))/length(X));

    结语：

            此章节为作者为准备考试所复习，暂时结束，大致的经典优化算法就是这些，后续遇到更好的智能优化算法还会继续更新。

  • 相关阅读:
    spi协议精讲
    Vue项目初始化搭建JDK+Node.js+npm+vue-cli+创建Vue项目
    Python：不要在循环中删除元素
    达梦数据库使用和常见问题
    期末作业C#实现学生宿舍管理系统
    MySQL执行计划
    linux设置ssh免密登录
    MCE | 神经元为胰腺癌细胞提供营养
    Html-背景属性
    Dapr v1.11 版本已发布
  • 原文地址：https://blog.csdn.net/gr1423428723/article/details/138170753