算法实例(2)_遗传算法

最新推荐文章于 2024-08-14 20:58:45 发布

半只小闲鱼

最新推荐文章于 2024-08-14 20:58:45 发布

阅读量1.1k

点赞数 2

分类专栏：算法数学建模文章标签：算法

本文链接：https://blog.csdn.net/qq_52689354/article/details/126250465

版权

算法同时被 2 个专栏收录

19 篇文章 0 订阅

订阅专栏

数学建模

15 篇文章 0 订阅

订阅专栏

遗传算法

在开始之前，我们先来了解下遗传算法中的几个概念

概念1：基因和染色体

在遗传算法中，我们首先需要将要解决的问题映射成一个数学问题，也就是所谓的“数学建模”，那么这个问题的一个可行解即被称为一条“染色体”。一个可行解一般由多个元素构成，那么这每一个元素就被称为染色体上的一个“基因”。

比如说，对于如下函数而言，[1,2,3]、[1,3,2]、[3,2,1]均是这个函数的可行解（代进去成立即为可行解），那么这些可行解在遗传算法中均被称为染色体。

3x+4y+5z<100
这些可行解一共有三个元素构成，那么在遗传算法中，每个元素就被称为组成染色体的一个基因。

概念2：适应度函数

在自然界中，似乎存在着一个上帝，它能够选择出每一代中比较优良的个体，而淘汰一些环境适应度较差的个人。那么在遗传算法中，如何衡量染色体的优劣呢？这就是由适应度函数完成的。适应度函数在遗传算法中扮演者这个“上帝”的角色。

遗传算法在运行的过程中会进行N次迭代，每次迭代都会生成若干条染色体。适应度函数会给本次迭代中生成的所有染色体打个分，来评判这些染色体的适应度，然后将适应度较低的染色体淘汰掉，只保留适应度较高的染色体，从而经过若干次迭代后染色体的质量将越来越优良。

概念3：交叉

遗传算法每一次迭代都会生成N条染色体，在遗传算法中，这每一次迭代就被称为一次“进化”。那么，每次进化新生成的染色体是如何而来的呢？——答案就是“交叉”，你可以把它理解为交配。

交叉的过程需要从上一代的染色体中寻找两条染色体，一条是爸爸，一条是妈妈。然后将这两条染色体的某一个位置切断，并拼接在一起，从而生成一条新的染色体。这条新染色体上即包含了一定数量的爸爸的基因，也包含了一定数量的妈妈的基因。

在这里插入图片描述
那么，如何从上一代染色体中选出爸爸和妈妈的基因呢？这不是随机选择的，一般是通过轮盘赌算法完成。

在每完成一次进化后，都要计算每一条染色体的适应度，然后采用如下公式计算每一条染色体的适应度概率。那么在进行交叉过程时，就需要根据这个概率来选择父母染色体。适应度比较大的染色体被选中的概率就越高。这也就是为什么遗传算法能保留优良基因的原因。

染色体i被选择的概率 = 染色体i的适应度 / 所有染色体的适应度之和

概念4：变异

交叉能保证每次进化留下优良的基因，但它仅仅是对原有的结果集进行选择，基因还是那么几个，只不过交换了他们的组合顺序。这只能保证经过N次进化后，计算结果更接近于局部最优解，而永远没办法达到全局最优解，为了解决这一个问题，我们需要引入变异。

变异很好理解。当我们通过交叉生成了一条新的染色体后，需要在新染色体上随机选择若干个基因，然后随机修改基因的值，从而给现有的染色体引入了新的基因，突破了当前搜索的限制，更有利于算法寻找到全局最优解。

概念5：复制

每次进化中，为了保留上一代优良的染色体，需要将上一代中适应度最高的几条染色体直接原封不动地复制给下一代。

假设每次进化都需生成N条染色体，那么每次进化中，通过交叉方式需要生成N-M条染色体，剩余的M条染色体通过复制上一代适应度最高的M条染色体而来。

遗传算法的流程
通过上述概念，相信遗传算法的大致原理你已经了解，下面我们将这些概念串联起来，介绍遗传算法的执行流程。
在算法初始阶段，它会随机生成一组可行解，也就是第一代染色体。
然后采用适应度函数分别计算每一条染色体的适应程度，并根据适应程度计算每一条染色体在下一次进化中被选中的概率(这个上面已经介绍，这里不再赘述)。
上面都是准备过程，下面正式进入“进化”过程。

通过“交叉”，生成N-M条染色体；
再对交叉后生成的N-M条染色体进行“变异”操作；
然后使用“复制”的方式生成M条染色体；
到此为止，N条染色体生成完毕！紧接着分别计算N条染色体的适应度和下次被选中的概率。

这就是一次进化的过程，紧接着进行新一轮的进化。

究竟需要进化多少次？

每一次进化都会更优，因此理论上进化的次数越多越好，但在实际应用中往往会在结果精确度和执行效率之间寻找一个平衡点，一般有两种方式。

1. 限定进化次数

在一些实际应用中，可以事先统计出进化的次数。比如，你通过大量实验发现：不管输入的数据如何变化，算法在进化N次之后就能够得到最优解，那么你就可以将进化的次数设成N。

然而，实际情况往往没有那么理想，往往不同的输入会导致得到最优解时的迭代次数相差甚远，这是你可以考虑采用第二种方式。

2. 限定允许范围

如果算法要达到全局最优解可能要进过很多很多很多次的进化，这极大影响系统的性能。那么我们就可以在算法的精确度和系统效率之间寻找一个平衡点。我们可以事先设定一个可以接收的结果范围，当算法进行X次进化后，一旦发现了当前的结果已经在误差范围之内了，那么就终止算法。

但这种方式也有个缺点，有些情况下可能稍微进化几次就进入了误差允许范围，但有些情况下需要进化很多很多很多很多次才能进入误差允许范围。这种不确定性导致算法的执行效率不可控。

所以，究竟选择何种方式来控制算法的迭代次数，这需要你根据具体的业务场景合理地选择。这里无法给出普世的方式，需要你自己在真实的实践中找到答案。

采用遗传算法解决负载均衡调度问题
算法都是用来解决实际问题的，到此为止，我想你对遗传是算法已经有了个全面的认识，下面我们就用遗传算法来解决一个实际问题——负载均衡调度问题。

假设有N个任务，需要负载均衡器分配给M个服务器节点去处理。每个任务的任务长度、每台服务器节点(下面简称“节点”)的处理速度已知，请给出一种任务分配方式，使得所有任务的总处理时间最短。

数学建模

拿到这个问题后，我们首先需要将这个实际问题映射成遗传算法的数学模型。

任务长度矩阵（简称：任务矩阵）
我们将所有任务的任务长度用矩阵tasks表示，如：

Tasks={2,4,6,8}
那么，tasks[i]中的i表示任务的编号，而tasks[i]表示任务i的任务长度。

节点处理速度矩阵（简称：节点矩阵）
我们将所有服务器节点的处理速度用矩阵nodes表示，如：

Nodes={2,1}
那么，nodes[j]中的j表示节点的编号，而nodes[j]表示节点j的处理速度。

任务处理时间矩阵

当任务矩阵Tasks和节点矩阵Nodes确定下来之后，那么所有任务分配给所有节点的任务处理时间都可以确定了，我们用矩阵timeMatrix表示，它是一个二维数组：

timeMatrix[i][j]表示将任务i分配给节点j处理所需的时间，它通过如下公式计算：

timeMatrix[i][j] = tasks[i]/nodes[j]

染色体
通过上文我们知道，每次进化都会产生N条染色体，每一条染色体都是当前问题的一个可行解，可行解由多个元素构成，每个元素称为染色体的一个基因。下面我们就用一个染色体矩阵来记录算法每次进化过程中的可行解。

一条染色体的构成如下：

chromosome={1,2,3,4}

一条染色体就是一个一位数组，一位数组的下标表示任务的编号，数组的值表示节点的编号。那么chromosome[i]=j的含义就是：将任务i分配给了节点j。

上面的例子中，任务集合为Tasks={2,4,6,8}，节点集合为Nodes={2,1}，那么染色体chromosome={3,2,1,0}的含义是：

将任务0分配给3号节点
将任务1分配给2号节点
将任务2分配给1号节点
将任务3分配给0号节点

适应度矩阵

通过上文可知，在遗传算法中扮演者“上帝”角色的是适应度函数，它会评判每一条染色体的适应度，并保留适应度高的染色体、淘汰适应度差的染色体。那么在算法实现时，我们需要一个适应度矩阵，记录当前N条染色体的适应度，如下所示：

adaptability={0.6, 2, 3.2, 1.8}

adaptability数组的下标表示染色体的编号，而adaptability[i]则表示编号为i的染色体的适应度。

在负载均衡调度这个实例中，我们将N个任务执行总时长作为适应度评判的标准。当所有任务分配完后，如果总时长较长，那么适应度就越差；而总时长越短，则适应度越高。

选择概率矩阵

通过上文可知，每次进化过程中，都需要根据适应度矩阵计算每一条染色体在下一次进化中被选择的概率，这个矩阵如下所示：

selectionProbability={0.1, 0.4, 0.2, 0.3}

矩阵的下标表示染色体的编号，而矩阵中的值表示该染色体对应的选择概率。其计算公式如下：

selectionProbability[i] = adaptability[i] / 适应度之和

算法一

上述一切知识点铺垫完成之后，接下来我们就可以上代码了，相信Talk is cheap, show you the code!

/**
 * 遗传算法
 * @param iteratorNum 迭代次数
 * @param chromosomeNum 染色体数量
 */
function gaSearch(iteratorNum, chromosomeNum) {
    // 初始化第一代染色体
    var chromosomeMatrix = createGeneration();

    // 迭代繁衍
    for (var itIndex=1; itIndex<iteratorNum; itIndex++) {
        // 计算上一代各条染色体的适应度
        calAdaptability(chromosomeMatrix);

        // 计算自然选择概率
        calSelectionProbability(adaptability);

        // 生成新一代染色体
        chromosomeMatrix = createGeneration(chromosomeMatrix);

    }
}

代码一来，一切都清晰了，似乎不需要过多的解释了。上面是遗传算法最主要的框架，其中的一些细节封装在了一个个子函数中
完整的代码在这位大佬这里：

https://github.com/bz51/GeneticAlgorithm/blob/master/GA.js

结果展示

在这里插入图片描述

算法二

Step1:d定义>

import numpy as np


def schaffer(p):
    '''
    This function has plenty of local minimum, with strong shocks
    global minimum at (0,0) with value 0
    '''
    x1, x2 = p
    x = np.square(x1) + np.square(x2)
    return 0.5 + (np.square(np.sin(x)) - 0.5) / np.square(1 + 0.001 * x)

Step2：执行遗传算法

from sko.GA import GA

ga = GA(func=schaffer, n_dim=2, size_pop=50, max_iter=800, prob_mut=0.001, lb=[-1, -1], ub=[1, 1], precision=1e-7)
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)

演示代码：

import pandas as pd
import matplotlib.pyplot as plt

Y_history = pd.DataFrame(ga.all_history_Y)
fig, ax = plt.subplots(2, 1)
ax[0].plot(Y_history.index, Y_history.values, '.', color='red')
Y_history.min(axis=1).cummin().plot(kind='line')
plt.show()

在这里插入图片描述

TSP遗传算法(旅行推销员问题)

旅行商问题

旅行商问题(traveling salesman problem,TSP)可描述为:已知N个城市之间的相互距离,现有一个商人必须遍访这N个城市,并且每个城市只能访问一次,最后又必须返回出发城市。如何安排他对这些城市的访问次序,使其旅行路线总长度最短。

旅行商问题是一个典型的组合优化问题,其可能的路径数目是与城市数目N呈指数型增长的，一般很难精确地求出其最优解,因而寻找其有效的近似求解算法具有重要的理论意义,特别是当N的数目很大时,用常规的方法求解计算量太大。对庞大的搜索空间中寻求最优解,对于常规方法和现有的计算工具而言,存在着诸多的计算困难。使用遗传算法的搜索能力可以很容易地解决这类问题。另一方面,很多实际应用问题,经过简化处理后,均可化为旅行商问题,因而对旅行商问题求解方法的研究具有重要的应用价值。
TSP问题编码
设D={dij}是由城市i和城市j之间的距离组成的距离矩阵,旅行商问题就是求出一条通过所有城市且每个城市只通过一次的具有最短距离的回路。
在旅行商问题的各种求解方法中,描述旅行路线的方法主要有如下两种:(1)巡回旅行路线经过的连接两个城市的路线的顺序排列;(2)巡回旅行路线所经过的各个城市的顺序排列。大多数求解旅行商问题的遗传算法是以后者为描述方法的,它们都采用所遍历城市的顺序来表示各个个体的编码串,其等位基因为N个整数值或N个记号.
以城市的遍历次序作为遗传算法的编码,目标函数取路径长度。在群体初始化、交叉操作和变异操作中考虑TSP问题的合法性约束条件(即对所有的城市做到不重不漏)

TSP问题的遗传算法设计

1.参数编码和初始群体设定

一般来说遗传算法对解空间的编码大多采用二进制编码形式,但对于TSP-类排序问题，采用对访问城市序列进行排列组合的方法编码，即某个巡回路径的染色体个体是该巡回路径的城市序列。

针对TSP问题,编码规则通常是取N进制编码,即每个基因仅从1到N的整数里面取一个值,

每个个体的长度为N,N为城市总数。

定义一个s行t列的pop矩阵来表示群体, t为城市个数＋1,即N+1,s为样本中个体数目。

针对30个城市的TSP问题, t取值31,即矩阵每一行的前30个元素表示经过的城市编号,最后一个元素表示经过这些城市要走的距离。
参数编码和初始群体设定程序为:

pop=zeros(s,t);
for i=1:s
pop(i,1:t-1)=randperm(t-1);
end

2.适应度函数设计

在TSP的求解中,用距离的总和作为适应度函数,来衡量求解结果是否最优。将pop矩阵中每一行表示经过的距离的最后一个元素作为适应度函数。
两个城市m和n间的距离为:

function [d]=juli(m,n)
x=[87 91 83 71 64 68 83 87 74 71 58 54 51 37 41 2 7 22 25 18 4 13 18 24 25 41 45 44 58 62];
y=[7 38 46 44 60 58 69 76 78 71 69 62 67 84 94 99 64 60 62 54 50 40 40 42 38 26 21 35 35 32];
%m=m+1;
%n=n+1;
d=sqrt((x(m)-x(n))^2+(y(m)-y(n))^2);
根据t的定义,两两城市组合数共有t一2组,则目标函数为:



自适应度函数取目标函数的倒数,即: 



function [pop]=qiujuli(pop)
[s,t]=size(pop);
for i=1:1:s
   dd=0;
   for j=1:1:t-2
      dd=dd+chap10_1calculate(pop(i,j),pop(i,j+1));
   end
%   dd=dd+ga_tsp_juli(pop(i,1),pop(i,t-1));
   pop(i,t)=dd;
end
### 第三步:计算选择算子
选择就是从群体中选择优胜个体、淘汰劣质个体的操作,它是建立在群体中个体适应度评估基础上。仿真中采用最优保存方法，即将群体中适应度最大的c个个体直接替换适应度最小的c个个体。仿真中,取c=15。

function [pop]=select(pop,k)
 
[s,t]=size(pop);
m11=(pop(:,t));
%m11
m11=m11';
mmax=zeros(1,k);
mmin=zeros(1,k);
 
num=1;
while num<k+1
   [a,mmax(num)]=max(m11);
   m11(mmax(num))=0;
   num=num+1;
end
 
num=1;
while num<k+1
   [b,mmin(num)]=min(m11);
   m11(mmin(num))=a;
   num=num+1;
end
 
for i=1:k
   pop(mmax(i),:)=pop(mmin(i),:);
end

第四步:计算交叉算子

交叉算子在遗传算法中起着核心的作用,它是指将个体进行两两配对,并以交叉概率p。将配对的父代个体加以替换重组而生成新个体的操作。仿真中,取p=0.90。
有序交叉法实现的步骤是:

随机选取两个交叉点crosspoint(1)和 crosspoint(2);
两后代先分别按对应位置复制双亲X1和X2匹配段中的两个子串A1和B1;
在对应位置交换双亲匹配段以外的城市,如果交换后,后代X1中的某一城市a与子串A1中的城市重复﹐则将该城市取代子串B1与Al中的城市a具有相同位置的新城市,直到与子串A1中的城市均不重复为止,对后代X2也采用同样方法

function [pop]=cross(pop)
[s,t]=size(pop);
 
pop1=pop;
 
for i=1:2:s
   m=randperm(t-3)+1;
   crosspoint(1)=min(m(1),m(2));
   crosspoint(2)=max(m(1),m(2));
 
%   middle=pop(i,crosspoint(1)+1:crosspoint(2));
%   pop(i,crosspoint(1)+1:crosspoint(2))=pop(i+1,crosspoint(1)+1:crosspoint(2));
%   pop(i+1,crosspoint(1)+1:crosspoint(2))=middle;
   
   for j=1:crosspoint(1)
      while find(pop(i,crosspoint(1)+1:crosspoint(2))==pop(i,j))
         zhi=find(pop(i,crosspoint(1)+1:crosspoint(2))==pop(i,j));
         y=pop(i+1,crosspoint(1)+zhi);
         pop(i,j)=y;
      end
   end
   
   for j=crosspoint(2)+1:t-1
      while find(pop(i,crosspoint(1)+1:crosspoint(2))==pop(i,j))
         zhi=find(pop(i,crosspoint(1)+1:crosspoint(2))==pop(i,j));
         y=pop(i+1,crosspoint(1)+zhi);
         pop(i,j)=y;
   end
end
end
[pop]=chap10_1dis(pop);
 
for i=1:s
    if pop1(i,t)<pop(i,t)
       pop(i,:)=pop1(i,:);
    end
end

第五步:计算变异算子

变异操作是以变异概率pm对群体中个体串某些基因位上的基因值作变动,若变异后子代的适应度值更加优异,则保留子代染色体,否则,仍保留父代染色体。仿真中,取pm=0.20。
这里采用倒置变异法:假设当前个体X为(1374 8 05 9 6 2),如果当前随机概率值小于pm,则随机选择来自同一-个体的两个点mutatepoint(1)和 mutatepoint(2),然后倒置该两点的中间部分,产生新的个体。
例如,假设随机选择个体X的两个点“7”和“9”,则倒置该两个点的中间部分,即将“4805”变为“5084”,产生新的个体X为(1 3 7 5 0 8 4 9 6 2)。

function [pop]=cross(pop)
[s,t]=size(pop);
 
pop1=pop;
 
for i=1:2:s
   m=randperm(t-3)+1;
   crosspoint(1)=min(m(1),m(2));
   crosspoint(2)=max(m(1),m(2));
 
%   middle=pop(i,crosspoint(1)+1:crosspoint(2));
%   pop(i,crosspoint(1)+1:crosspoint(2))=pop(i+1,crosspoint(1)+1:crosspoint(2));
%   pop(i+1,crosspoint(1)+1:crosspoint(2))=middle;
   
   for j=1:crosspoint(1)
      while find(pop(i,crosspoint(1)+1:crosspoint(2))==pop(i,j))
         zhi=find(pop(i,crosspoint(1)+1:crosspoint(2))==pop(i,j));
         y=pop(i+1,crosspoint(1)+zhi);
         pop(i,j)=y;
      end
   end
   
   for j=crosspoint(2)+1:t-1
      while find(pop(i,crosspoint(1)+1:crosspoint(2))==pop(i,j))
         zhi=find(pop(i,crosspoint(1)+1:crosspoint(2))==pop(i,j));
         y=pop(i+1,crosspoint(1)+zhi);
         pop(i,j)=y;
   end
end
end
[pop]=chap10_1dis(pop);
 
for i=1:s
    if pop1(i,t)<pop(i,t)
       pop(i,:)=pop1(i,:);
    end
end

第六步:主函数

在主函数将以上步骤串起来，这里用到了两个数据集，30和48城市的

clear all;
close all;
 
% t=31;     %Number of Cities is t-1
t=49
s=500;    %Number of Samples
 
pc=0.90;
pm=0.20;
 
pop=zeros(s,t);
 
for i=1:s
   pop(i,1:t-1)=randperm(t-1);
end
 
for k=1:1:5000
if mod(k,10)==1
    k
end
 
pop=chap10_1dis(pop);
 
c=15;
pop=chap10_1select(pop,c);
 
p=rand;
if p>=pc
   pop=chap10_1cross(pop);
end
if p>=pm
   pop=chap10_1mutate(pop);
end
 
end
pop
 
min(pop(:,t))
J=pop(:,t);
fi=1./J;
 
[Oderfi,Indexfi]=sort(fi);   % Arranging fi small to bigger
BestS=pop(Indexfi(s),:);     % Let BestS=E(m), m is the Indexfi belong to max(fi)
 
I=BestS;
 
% x=[87 91 83 71 64 68 83 87 74 71 58 54 51 37 41 2 7 22 25 18 4 13 18 24 25 41 45 44 58 62];
% y=[7 38 46 44 60 58 69 76 78 71 69 62 67 84 94 99 64 60 62 54 50 40 40 42 38 26 21 35 35 32];
 
%x=[87 58 91 83 62 71 64 68 83 87 74 71 58 54 51 37 41 2 7 22 25 18 4 13 18 24 25 41 45 44];
%y=[7 35 38 46 32 44 60 58 69 76 78 71 69 62 67 84 94 99 64 60 62 54 50 40 40 42 38 26 21 35];
 
att48=[1 6734 1453
2 2233 10
3 5530 1424
4 401 841
5 3082 1644
6 7608 4458
7 7573 3716
8 7265 1268
9 6898 1885
10 1112 2049
11 5468 2606
12 5989 2873
13 4706 2674
14 4612 2035
15 6347 2683
16 6107 669
17 7611 5184
18 7462 3590
19 7732 4723
20 5900 3561
21 4483 3369
22 6101 1110
23 5199 2182
24 1633 2809
25 4307 2322
26 675 1006
27 7555 4819
28 7541 3981
29 3177 756
30 7352 4506
31 7545 2801
32 3245 3305
33 6426 3173
34 4608 1198
35 23 2216
36 7248 3779
37 7762 4595
38 7392 2244
39 3484 2829
40 6271 2135
41 4985 140
42 1916 1569
43 7280 4899
44 7509 3239
45 10 2676
46 6807 2993
47 5185 3258
48 3023 1942
];
x=att48(1:end,2);
y=att48(1:end,3);
 
for i=1:1:t-1
	x1(i)=x(I(i));
	y1(i)=y(I(i));
end
x1(t)=x(I(1));
y1(t)=y(I(1));
 
figure(1);
plot(x1,y1,'-ob');

运行结果
30个城市TSP：迭代100次，ans = 672.5364

30个城市TSP：迭代1000次，ans = 433.2609
在这里插入图片描述

简化版代码

第1步：定义您的问题。准备点坐标和距离矩阵。
在这里，我随机生成数据作为演示：

import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt

num_points = 50

points_coordinate = np.random.rand(num_points, 2)  # generate coordinate of points
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')


def cal_total_distance(routine):
    '''The objective function. input routine, return total distance.
    cal_total_distance(np.arange(num_points))
    '''
    num_points, = routine.shape
    return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])

第2步：执行GA

from sko.GA import GA_TSP

ga_tsp = GA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=50, max_iter=500, prob_mut=1)
best_points, best_distance = ga_tsp.run()

Step3：绘制结果：

fig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')
ax[1].plot(ga_tsp.generation_best_Y)
plt.show()