精华内容
下载资源
问答
  • 用于带障碍的路径规划仿真及可视化,引入了geatpy 进化算法作为示例,可引入其它方法对比其目标函数值。
  • 算法学习之遗传算法路径规划python代码实现)

    千次阅读 多人点赞 2020-05-22 15:47:45
    遗传算法学习——使用python路径规划一、引言二、算法伪代码三、算法流程以及代码实现1、地图创建四、代码工程文档 一、引言   机器人控制过程中,路径规划是其中重要的一环。因此本文采用遗传算法对机器人移动...

    一、引言

      机器人控制过程中,路径规划是其中重要的一环。因此本文采用遗传算法对机器人移动时的路径进行优化,最终选取出一条最优路径。采用栅格图法为机器人移动的地图进行建模,并根据遗传算法机理,在python上仿真实现了最优路径的选取。
      在查阅资料时发现遗传算法做路径规划仿真大多采用matlab软件,而python实现的例子几乎没有。再加上本人正好在学习python,就参考matlab代码做了python的仿真。
    本文重点参考文献链接:基于遗传算法的移动机器人路径规划,在此感谢链接作者。

    二、算法伪代码

      算法伪代码如下:
        step1:建立地图
        step2:初始化种群
        step3:计算个体的适应度值
        step4:选择适应度合适的个体进入下一代
        step5:交叉
        step6:变异
        step7:更新种群,若出现最优路径或达到迭代次数,则转至step8,否则转step3
        step8:输出最优的个体作为最优解
      本文会按照算法的伪代码一步一步讲解。

    三、算法流程以及代码实现

    1、地图创建

      本文采取栅格图法创建地图空间,便于确定直角坐标系。
      如图1所示,图中共有25个栅格,每个栅格表示一段路径,其中黑色方块表示的是障碍物,白色方块表示的是自由栅格,即可行走。
      以图中为例,在构建栅格地图时,首先以地图左下角第一个栅格为坐标原点建立直角坐标系。因此每一个栅格可以用(x,y)的坐标形式表示,比如左下角第一个栅格可以表示为(1,1)。并且从左下角开始每一个栅格进行编号N,编号从0开始,设定编号0为起点,24为终点。说明下,这种编号方式就是遗传算法中的编码操作,会使相对复杂的问题结构,变得简单。比如现在以坐标形式表示一条从(1,1)到(5,5)的路径:(1,1)->(2,2)->(3,2)->(4,2)->(4,3)->(5,4)->(5,5)。这种方法在编程中需要使用二维数组来存放这么一条路径,但通过编码之后,这样的路径可以简化为:(0,6,7,8,13,19,24),只需要一个向量便可以表示。其中坐标和编码之间的转换关系为:
    x=N/col+1x = N/col+1
    y=N%col+1y = N\%col+1

    其中col为栅格图的列数,N为编码,x,y分别为横纵坐标,也可读作x行y列(需要注意,行数从下往上数,列数从左往右数)。
    图1 栅格图编码
      以上方法以一种通用的方法,本案例代码中由于使用python编程,构建地图用的是列表的方式,因此行数从上往下,列数从左往右,并且由于列表的索引从0开始,因此坐标原点也是用的(0,0)如下图所示,光标选中的左上角元素坐标为(0,0),编码为0,编码从左往右一次是0,1,2…直到右下角的99。
    带有障碍的栅格地图
      python实现画栅格地图的代码如下:

    import random
    class Map():
        '''
        :param:地图类
        :param:使用时需要传入行列两个参数 再实例化
        '''
    
        def __init__(self,row,col):
            '''
            :param:row::行
            :param:col::列
            '''
            self.data = []
            self.row = row
            self.col = col
        def map_init(self):
            self.data = [[0 for i in range(self.col)] for j in range(self.row)]
            # for i in range(self.row):
            #     for j in range(self.col):
            #         print(self.data[i][j],end=' ')
            #     print('')
        def map_Obstacle(self,num):
            '''
            :param:num:地图障碍物数量
            :return:返回包含障碍物的地图数据
            '''
            self.num = num
            for i in range(self.num):#生成小于等于num个障碍
                self.data[random.randint(0,self.row-1)][random.randint(0,self.col-1)] = 1
            if self.data[0][0] == 1:        #判断顶点位置是否是障碍  若是 修改成可通行
                self.data[0][0] = 0
    
            if self.data[self.row-1][0] == 1:
                self.data[self.row-1][0] = 0
    
            if self.data[0][self.col-1] == 1:
                self.data[0][self.col - 1] = 0
    
            if self.data[self.row-1][self.col - 1] == 1:
                self.data[self.row - 1][self.col - 1] = 0
            for i in range(self.row):       #显示出来
                for j in range(self.col):
                    print(self.data[i][j],end=' ')
                print('')
            return self.data
    

      使用时代码如下:

    map = Map(10,10)#10行 10列
    map.map_init()#生成初始的10行10列 的0矩阵
    arr = map.map_Obstacle(20)#这里传入的参数是在10*10 的0矩阵上随机生成20个障碍物,因为是随机生成,实际的障碍物数量<=20 返回带有障碍的数据给arr 
    

    到这里,已经生成了一个带有障碍的栅格地图,现在进行step2

    2、种群初始化

      种群的概念是个体的集合。在这之前需说明遗传算法是从问题解的编码组开始,而非单个解开始搜索。即种群的每一个个体都是当前问题的一个可行解,使用遗传算法的过程,是从这些可行解(个体)中找出最优解(适应度最强的个体)。对于本文来讲,种群的初始化就是找到若干组可以从起点到达终点的路径,每一组路径就是一个个体,种群的大小就是个体的数量。遗传算法求解最优路径过程中,种群的初始化是一个难点,也是最难的一步,既要保证路径的可行性,还要保证路径的连续性(即不能跳跃行走)。在遗传算法路径规划中,种群的初始化过程是最重要的,也是最难的。下面开始介绍本文用到的方法。
      本文实现种群初始化分为两步,第一步为找到一条间断的路径,即从栅格图的每一行中随机取出一个不是障碍物的元素。比如在图1中取(0,6,10,16,24),在每一行中都取了一个元素,并且这些元素都是自由栅格。第二步为将间断的路径连续化,在这步中,需要先从路径的第一个栅格判断此条路径与之相邻的下一个栅格是否为连续,判断公式为:
    D=Max{abs(x1x2),abs(y1y2)}D = Max\{abs(x1-x2),abs(y1-y2)\}
    若D=1,则说明路径连续,反之不连续。对于不连续的两个坐标,采用中点法计算两个栅格的中点栅格,公式如下:
    中点法插点
    其中(X_new,Y_new)是新插入栅格的坐标。
     1.若新栅格为障碍物栅格:则以上下左右(这个顺序可以任意,也可以改为8个方向搜索)顺序取新栅格的相邻栅格,并判断此栅格是否已经在路径中。
      a).若此栅格为无障碍栅格且不在路径中则插入路径中;
      b).若遍历上下左右四个栅格后没有满足条件的栅格则删除这条路径;
     2.若新栅格为无障碍物栅格,则插入两个到不连续栅格中间。
     3.继续判断新插入的栅格和新插入的栅格的前一个栅格是否连续,若不连续则循环以上步骤,直到两个栅格连续。当两个栅格连续后取下一个栅格,循环以上步骤,直到整条路径连续。
     4.按照以上步骤生成了一条初始路径之后,会出现走回头路的现象,比如按照图1生成了以下一条路径:(0,6,1,2,6,7,8,13,19,24)可以看出,6号这个节点一共走了两次,这样的话这条路径就会有冗余部分,也就是我认为的走了回头路。这种情况下,为了使算法加快收敛速度,就要保证初始种群比较优良,也就是初始的路径长度就短,那么就有必要滤除这段重复的路段,这里专门设计了一个滤除一段路径中两个相同元素之间内容的算法。代码在我的另一篇博文里面,这里附上链接:python实现将列表中重复元素之间的内容全部滤除
      这里初始化种群的步骤十分繁琐,因此直接上初始化代码,代码上有详细注释,代码如下:

    #step2 初始化种群 也是最难的一步
    class Population():
        def __init__(self,row,col,num,NP):
            self.row = row
            self.col = col
            self.num = num
            self.NP = NP
            self.p_start = 0  # 起始点序号  10进制编码
            self.p_end = self.row * self.col - 1  # 终点序号 10进制编码
            self.xs = (self.p_start) // (self.col)  # 行
            self.ys = (self.p_start) % (self.col)  # 列
            self.xe = (self.p_end) // (self.col)  # 终点所在的行
            self.ye = (self.p_end) % (self.col)
            self.map_start = Map(self.row, self.col)  # Map的实例化 主函数中可以不再初始化地图
            self.map_start.map_init()  # map初始化 生成0矩阵
            self.map = self.map_start.map_Obstacle(self.num)  # 生成带障碍的map
    
            self.can = []  # 这个列表存放整个地图中不是障碍物的点的坐标 按行搜索
            self.popu = [[0 for i in range(self.col)]
                         for j in range(self.NP)]  # 种群的列表,包含NP个间断可行解,即从起点到终点的路径
            # self.popu = []#存放可行解  可行路径
            self.end_popu = []	#存放最终的连续NP条路径
    
        def Population_Init(self):
            '''
            :return:返回NP条连续可行解  即初始路径
            '''
            for i in range(self.NP):       #添加NP条路径
                j = 0
                for xk in range(0,self.row):  #多少行
                    self.can = []             #清空can列表   用来缓存当前行的可行点
                    for yk in range(0,self.col):   #从起点0 开始的列开始搜寻到终点的列  也就是此行中的自由点
                        num = (yk) + (xk) * self.col
                        if self.map_start.data[xk][yk] == 0:
                            self.can.append(num)
                    # print(self.can,'自由点\n')
                    length = len(self.can)     #此行的长度  随机生成指针取出坐标
                    # print(length)
                    self.popu[i][j] = (self.can[random.randint(0,length-1)])
                    j += 1#j从1 开始  是因为 0 的元素时起点编号
                self.popu[i][0] = self.p_start  # 每一行的第一个元素设置为起点序号
                self.popu[i][-1] = self.p_end
    
                temp = self.Generate_Continuous_Path(self.popu[i])#将返回的一条路经 添加进end中
                if temp != []:      #将返回的空列表路径滤除
                    temp = fiter.fiter(temp)    #滤除重复节点路径
                    self.end_popu.append(temp)
                # print(self.end_popu, end='\n')
            # print('测试1',self.popu,end='\n')
            return self.end_popu
        # @staticmethod
        def Generate_Continuous_Path(self,old_popu):#生成连续的路径个体
            '''
            :param old_popu: 未进行连续化的一条路径
            :return:        返回已经连续化的一条路径
            '''
            self.new_popu = old_popu    #传进来的参数就是一行的数组  一条路径
            self.flag = 0
            self.lengh = len(self.new_popu)  #先将第一条路经的长度取出来
            i = 0
            # print("lengh =",self.lengh )
            while i!= self.lengh-1:       #i 不等于 当前行的长度减去1  从0计数  这里有问题 待修改
                x_now = (self.new_popu[i]) // (self.col)  # 行 解码  算出此条路经的第i个元素的直角坐标
                y_now = (self.new_popu[i]) % (self.col)  # 列
                x_next =  (self.new_popu[i+1]) // (self.col) #计算此条路经中下一个点的坐标
                y_next =  (self.new_popu[i+1]) % (self.col)
                #最大迭代次数
                max_iteration = 0
    
                #判断下一个点与当前点的坐标是否连续 等于1 为连续
                while max(abs(x_next - x_now), abs(y_next - y_now)) != 1:
                    x_insert = math.ceil((x_next + x_now) // 2)      #行
                    y_insert = math.ceil((y_next + y_now) // 2) #ceil向上取整数     #列
                    # print("x_insert = ",x_insert,"\ty_insert = ",y_insert)
                    flag1 = 0
    
                    if self.map_start.data[x_insert][y_insert] == 0:  #插入的这个坐标为0 可走
                        num_insert = (y_insert) + (x_insert) * self.col #计算出插入坐标的编码
                        self.new_popu.insert(i+1,num_insert)
                        # print(self.new_popu)
                        # print(num_insert)
                    else:#插入的栅格为障碍  判断插入的栅格上下左右是否为障碍,以及是否在路径中,若不是障碍且不在路径中,就插入
                        #判断下方
                        if (x_insert + 1 < self.row)and flag1 == 0:       #保证坐标是在地图上的
                            if ((self.map_start.data[x_insert+1][y_insert] == 0)#下方不是障碍物
                                and (((y_insert) + (x_insert+1) * self.col) not in self.new_popu)):#编号不在已知路径中
                                num_insert = (y_insert) + (x_insert+1) * self.col  #计算下方的编号
                                self.new_popu.insert(i + 1, num_insert) #插入编号
                                flag1 = 1       #设置标志位 避免下面重复插入
    
                                # print('下方插入',num_insert)
                        #判断右方
                        if (y_insert + 1 < self.col)and flag1 == 0:  # 保证坐标是在地图上的 并且前面没有插入
                            if ((self.map_start.data[x_insert][y_insert+1] == 0)#右方不是障碍物
                                and (((y_insert+1) + (x_insert) * self.col) not in self.new_popu)):#编号不在已知路径中
                                num_insert = (y_insert+1) + (x_insert) * self.col  #计算右方的编号
                                self.new_popu.insert(i + 1, num_insert) #插入编号
                                flag1 = 1  # 设置标志位 避免下面重复插入
                                # print('右方插入',num_insert)
                        #判断上方
                        if (x_insert - 1 > 0) and flag1 == 0:  # 保证坐标是在地图上的
                            if ((self.map_start.data[x_insert-1][y_insert] == 0)#右方不是障碍物
                                and (((y_insert) + (x_insert-1) * self.col) not in self.new_popu)):#编号不在已知路径中
                                num_insert = (y_insert) + (x_insert-1) * self.col  #计算右方的编号
                                self.new_popu.insert(i + 1, num_insert) #插入编号
                                flag1 = 1  # 设置标志位 避免下面重复插入
                                # print('上方插入',num_insert)
                        #判断左方
                        if (y_insert - 1 > 0)and flag1 == 0:  # 保证坐标是在地图上的
                            if ((self.map_start.data[x_insert][y_insert-1] == 0)#右方不是障碍物
                                and (((y_insert-1) + (x_insert) * self.col) not in self.new_popu)):#编号不在已知路径中
                                num_insert = (y_insert-1) + (x_insert) * self.col  #计算右方的编号
                                self.new_popu.insert(i + 1, num_insert) #插入编号
                                flag1 = 1  # 设置标志位 避免下面重复插入
                                # print('左方插入',num_insert)
                        if flag1 == 0:  #如果前面没有插入新点  说明这条路径不对 删除
                            self.new_popu = []
                            break
                    x_next = num_insert//self.col
                    y_next = num_insert%self.col
                    # x_next = x_insert
                    # y_next = y_insert
                    max_iteration += 1#迭代次数+1
                    if max_iteration > 20:
                        self.new_popu = []  #超出迭代次数 说明此条路经可能无法进行连续   删除路径
                        break
                if self.new_popu == []:
                    break
                self.lengh = len(self.new_popu)
                i = i+1
            # print(self.new_popu,'连续')
            return  self.new_popu#返回的是一条路径
    

      这里需要注意的是,在种群初始化对象的__init__方法中已经初始化了地图数据,因此主函数只需要对种群初始化对象实例化就行,使用方法如下:

    population = Population(10,10,20,100)   #实例化对象 并生成初始带随机障碍的地图
    popu = population.Population_Init()  #种群初始化  得到np条初始可行路径
    

    小结

      种群初始化是最难的一步,因为需要判断的地方非常多,比如连续化路径时,还要判断是否超出了地图边界,因为地图是列表数据,一旦超出索引,程序就会报错。还有许多需要注意的地方,这些在程序中都有注释,此外前文提到如果插入的点是障碍物,那么就在其上下左右四个方向搜寻自由点,其实也可以改为在其周围8个方向搜寻自由点,毕竟对角线的方向也是可以走的。
      然而出现的问题还是有的,使用中点插值法使路径连续化会产生一些错误,有的时候会穿过障碍行驶。针对这个问题,在下一步计算适应度函数的过程中,使用了惩罚函数来惩罚这些穿越障碍物的路径,但这个问题是这种方法本身的局限,靠惩罚函数是不能彻底解决的。因此在一些文献中,初始化种群还会采用RRT随机树的方法,这个在之后会改进。
      到这里初始化种群结束,进行下一步,计算适应度函数。

    3、适者生存之适应度函数

      这里就进入到了达尔文进化论的范畴了,其实并没有那么玄乎.适应度函数就是衡量一个问题的解是否是最优的标准,比如一个适应度函数:
    f(x)=x2+1f(x) = x^2+1
      要求f(x)的值取得最小时,x的值是最优的,那么f(x)就是适应度函数的值,x取得每一个值都是适应度函数的解,问题就是这个解是不是最优的。归根到底,遗传算法其实就是解决这样一个问题的方法。
      而本文选取的适应度函数标准为使得路径最短。代码中给定相邻的栅格图为1个单位长度,则对角的长度为1.4(根号2)。其实除了路径最短这一标准外,还有路径平滑度,机器人能耗等标准,由于能力有限,故本文只采用路径最短作为标准。
      根据以上的解释,求适应度函数值就可以理解为求一条路径从起点到终点的距离。代码如下:

    #step3 计算适应度函数
    def calvalue(popu,col):
        '''
        :param popu: 传入种群信息
        :param col: 这个参数是地图的列数 ,因为要解码用到
        :return: 返回的是每一条路径长度组成的列表
        '''
        hang = len(popu)#计算行
        value = [] #存放计算出来的路径长度值
        for i in range(hang):#从第0行开始 到最后一行
            value.append(0)
            single_popu = popu[i] #把当前行的元素赋给 single_popu 之后操作
            single_lengh = len(single_popu)#将当前行的长度拿出来
            for j in range(single_lengh-1):#从第一个元素计算到倒数第一个元素
                x_now = (single_popu[j]) // (col)  # 行 解码  算出此条路经的第i个元素的直角坐标
                y_now = (single_popu[j]) % (col)  # 列
                x_next = (single_popu[j + 1]) // (col)  # 计算此条路经中下一个点的坐标
                y_next = (single_popu[j + 1]) % (col)
                if abs(x_now - x_next) + abs(y_now - y_next) == 1:#路径上下左右连续 不是对角的 则路径长度为1
                    value[i] = value[i]+1
                elif max(abs(x_now - x_next),abs(y_now - y_next))>=2:#惩罚函数 若跳跃或者穿过障碍
                    value[i] = value[i] + 100
                else:
                    value[i] = value[i]+1.4  #对角长度为根号2  即1.4
        return value
    

    小结

      在计算适应度值的时候,传入的种群路径有的其实并不连续,因为跨过了障碍,因此在函数中用以下代码进行惩罚:

    elif max(abs(x_now - x_next),abs(y_now - y_next))>=2:#惩罚函数 若跳跃或者穿过障碍
    	value[i] = value[i] + 100
    

    就是将路径长度直接加长,这样的话在下一步"选择"操作中就大概率的会将这样的路径过滤掉,但不能保证全部过滤。

    4、物竞天择之选择

      这里就是根据适应度值进行选择比较优良的个体存活到下一代,也就是物竞天择。本文中采取了两种方法进行选择优良个体进入下一代。第一种方法是参考文献中所采用的轮盘赌方法;第二种方法是将适应度值从小到大排序,假设种群中强者的概率是50%,即0.5,并且认为强者都能存活(这里强者的概率其实就是对强者筛选过了),即都能进入下一代。那么剩下的50%就是弱者,而对于弱者,也有一个弱者存活概率比如0.3,之后采取随机存活的方式,产生一个0-1之间的随机小数,若0.3大于这个小数 ,那么认为可以存活,否则直接舍弃。存活下来的强者和弱者构成了一个新的种群,进行下一步的交叉操作。首先来看参考文献中的轮盘赌方法进行筛选下一代,代码如下:

    #step4 选择
    def selection(pop,value):
        '''
        :param pop:种群
        :param value:适应度值列表
        :return:返回新的种群
        '''
    
        ###原来的方法会丢失种群数量
        now_value=[]#做倒数后的适应度值列表
        P_value = []  #存放适应度值占总比概率的列表
        random_deci = []#存放随机的小数
        new_popu = []       #选择后的种群
        sum_value = 0   #存放适应度的总和  计算轮盘赌的概率
        lengh = len(pop)#首先计算种群有多少行,也就是有多少初始路径
        for i in range(lengh):
            new_popu.append([]) #原始种群有多少个个体  新的种群就有多少
        for i in value:     #由于适应度值是距离  将需要的就是距离最短的  因此求倒数来排序
            now_value.append(1/i)  #将适应度取倒数   倒数越小适应度越小,路径越长  就可以按概率抛弃
            sum_value += (1/i)
        for i in now_value:#将每一个适应度值取出来  做比例  得出每个适应度值的占比概率 存到列表中
            P_value.append(i/sum_value)
        P_value = np.cumsum(P_value)#累加和 并排序 从小到大
        for i in range(lengh):
            random_deci.append(random.random())#产生 i个随即小数 存到列表中
        random_deci = sorted(random_deci)#从小到大排序
        fitin = 0
        newin = 0
        while(newin<lengh): #遍历每一行
            if random_deci[newin] < P_value[fitin]:
                new_popu[newin] = pop[fitin]#这里的代码有点难理解
                newin += 1
            else:
                fitin += 1
        return new_popu
    

    再看第二种方法,代码如下:

    #step4 选择
    def selection(pop,value):
        '''
        :param pop:种群
        :param value:适应度值列表
        :return:返回新的种群
        '''
        ####原来的方法会丢失种群数量
        retain_rate = 0.6#适应度最优的保留几率
        random_rate = 0.3#适应度差的保留几率
        dict_my = dict(zip(value, pop))     #将路径距离和路径打包字典
    
        new_popu = []
        sort_dis = sorted(dict_my)   #将字典按照键排序  也就是按照距离短-长排序
        retain_lengh = int(len(pop)*retain_rate)      #适应度强的按照存活概率 保留下
        temp = sort_dis[:retain_lengh]      #缓存下保留的距离,待会将距离对应的字典值取出来
        # print(temp,'优秀保留')
        for i in sort_dis[retain_lengh:]:   #距离长的按照随机概率保留
            if random.random() < random_rate:
                temp.append(i)
        #temp现在存放的就是进入下一代的种群信息,他是距离信息,下一步通过字典将种群提取出来
        for i in temp:
            new_popu.append(dict_my[i]) #至此种群选取完
        return new_popu
    

    小结

      这两种方法都能够选取下一代的种群,但都会出现一个问题:在选择的过程中,种群的规模会慢慢变小,有时侯不等迭代结束,种群就所剩无几(即可行路径没剩几条了)虽然能够选出比较优的路径解,但总感觉哪里不对以及路径不是最优的,因此后续还要再查资料,解决这个问题。

    5、遗传学之交叉

      交叉的方法有许多种,由于这里采用十进制编码方式,因此单点交叉方法比较简单实用。
      代码中的思想为:给定一个交叉概率pc,然后再产生0-1之间的一个随机数,并和交叉概率pc比较,若pc大于这个随机小数,则进行交叉操作。
      交叉操作具体为:由于路径的不同而导致个体长度的不同,因此,选择相邻的两个个体,其交叉点的位置不一定相同,这里选取相同编码处进行交叉(起点和终点除外),以保证路径的连续性。具体的交叉操作是找出两条路径中所有相同的点,然后随机选择其中的一个点,将之后的路径进行交叉操作。比如一条连续路径为(0,6,7,8,13,18,23,24),另一条连续路径为(0,1,6,7,8,13,19,24),重复的序号(6,7,8,13),则随机选取一个序号,如13,将13后面的序号交叉,即可得到两条新的路径(0,6,7,8,13,19,24)和(0,1,6,7,8,13,18,23,24),其中(0,6,7,8,13,19,24)其实是肉眼看到的最优的路径,这就是交叉的操作。
    代码如下:

    #step5 交叉   拟采用单点交叉
    def cross(parents,pc):
        '''
        :param parents: 交叉的父类
        :param pc:   交叉概率
        :return:
        '''
        children = []  #首先创建一个子代 空列表 存放交叉后的种群
        single_popu_index_list = []#存放重复内容的指针
        lenparents = len(parents)  #先提取出父代的个数  因为要配对 奇数个的话会剩余一个
        parity = lenparents % 2 #计算出长度奇偶性  parity= 1 说明是奇数个  则需要把最后一条个体直接加上 不交叉
        for i in range(0,lenparents-1,2):       #每次取出两条个体 如果是奇数个则长度要减去 一  range函数不会取最后一个
            single_now_popu = parents[i]   #取出当前选中的两个父代中的第一个
            single_next_popu = parents[i+1]#取出当前选中的两个父代中的第二个
            children.append([]) #子代添加两行  稍后存取新的种群
            children.append([])
            index_content = list(set(single_now_popu).intersection(set(single_next_popu))) #第一条路经与第二条路经重复的内容
            num_rep = len(index_content)          #重复内容的个数
            if random.random() < pc and num_rep>=3:
                content = index_content[random.randint(0,num_rep-1)]   #随机选取一个重复的内容
                now_index = single_now_popu.index(content)  #重复内容在第一个父代中的索引
                next_index = single_next_popu.index(content)#重复内容在第二个父代中的索引
                children[i] = single_now_popu[0:now_index + 1] + single_next_popu[next_index + 1:]
                children[i+1] = single_next_popu[0:next_index + 1] + single_now_popu[now_index + 1:]
            else:
                children[i] = parents[i]
                children[i+1] = parents[i+1]
        if parity == 1:     #如果是个奇数  为了保证种群规模不变 需要加上最后一条
            children.append([]) #子代在添加一行
            children[-1] = parents[-1] #直接遗传给下一代
        return children
    

    这里如果pc概率小于随机小数,不进行交叉操作的话,就把原来的路径信息直接赋给下一代,因为要保证交叉操作不会减小种群规模。还有一点要注意,如果父代种群有奇数条路径,可以看到,代码中是每次选择相邻的两条路径进行交叉,只能交叉偶数个父代,那么最后会剩余一个父代没法交叉,这里就是直接遗传给下一代了。

        if parity == 1:     #如果是个奇数  为了保证种群规模不变 需要加上最后一条
            children.append([]) #子代在添加一行
            children[-1] = parents[-1] #直接遗传给下一代
    

    小结

      交叉操作是肯定不会减小种群的规模的,并且还会使种群多样化,在写这篇博客时,忽然想到,如果交叉时选择的父代是随机选择,指定交叉次数,会不会能使种群的规模达到我所指定的规模?还有一个之前就想到的问题,如果相同的元素除了起始点和交叉点之外大于等于2个,那么父代交叉完一次后,子代可能还会有重复的元素,子代是否能再次交叉一次?直至两个子代没有相同元素。这些操作能够使种群更加多样化,并且有可能解决种群规模在选择操作时减小的问题。这里记下,之后检验。

    6、遗传学之变异

      遗传学基因突变是指基因中的某一段发生变化,在此路径规划中的变异就是指一条路径的某一段发生了变化,在代码中的思想为:给定一个变异概率pm,然后产生0-1之间的一个随机数,并和变异概率pm比较,若pm大于随机小数则进行变异操作。
      变异方法是随机选取路径中除起点和终点以外的两个栅格,去除这两个栅格之间的路径,然后以这两个栅格为相邻点,使用初始化种群中的路径连续化步骤将这两个点进行连续操作。此时有可能无法产生连续的路径,则需要重新选择两个栅格执行以上操作,直到完成变异操作。这种重新生成的连续路径就是变异之后的连续路径。
    变异代码如下:

    #step6 变异
    def mutation(children,pm):
        '''
        :param children: 子代种群
        :param pm: 变异概率
        :return: 返回变异后的新种群
        '''
    
        row = len(children)   #子代有多少行   即多少条路经
        new_popu = []
        for i in range(row):#提取出来每一行
    
            single_popu = children[i]
            if random.random()<pm:#小于变异概率   就变异
                col = len(single_popu)#每一行的长度 即列数  也就是这条路径有多少节点
                first = random.randint(1,col-2) #随机选取两个指针
                second = random.randint(1,col-2)
                if first != second :    #判断两个指针是否相同  不相同的话把两个指针中间的部分删除 在进行连续化
                    #判断一下指针大小  便于切片
                    if(first<second):
                        single_popu = single_popu[0:first]+single_popu[second+1:]
                    else :
                        single_popu = single_popu[0:second] + single_popu[first+1:]
                temp = population.Generate_Continuous_Path(single_popu)#连续化
                if temp!= []:
                    new_popu.append(temp)
            else:       #不变异的话就将原来的个体直接赋值
                new_popu.append(single_popu)
        return new_popu
    

    小结

      额。。。这个代码贴上来之后,我好像发现了种群规模减小的原因了。原因就是:此时有可能无法产生连续的路径,则需要重新选择两个栅格执行以上操作,直到完成变异操作。这句话我在程序中忽略了,对于无法生成路径的变异个体,我直接舍弃了,所以会导致变异丢失种群数量。在程序中应该判断如下:

    if temp == []:#路径为空,说明没有生成连续路径
    	'''
    	这段里面应该在进行循环变异知道变异成功
    	'''
    	pass
    

    7、更新种群以及输出结果

      至此,选择、交叉、变异全部完成,种群更新完毕。若没有出现最优路径或者超出迭代次数,则继续进行之前的操作。若是出现最优路径或迭代完毕则从种群中找到适应度值最小的个体进行输出,并将栅格地图中的编码用’*‘代替。画出清晰的路径图。
    这段代码在主函数中呈现:

    if __name__ == '__main__':
        print('原始地图')
        population = Population(10,10,20,100)   #实例化对象 并生成初始带随机障碍的地图
    
        popu = population.Population_Init()  #种群初始化  得到np条初始可行路径
        # print(popu,'连续路径')#打印出np条可行路径
    
        for i in range(200):    #迭代200代
            lastpopu = popu #上次的种群先保存起来
            value = calvalue(popu,population.col) #计算适应度值
            # print(value,'适应度值')#打印出np条路径的适应度值
    
            new = selection(popu,value)#选择 按适应度返回新的种群
            # print(new,'新种群')
            # value = calvalue(new,population.col) #计算适应度值
            # print(value,'新种群适应度值')#打印出np条路径的适应度值
            child = cross(new,0.8)  #交叉  产生子代
            # print(child,'交叉子代')
            # value = calvalue(child,population.col) #计算适应度值
            # print(value,'子代适应度值')#打印出np条路径的适应度值
            popu = mutation(child,0.8)#变异 产生子代 并更新原始种群
    
            if popu == []:  #如果本次的种群成空了  则把上次的种群信息拿出来,并迭代结束
                popu = lastpopu
                break
            # print('第',i,'次迭代后的种群为:',popu)
        if popu == []:  #迭代完成后先判断 若迭代完成了 但是种群路径还是空 说明可能是没有路径
            print('无路径')
        else:
            value = calvalue(popu,population.col) #计算适应度值
            minnum = value[0]
            for i in range(len(value)):#找出最小的适应度函数值
                if value[i] < minnum:#小于最小适应度  则替代
                    minnum = value[i]
            popu = popu[value.index(minnum)]#取出最小的适应度值的个体
    
            for i in popu:
                x = (i) // (population.map_start.col)  # 行 解码  算出此条路经的第i个元素的直角坐标
                y = (i) % (population.map_start.col)  # 列
                population.map_start.data[x][y] = '*'   #将路径用*表示
            print('\n规划地图')
            for i in range(population.map_start.row):   #显示路径
                for j in range(population.map_start.col):
                    print(population.map_start.data[i][j],end=' ')
                print('')
            print('最短路径值为:',minnum)
            print('最短路径为:',popu)
    

    结果如下图,上面的数组为生成的带障碍的原始地图,下面的为规划好路径的地图,最后输出最优路径的路径长度以及路径编码号。
    结果图

    四、代码工程文档

      理论再多,不如直接上可以用的代码,下面的代码为整体的工程代码,直接可Ctrl+C->Ctrl+V运行:

    # -*- coding: UTF-8 -*-
    '''*******************************
    @ 开发人员:Mr.Zs
    @ 开发时间:2020/5/18 19:22
    @ 开发环境:PyCharm
    @ 项目名称:算法类总工程->遗传算法路径规划V1.0.py
    ******************************'''
    import random   #生成随机整数
    import numpy as np #生成随机小数
    import math #用于计算除法 取整等运算
    print(r'''
    遗传算法是对参数集合的编码而非针对参数本身开始进化
    遗传算法是从问题解的编码组开始,而非单个解开始搜索
    
    step1:建立地图
    step2:初始化种群(随机生成若干条从起点能够到达终点的路径(可行解),每一条可行路径为一个个体)
    step3:计算个体的适应度值
    step4:选择适应度合适的个体进入下一代
    step5:交叉
    step6:变异
    step7:更新种群,若没有出现最优个体,则转至step3
    step8:输出最优的个体作为最优解
    参考文献:https://blog.csdn.net/qq_40870689/article/details/86916910?ops_request_misc=&request_id=&biz_id=102&utm_term=%E9%81%97%E4%BC%A0%E7%AE%97%E6%B3%95%E8%B7%AF%E5%BE%84%E8%A7%84%E5%88%92&utm_medium=distribute.pc_search_result.none-task-blog-2~all~sobaiduweb~default-1-86916910
    ''')
    
    #这个对象专门滤除路径中的回头路
    class Fiter:
        def __init__(self):
            self.b = 1  # 标志位
    
        def function(self, a):  # 定义一个函数
            for i in a:  # 遍历列表中的内容
                a = a[a.index(i) + 1:]  # 把当前内容索引的后面的内容剪切下来  因为前面的已经比对过了
                if i in a:  # 如果当前内容与后面有重复
                    return i, 1  # 返回当前重复的内容 以及标志位1
                else:  # 没有重复就不用管  继续for循环
                    pass
            return 0, 0  # 全部遍历完  没有重复的就返回0 这里返回两个0 是因为返回的数量要保持一致
    
        def fiter(self, a):
            while (self.b == 1):  # 标志位一直是 1  则说明有重复的内容
                (i, self.b) = self.function(a)  # 此时接受函数接收 返回值 i是重复的内容 b是标志位
                c = [j for j, x in enumerate(a) if x == i]  # 将重复内容的索引全部添加进c列表中
                a = a[0:c[0]] + a[c[-1]:]  # a列表切片在重组
            return a
    fiter = Fiter()#实例化对象
    #将地图上的点抽象话,可以用x表示横坐标 y表示纵坐标
    class Point:
        def __init__(self, x, y):
            self.x = x
            self.y = y
    
        def __eq__(self, other): #函数重载  #判断两个坐标  的  值 是否一样
            if((self.x == other.x )and (self.y == other.y)):
                return  True
            else:
                return False
        def __ne__(self, other):
            pass
    
    #step1
    class Map():
        '''
        :param:地图类
        :param:使用时需要传入行列两个参数 再实例化
        '''
    
        def __init__(self,row,col):
            '''
            :param:row::行
            :param:col::列
            '''
            self.start_point = Point(0,0)
            self.end_point = Point(row-1,col-1)
            self.data = []
            self.row = row
            self.col = col
        def map_init(self):
            '''
            :param:创建栅格地图row行*col列 的矩阵
            '''
            self.data = [[0 for i in range(self.col)] for j in range(self.row)]
            # for i in range(self.row):
            #     for j in range(self.col):
            #         print(self.data[i][j],end=' ')
            #     print('')
        def map_Obstacle(self,num):
            '''
            :param:num:地图障碍物数量
            :return:返回包含障碍物的地图数据
            '''
            self.num = num
            for i in range(self.num):#生成小于等于num个障碍
                self.data[random.randint(0,self.row-1)][random.randint(0,self.col-1)] = 1
            if self.data[0][0] == 1:        #判断顶点位置是否是障碍  若是 修改成可通行
                self.data[0][0] = 0
    
            if self.data[self.row-1][0] == 1:
                self.data[self.row-1][0] = 0
    
            if self.data[0][self.col-1] == 1:
                self.data[0][self.col - 1] = 0
    
            if self.data[self.row-1][self.col - 1] == 1:
                self.data[self.row - 1][self.col - 1] = 0
            for i in range(self.row):       #显示出来
                for j in range(self.col):
                    print(self.data[i][j],end=' ')
                print('')
            return self.data
    #step2 初始化种群 也是最难的一步
    class Population():
        def __init__(self,row,col,num,NP):
            self.row = row
            self.col = col
            self.num = num
            self.NP = NP
            self.p_start = 0  # 起始点序号  10进制编码
            self.p_end = self.row * self.col - 1  # 终点序号 10进制编码
            self.xs = (self.p_start) // (self.col)  # 行
            self.ys = (self.p_start) % (self.col)  # 列
            self.xe = (self.p_end) // (self.col)  # 终点所在的行
            self.ye = (self.p_end) % (self.col)
            self.map_start = Map(self.row, self.col)  # Map的实例化 主函数中可以不再初始化地图
            self.map_start.map_init()  # map初始化 生成0矩阵
            self.map = self.map_start.map_Obstacle(self.num)  # 生成带障碍的map
    
            self.can = []  # 这个列表存放整个地图中不是障碍物的点的坐标 按行搜索
            self.popu = [[0 for i in range(self.col)]
                         for j in range(self.NP)]  # 种群的列表,包含NP个间断可行解,即从起点到终点的路径
            # self.popu = []#存放可行解  可行路径
            self.end_popu = []
    
        def Population_Init(self):
            '''
            :return:无返回值,作用是找出NP条不连续可行解
            '''
    
    
            for i in range(self.NP):       #添加NP条路径
                j = 0
                for xk in range(0,self.row):  #多少行
                    self.can = []             #清空can列表   用来缓存当前行的可行点
                    for yk in range(0,self.col):   #从起点0 开始的列开始搜寻到终点的列  也就是此行中的自由点
                        num = (yk) + (xk) * self.col
                        if self.map_start.data[xk][yk] == 0:
                            self.can.append(num)
                    # print(self.can,'自由点\n')
                    length = len(self.can)     #此行的长度  随机生成指针取出坐标
                    # print(length)
                    self.popu[i][j] = (self.can[random.randint(0,length-1)])
                    j += 1#j从1 开始  是因为 0 的元素时起点编号
                self.popu[i][0] = self.p_start  # 每一行的第一个元素设置为起点序号
                self.popu[i][-1] = self.p_end
    
                temp = self.Generate_Continuous_Path(self.popu[i])#将返回的一条路经 添加进end中
                if temp != []:      #将返回的空列表路径滤除
                    temp = fiter.fiter(temp)    #滤除重复节点路径
                    self.end_popu.append(temp)
                # print(self.end_popu, end='\n')
            # print('测试1',self.popu,end='\n')
            return self.end_popu
        # @staticmethod
        def Generate_Continuous_Path(self,old_popu):#生成连续的路径个体
            '''
            :param old_popu: 未进行连续化的一条路径
            :return:        无返回                   # 已经连续化的一条路径
            '''
            self.new_popu = old_popu    #传进来的参数就是一行的数组  一条路径
            self.flag = 0
            self.lengh = len(self.new_popu)  #先将第一条路经的长度取出来
            i = 0
            # print("lengh =",self.lengh )
            while i!= self.lengh-1:       #i 不等于 当前行的长度减去1  从0计数  这里有问题 待修改
                x_now = (self.new_popu[i]) // (self.col)  # 行 解码  算出此条路经的第i个元素的直角坐标
                y_now = (self.new_popu[i]) % (self.col)  # 列
                x_next =  (self.new_popu[i+1]) // (self.col) #计算此条路经中下一个点的坐标
                y_next =  (self.new_popu[i+1]) % (self.col)
                #最大迭代次数
                max_iteration = 0
    
                #判断下一个点与当前点的坐标是否连续 等于1 为连续
                while max(abs(x_next - x_now), abs(y_next - y_now)) != 1:
                    x_insert = math.ceil((x_next + x_now) // 2)      #行
                    y_insert = math.ceil((y_next + y_now) // 2) #ceil向上取整数     #列
                    # print("x_insert = ",x_insert,"\ty_insert = ",y_insert)
                    flag1 = 0
    
                    if self.map_start.data[x_insert][y_insert] == 0:  #插入的这个坐标为0 可走
                        num_insert = (y_insert) + (x_insert) * self.col #计算出插入坐标的编码
                        self.new_popu.insert(i+1,num_insert)
                        # print(self.new_popu)
                        # print(num_insert)
                    else:#插入的栅格为障碍  判断插入的栅格上下左右是否为障碍,以及是否在路径中,若不是障碍且不在路径中,就插入
                        #判断下方
                        if (x_insert + 1 < self.row)and flag1 == 0:       #保证坐标是在地图上的
                            if ((self.map_start.data[x_insert+1][y_insert] == 0)#下方不是障碍物
                                and (((y_insert) + (x_insert+1) * self.col) not in self.new_popu)):#编号不在已知路径中
                                num_insert = (y_insert) + (x_insert+1) * self.col  #计算下方的编号
                                self.new_popu.insert(i + 1, num_insert) #插入编号
                                flag1 = 1       #设置标志位 避免下面重复插入
    
                                # print('下方插入',num_insert)
                        #判断右方
                        if (y_insert + 1 < self.col)and flag1 == 0:  # 保证坐标是在地图上的 并且前面没有插入
                            if ((self.map_start.data[x_insert][y_insert+1] == 0)#右方不是障碍物
                                and (((y_insert+1) + (x_insert) * self.col) not in self.new_popu)):#编号不在已知路径中
                                num_insert = (y_insert+1) + (x_insert) * self.col  #计算右方的编号
                                self.new_popu.insert(i + 1, num_insert) #插入编号
                                flag1 = 1  # 设置标志位 避免下面重复插入
                                # print('右方插入',num_insert)
                        #判断上方
                        if (x_insert - 1 > 0) and flag1 == 0:  # 保证坐标是在地图上的
                            if ((self.map_start.data[x_insert-1][y_insert] == 0)#右方不是障碍物
                                and (((y_insert) + (x_insert-1) * self.col) not in self.new_popu)):#编号不在已知路径中
                                num_insert = (y_insert) + (x_insert-1) * self.col  #计算右方的编号
                                self.new_popu.insert(i + 1, num_insert) #插入编号
                                flag1 = 1  # 设置标志位 避免下面重复插入
                                # print('上方插入',num_insert)
                        #判断左方
                        if (y_insert - 1 > 0)and flag1 == 0:  # 保证坐标是在地图上的
                            if ((self.map_start.data[x_insert][y_insert-1] == 0)#右方不是障碍物
                                and (((y_insert-1) + (x_insert) * self.col) not in self.new_popu)):#编号不在已知路径中
                                num_insert = (y_insert-1) + (x_insert) * self.col  #计算右方的编号
                                self.new_popu.insert(i + 1, num_insert) #插入编号
                                flag1 = 1  # 设置标志位 避免下面重复插入
                                # print('左方插入',num_insert)
                        if flag1 == 0:  #如果前面没有插入新点  说明这条路径不对 删除
                            self.new_popu = []
                            break
                    x_next = num_insert//self.col
                    y_next = num_insert%self.col
                    # x_next = x_insert
                    # y_next = y_insert
                    max_iteration += 1#迭代次数+1
                    if max_iteration > 20:
                        self.new_popu = []  #超出迭代次数 说明此条路经可能无法进行连续   删除路径
                        break
                if self.new_popu == []:
                    break
                self.lengh = len(self.new_popu)
                i = i+1
            # print(self.new_popu,'连续')
            return  self.new_popu#返回的是一条路径
    #step3 计算适应度函数
    def calvalue(popu,col):
        '''
        :param popu: 传入种群信息
        :param col: 这个参数是地图的列数 ,因为要解码用到
        :return:    返回的是每一条路径长度组成的列表
        '''
        hang = len(popu)#计算行
        value = [] #存放计算出来的路径长度值
        for i in range(hang):#从第0行开始 到最后一行
            value.append(0)
            single_popu = popu[i] #把当前行的元素赋给 single_popu 之后操作
            single_lengh = len(single_popu)#将当前行的长度拿出来
            for j in range(single_lengh-1):#从第一个元素计算到倒数第一个元素
                x_now = (single_popu[j]) // (col)  # 行 解码  算出此条路经的第i个元素的直角坐标
                y_now = (single_popu[j]) % (col)  # 列
                x_next = (single_popu[j + 1]) // (col)  # 计算此条路经中下一个点的坐标
                y_next = (single_popu[j + 1]) % (col)
                if abs(x_now - x_next) + abs(y_now - y_next) == 1:#路径上下左右连续 不是对角的 则路径长度为1
                    value[i] = value[i]+1
                elif max(abs(x_now - x_next),abs(y_now - y_next))>=2:#惩罚函数 若跳跃或者穿过障碍
                    value[i] = value[i] + 100
                else:
                    value[i] = value[i]+1.4  #对角长度为根号2  即1.4
        return value
    #step4 选择
    def selection(pop,value):
        '''
        :param pop:种群
        :param value:适应度值列表
        :return:返回新的种群
        '''
    
        ###原来的方法会丢失种群数量
        now_value=[]#做倒数后的适应度值列表
        P_value = []  #存放适应度值占总比概率的列表
        random_deci = []#存放随机的小数
        new_popu = []       #选择后的种群
        sum_value = 0   #存放适应度的总和  计算轮盘赌的概率
        lengh = len(pop)#首先计算种群有多少行,也就是有多少初始路径
        for i in range(lengh):
            new_popu.append([]) #原始种群有多少个个体  新的种群就有多少
        for i in value:     #由于适应度值是距离  将需要的就是距离最短的  因此求倒数来排序
            now_value.append(1/i)  #将适应度取倒数   倒数越小适应度越小,路径越长  就可以按概率抛弃
            sum_value += (1/i)
        for i in now_value:#将每一个适应度值取出来  做比例  得出每个适应度值的占比概率 存到列表中
            P_value.append(i/sum_value)
        P_value = np.cumsum(P_value)#累加和 并排序 从小到大
        for i in range(lengh):
            random_deci.append(random.random())#产生 i个随即小数 存到列表中
        random_deci = sorted(random_deci)#从小到大排序
        fitin = 0
        newin = 0
        while(newin<lengh): #遍历每一行
            if random_deci[newin] < P_value[fitin]:
                new_popu[newin] = pop[fitin]
                newin += 1
            else:
                fitin += 1
        return new_popu
        # ####原来的方法会丢失种群数量
        # retain_rate = 0.6#适应度最优的保留几率
        # random_rate = 0.3#适应度差的保留几率
        # dict_my = dict(zip(value, pop))     #将路径距离和路径打包字典
        #
        # new_popu = []
        # sort_dis = sorted(dict_my)   #将字典按照键排序  也就是按照距离短-长排序
        # retain_lengh = int(len(pop)*retain_rate)      #适应度强的按照存活概率 保留下
        # temp = sort_dis[:retain_lengh]      #缓存下保留的距离,待会将距离对应的字典值取出来
        # # print(temp,'优秀保留')
        # for i in sort_dis[retain_lengh:]:   #距离长的按照随机概率保留
        #     if random.random() < random_rate:
        #         temp.append(i)
        # #temp现在存放的就是进入下一代的种群信息,他是距离信息,下一步通过字典将种群提取出来
        # for i in temp:
        #     new_popu.append(dict_my[i]) #至此种群选取完
        # return new_popu
    #step5 交叉   拟采用单点交叉
    def cross(parents,pc):
        '''
        :param parents: 交叉的父类
        :param pc:   交叉概率
        :return:
        '''
        children = []  #首先创建一个子代 空列表 存放交叉后的种群
        single_popu_index_list = []#存放重复内容的指针
        lenparents = len(parents)  #先提取出父代的个数  因为要配对 奇数个的话会剩余一个
        parity = lenparents % 2 #计算出长度奇偶性  parity= 1 说明是奇数个  则需要把最后一条个体直接加上 不交叉
        for i in range(0,lenparents-1,2):       #每次取出两条个体 如果是奇数个则长度要减去 一  range函数不会取最后一个
            single_now_popu = parents[i]   #取出当前选中的两个父代中的第一个
            single_next_popu = parents[i+1]#取出当前选中的两个父代中的第二个
            children.append([]) #子代添加两行  稍后存取新的种群
            children.append([])
            index_content = list(set(single_now_popu).intersection(set(single_next_popu))) #第一条路经与第二条路经重复的内容
            num_rep = len(index_content)          #重复内容的个数
            if random.random() < pc and num_rep>=3:
                content = index_content[random.randint(0,num_rep-1)]   #随机选取一个重复的内容
                now_index = single_now_popu.index(content)  #重复内容在第一个父代中的索引
                next_index = single_next_popu.index(content)#重复内容在第二个父代中的索引
                children[i] = single_now_popu[0:now_index + 1] + single_next_popu[next_index + 1:]
                children[i+1] = single_next_popu[0:next_index + 1] + single_now_popu[now_index + 1:]
            else:
                children[i] = parents[i]
                children[i+1] = parents[i+1]
        if parity == 1:     #如果是个奇数  为了保证种群规模不变 需要加上最后一条
            children.append([]) #子代在添加一行
            children[-1] = parents[-1]
        return children
    #step6 变异
    def mutation(children,pm):
        '''
        :param children: 子代种群
        :param pm: 变异概率
        :return: 返回变异后的新种群
        '''
    
        row = len(children)   #子代有多少行   即多少条路经
        new_popu = []
        for i in range(row):#提取出来每一行
    
            single_popu = children[i]
            if random.random()<pm:#小于变异概率   就变异
                col = len(single_popu)#每一行的长度 即列数  也就是这条路径有多少节点
                first = random.randint(1,col-2) #随机选取两个指针
                second = random.randint(1,col-2)
                if first != second :    #判断两个指针是否相同  不相同的话把两个指针中间的部分删除 在进行连续化
                    #判断一下指针大小  便于切片
                    if(first<second):
                        single_popu = single_popu[0:first]+single_popu[second+1:]
                    else :
                        single_popu = single_popu[0:second] + single_popu[first+1:]
                temp = population.Generate_Continuous_Path(single_popu)#连续化
                if temp!= []:
                    new_popu.append(temp)
            else:       #不变异的话就将原来的个体直接赋值
                new_popu.append(single_popu)
        return new_popu
    if __name__ == '__main__':
        print('原始地图')
        population = Population(10,10,20,100)   #实例化对象 并生成初始带随机障碍的地图
    
        popu = population.Population_Init()  #种群初始化  得到np条初始可行路径
        # print(popu,'连续路径')#打印出np条可行路径
    
        for i in range(200):    #迭代200代
            lastpopu = popu #上次的种群先保存起来
            value = calvalue(popu,population.col) #计算适应度值
            # print(value,'适应度值')#打印出np条路径的适应度值
    
            new = selection(popu,value)#选择 按适应度返回新的种群
            # print(new,'新种群')
            # value = calvalue(new,population.col) #计算适应度值
            # print(value,'新种群适应度值')#打印出np条路径的适应度值
            child = cross(new,0.8)  #交叉  产生子代
            # print(child,'交叉子代')
            # value = calvalue(child,population.col) #计算适应度值
            # print(value,'子代适应度值')#打印出np条路径的适应度值
            popu = mutation(child,0.8)#变异 产生子代 并更新原始种群
    
            if popu == []:  #如果本次的种群成空了  则把上次的种群信息拿出来,并迭代结束
                popu = lastpopu
                break
            # print('第',i,'次迭代后的种群为:',popu)
        if popu == []:  #迭代完成后先判断 若迭代完成了 但是种群路径还是空 说明可能是没有路径
            print('无路径')
        else:
            value = calvalue(popu,population.col) #计算适应度值
            minnum = value[0]
            for i in range(len(value)):#找出最小的适应度函数值
                if value[i] < minnum:#小于最小适应度  则替代
                    minnum = value[i]
            popu = popu[value.index(minnum)]#取出最小的适应度值的个体
    
            for i in popu:
                x = (i) // (population.map_start.col)  # 行 解码  算出此条路经的第i个元素的直角坐标
                y = (i) % (population.map_start.col)  # 列
                population.map_start.data[x][y] = '*'   #将路径用*表示
            print('\n规划地图')
            for i in range(population.map_start.row):   #显示路径
                for j in range(population.map_start.col):
                    print(population.map_start.data[i][j],end=' ')
                print('')
            print('最短路径值为:',minnum)
            print('最短路径为:',popu)
    

    结束语

      代码还有个待改进的地方,有时候明明有路径,但是却会输出“无路径”,可能是与路径连续化这个代码有关系

    问题解决

    1、解决种群规模随迭代次数增加而减小的问题

      解决方法与上文猜测一致,将变异操作的代码改为如下所示即可:

    #step6 变异
    def mutation(children,pm):
        '''
        :param children: 子代种群
        :param pm: 变异概率
        :return: 返回变异后的新种群
        '''
    
        row = len(children)   #子代有多少行   即多少条路经
        new_popu = []
        for i in range(row):#提取出来每一行
            temp = [] #定义一个缓存路径的空列表
            while(temp == []): #在这里加入死循环,一直等待变异成功
                single_popu = children[i]
                if random.random()<pm:#小于变异概率   就变异
                    col = len(single_popu)#每一行的长度 即列数  也就是这条路径有多少节点
                    first = random.randint(1,col-2) #随机选取两个指针
                    second = random.randint(1,col-2)
                    if first != second :    #判断两个指针是否相同  不相同的话把两个指针中间的部分删除 在进行连续化
                        #判断一下指针大小  便于切片
                        if(first<second):
                            single_popu = single_popu[0:first]+single_popu[second+1:]
                        else :
                            single_popu = single_popu[0:second] + single_popu[first+1:]
                    temp = population.Generate_Continuous_Path(single_popu)#连续化
                    if temp!= []:
                        new_popu.append(temp)
    
                else:       #不变异的话就将原来的个体直接赋值
                    new_popu.append(single_popu)
        return new_popu
    
    展开全文
  • 使用Python的scikit-opt包ÿ...设计遗传算法,求解以下运输规划问题,数据见excel文档。 <p><img alt="" height="588" src="https://img-ask.csdnimg.cn/upload/1618891187386.png" width="902" /></p>
  • 遗传算法路径规划

    2018-01-19 10:37:12
    用 MATLAB 实现基于遗传算法路径规划,只是一个9x9方格的小程序。
  • 自从上一篇博客详细讲解了Python遗传和进化算法工具箱及其在带约束的单目标函数值优化中的应用之后,我经过不断学习工具箱的官方文档以及对源码的研究,逐步更加地掌握如何利用遗传算法求解更多有趣的问题了。...

    前言

    自从上一篇博客详细讲解了Python遗传和进化算法工具箱及其在带约束的单目标函数值优化中的应用之后,我经过不断学习工具箱的官方文档以及对源码的研究,逐步更加地掌握如何利用遗传算法求解更多有趣的问题了。

    首先简单回顾一下Python高性能实用型遗传和进化算法工具箱的用法。对于一个优化问题,需要做两个步骤便可进行求解:Step1:自定义问题类;Step2:编写执行脚本调用Geatpy进化算法模板对问题进行求解。在上一篇博客曾“详细”介绍过具体的用法:https://blog.csdn.net/weixin_37790882/article/details/84034956,但完整的中文教程可以参考官方文档

    下面切入主题:

    本文的“最短路径问题”专指图的最短路径问题,而且只研究单目标的最短路径问题。(实际上最短路径问题还可以是多目标的),参考用书推荐:《网络模型与多目标遗传算法》。这本书我已经上传到此处,可以直接下载到:

    https://download.csdn.net/download/weixin_37790882/11700202

    图的最短路径问题从广义上说其实有很多种类型,比如:

    • 确定起点和终点的有向图和无向图的最短路径问题。
    • 旅行商最短路径问题。
    • 车辆配送中的最短路径问题。
    • 物流中心选址的最短路径问题。

    其中确定起点和终点的无向图的最短路径问题可以延伸为机器人路径规划问题。也就是说,机器人避开障碍物从起点走向终点的那种路径规划问题本质上是一种复杂一些的“确定起点和终点的无向图的最短路径问题”。

    本文主要讲解如何利用遗传算法求解“确定起点和终点的有向图最短路径问题”。对于无向图的,另外再写博客以机器人路径规划作为案例来展开叙述。

    正文

    问题举例:

    以《网络模型与多目标遗传算法》一书中的一个小案例为例,如图所示,当从结点1驶向结点10时,我们经常会考虑怎样选择路径使得花最短的距离达到目的地。此时不需要像旅行商问题(TSP)那样,此类问题不需要要求所有的地点都访问一遍,而只需要想办法用最短的距离从起点走到终点即可。

    书中讲述了两种遗传算法编码方式对上述问题进行求解。

    法一

    利用二进制编码。假设 Xij 是有向图中所有边的一个集合,那么其中一种可行的路径对应的染色体可以是:

    其中元素为1表示对应的边被“激活”,即最终路径上包含这条边。于是上面的染色体所代表的路径为:1 → 3 → 6 → 7 → 8 → 10。

    这种编码方式最为直观,但有个很大的弊端:在进化过程中染色体往往无法表达一个合法的路径。于是这种情况下需要给遗传算法添加很多约束,此时便增加了遗传算法寻优的难度。

    法二

    利用“基于优先级编码”。这种编码方式是Gen & Lin 在2005年提出的有利于很好地求解图的路径规划问题的编码方法。参考文献:(Lin L, Gen M. Priority-Based Genetic Algorithm for Shortest Path Routing Problem in OSPF[M]// Intelligent and Evolutionary Systems. 2009.)这篇文献可能比较难下载到,我将其上传到这里了,可以直接下载查看:

    https://download.csdn.net/download/weixin_37790882/11699834

    下面来详细讲解这种编码方式

    “基于优先级编码”是一种间接编码方式,这意味着染色体并不能直接表示路径,此时需要利用额外的数据来进行解码,解码后才表示一个路径。这种编码方式有个特点是染色体的每一位上的元素是互不相等的,这意味着这种编码方式具备“排列编码”的特征。(排列编码即比如从1到10的10个数中随机挑选若干个数组成的一个排列。)

    以一条染色体为例,看看“基于优先级编码”的染色体是如何表示有向图中的一条路径的:

    上面这条染色体是遗传算法中随机生成的(并非最优),结点优先级并不是指结点的访问先后顺序,而是结点的优先级,是给后面解码用的。为什么染色体长度是10?因为此时染色体每一位存储的是图的每个结点的优先级,因此染色体的长度需要和图的结点数一致。当然地,由于上面题目规定了是从结点1开始走,故设染色体长度为9也行。这里为了能和下标更直观地对应,就设染色体长度与图的结点数目一致。

    解码需要一个集合nodes,用于存储以各结点为起点的有向边的终点,即各个结点下一步可达的结点。本题的集合nodes为:nodes = [[], [2,3], [3,4,5], [5,6], [7,8], [4,6], [7,9], [8,9], [9,10], [10]]。因为python中列表的下标是从0数起的,而本题的结点是从1数起的,为了能直接对应,故上面的nodes的第0号元素设置为[],表示无意义。解析一下nodes的组成:第1号元素是[2,3],表示题目的有向图中1号结点下一步可达的结点是2和3。nodes的第2号元素是[3,4,5],表示2号结点下一步可达的结点是3,4,5。以此类推。

    于是上面的染色体[7, 3, 4, 6, 2, 5, 8, 10, 1, 9]的详细解码过程如下

    从1号结点出发,在nodes中下标为1的元素是[2,3],表示1号结点下一步可以去2号结点或3号结点。此时从染色体中找到这两个结点对应的优先级分别为3和4,如图所示:

    从中选出具有更高优先级的结点3作为结点1下一步需要访问的结点。

    紧接着,在nodes中下标为3的元素是[5,6],表示3号结点下一步可以去5号或6号结点。此时从染色体中找到这两个结点对应的优先级分别为2和5,如图所示:

    从中选出具有更高优先级的结点6作为结点1下一步需要访问的结点。

    以此类推,最终得到完整的访问路径为:

    1 → 3 → 6 → 7 → 8 → 10。

    有了解码得到路径,如何求个体的适应度?

    想要求个体的适应度,首先要定一个优化目标。本题是要让路径最短,于是我们便要根据访问结点求出路径长度,把这个作为优化目标。有了优化目标,便可利用基于目标函数值排序的适应度分配(ranking)求出适应度值。当然,这类单目标优化问题也可以直接让个体的适应度等于优化目标函数值(即路径长度)。而路径长度即为访问路径上各条有向边的权值之和。

    代码实现

    首先编写问题类,把待优化模型写在问题类中。然后编写执行脚本,调用“soea_SEGA_templet“(增强精英保留的遗传算法模板)进行进化优化。该算法模板的源码及算法流程详见

    https://github.com/geatpy-dev/geatpy/blob/master/geatpy/templates/soeas/GA/SEGA/soea_SEGA_templet.py

    由于本题比较简单,故用4个个体、10代的进化即可得到很好的结果。完整的实验代码如下:

    # -*- coding: utf-8 -*-
    import numpy as np
    import geatpy as ea
    
    class MyProblem(ea.Problem): # 继承Problem父类
        def __init__(self):
            name = 'Shortest_Path' # 初始化name(函数名称,可以随意设置)
            M = 1 # 初始化M(目标维数)
            maxormins = [1] # 初始化maxormins(目标最小最大化标记列表,1:最小化该目标;-1:最大化该目标)
            Dim = 10 # 初始化Dim(决策变量维数)
            varTypes = [1] * Dim # 初始化varTypes(决策变量的类型,元素为0表示对应的变量是连续的;1表示是离散的)
            lb = [0] * Dim # 决策变量下界
            ub = [9] * Dim # 决策变量上界
            lbin = [1] * Dim # 决策变量下边界
            ubin = [1] * Dim # 决策变量上边界
            # 调用父类构造方法完成实例化
            ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)
            # 设置每一个结点下一步可达的结点(结点从1开始数,因此列表nodes的第0号元素设为空列表表示无意义)
            self.nodes = [[], [2,3], [3,4,5], [5,6], [7,8], [4,6], [7,9], [8,9], [9,10], [10]]
            # 设置有向图中各条边的权重
            self.weights = {'(1, 2)':36, '(1, 3)':27, '(2, 4)':18, '(2, 5)':20, '(2, 3)':13, '(3, 5)':12, '(3, 6)':23,
                             '(4, 7)':11, '(4, 8)':32, '(5, 4)':16, '(5, 6)':30, '(6, 7)':12, '(6, 9)':38, '(7, 8)':20,
                             '(7, 9)':32, '(8, 9)':15, '(8, 10)':24, '(9, 10)':13}
        
        def decode(self, priority): # 将优先级编码的染色体解码得到一条从节点1到节点10的可行路径
            edges = [] # 存储边
            path = [1] # 结点1是路径起点
            while not path[-1] == 10: # 开始从起点走到终点
                currentNode = path[-1] # 得到当前所在的结点编号
                nextNodes = self.nodes[currentNode] # 获取下一步可达的结点组成的列表
                chooseNode = nextNodes[np.argmax(priority[np.array(nextNodes) - 1])] # 从NextNodes中选择优先级更高的结点作为下一步要访问的结点,因为结点从1数起,而下标从0数起,因此要减去1
                path.append(chooseNode)
                edges.append((currentNode, chooseNode))
            return path, edges
    
        def aimFunc(self, pop): # 目标函数
            pop.ObjV = np.zeros((pop.sizes, 1)) # 初始化ObjV
            for i in range(pop.sizes): # 遍历种群的每个个体,分别计算各个个体的目标函数值
                priority = pop.Phen[i, :]
                path, edges = self.decode(priority) # 将优先级编码的染色体解码得到访问路径及经过的边
                pathLen = 0
                for edge in edges:
                    key = str(edge) # 根据路径得到键值,以便根据键值找到路径对应的长度
                    if not key in self.weights:
                        raise RuntimeError("Error in aimFunc: The path is invalid. (当前路径是无效的。)", path)
                    pathLen += self.weights[key] # 将该段路径长度加入
                pop.ObjV[i] = pathLen # 计算目标函数值,赋值给pop种群对象的ObjV属性
    
    ## 执行脚本
    if __name__ == "__main__":
        problem = MyProblem()                      # 生成问题对象
        """=================================种群设置================================="""
        Encoding = 'P'                             # 编码方式
        NIND = 4                                   # 种群规模
        Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges, problem.borders) # 创建区域描述器
        population = ea.Population(Encoding, Field, NIND) # 实例化种群对象(此时种群还没被初始化,仅仅是完成种群对象的实例化)
        """===============================算法参数设置==============================="""
        myAlgorithm = ea.soea_SEGA_templet(problem, population) # 实例化一个算法模板对象
        myAlgorithm.MAXGEN = 10                    # 最大进化代数
        myAlgorithm.recOper = ea.Xovox(XOVR = 0.8) # 设置交叉算子
        myAlgorithm.mutOper = ea.Mutinv(Pm = 0.1)  # 设置变异算子
        myAlgorithm.logTras = 1  # 设置每隔多少代记录日志,若设置成0则表示不记录日志
        myAlgorithm.verbose = True  # 设置是否打印输出日志信息
        myAlgorithm.drawing = 1  # 设置绘图方式(0:不绘图;1:绘制结果图;2:绘制目标空间过程动画;3:绘制决策空间过程动画)
        """===========================调用算法模板进行种群进化======================="""
        [BestIndi, population] = myAlgorithm.run()  # 执行算法模板,得到最优个体以及最后一代种群
        BestIndi.save()  # 把最优个体的信息保存到文件中
        """==================================输出结果============================="""
        print('用时:%f 秒' % myAlgorithm.passTime)
        print('评价次数:%d 次' % myAlgorithm.evalsNum)
        print('最短路程为:%s'%(BestIndi.ObjV[0][0]))
        print('最佳路线为:')
        best_journey, edges = problem.decode(BestIndi.Phen[0])
        for i in range(len(best_journey)):
            print(int(best_journey[i]), end = ' ')
        print()
    

    运行结果如下:

    在有向图中表现为:

    后记

    值得注意的是:上面题目中的有向图并不存在回路,实际上,复杂的有向图往往会存在许多回路,此时需要进行一定的处理来避免陷入回路当中,即避免一直在回路上“打转”。处理方式有很多,例如在解码过程中对已经访问过的结点的有限度进行惩罚等等。这里暂时就不深入探讨了,待之后讲述无向图最短路径及机器人寻路问题时再展开叙述。

    最后回顾一下上一篇博客提到的”遗传算法套路“:

    该套路实现了具体问题、使用的算法以及所调用的相关算子之间的脱耦。而Geatpy工具箱已经内置了众多进化算法模板类以及相关的算子,直接调用即可。对于实际问题的求解,只需关心如何把问题写在自定义问题类中就好了。

    更多详细的教程可以详见:http://geatpy.com/index.php/geatpy%E6%95%99%E7%A8%8B/

    后续我将继续学习和挖掘该工具箱的更多深入的用法。希望这篇文章在帮助自己记录学习点滴之余,也能帮助大家!

    展开全文
  • 一、TSP问题 TSP问题(Travelling Salesman Problem)即旅行商问题,又译...下面使用遗传算法对此TSP问题求近似解,不了解此算法的同学请先移步此知乎链接: 如何通俗易懂地解释遗传算法?有什么例子? 最终效果如下图

    一、TSP问题

    TSP问题(Travelling Salesman Problem)即旅行商问题,又译为旅行推销员问题、货郎担问题,是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市,他必须选择所要走的路径,路径的限制是每个城市只能拜访一次,而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。简单示例如下图所示:
    在这里插入图片描述

    二、求解算法

    下面使用遗传算法对此TSP问题求近似解,不了解此算法的同学请先移步此知乎链接:
    如何通俗易懂地解释遗传算法?有什么例子?

    最终效果如下图

    在这里插入图片描述
    主要运行代码如下:

    # coding:utf:8
    import numpy as np
    import random
    import matplotlib.pyplot as plt
    from matplotlib.ticker import FormatStrFormatter
    import operator
    import time
    import pandas as pd
    
    from tsp1.mutations import select_best_mutaion
    
    
    def main():
        global p_mutation, max_generation
        generation = 1
    
        population_cur = init_population()
        fitness = get_fitness(population_cur)
    
        time_start = time.time()
    
        # 终止条件
        while generation < max_generation:
    
            # 父代最好的留1/4活下来
            population_next = select_sorted_population(fitness, population_cur, population_size // 4)
    
            # 杂交
            for i in range(population_size):
                p1, p2 = selection(fitness, 2)
                child1, child2 = crossover(population_cur[p1], population_cur[p2])
    
                # 变异
                if random.random() < p_mutation:
                    child1 = select_best_mutaion(child1, distmat)
                if random.random() < p_mutation:
                    child2 = select_best_mutaion(child2, distmat)
    
                population_next.append(child1)
                population_next.append(child2)
    
            # 选出下一代的种群
            population_next = select_sorted_population(get_fitness(population_next), population_next, population_size)
    
            # 找出精英记录下来
            pre_max_fitness, pre_max_individual = get_elite(fitness, population_cur)
            record(pre_max_fitness)
    
            # 换代
            population_cur = population_next
            generation += 1
            # 更新fitness
            fitness = get_fitness(population_cur)
    
        # 记录并画图
        final_fitness, final_individual = get_elite(fitness, population_cur)
        record(final_fitness)
    
        time_end = time.time()
        print('进化花费时间:', time_end - time_start)
        print('最后的路径距离(km):', get_distance(final_individual) * 111)
    
        plot(final_individual)
    
        return
    
    
    # 排序,并且返回length长的population
    def select_sorted_population(fitness, population, length):
        global population_size
        sort_dict = {}
        for i in range(len(population)):
            sort_dict[(fitness[i], 1 / fitness[i])] = i
    
        sorted_key = sorted(sort_dict.keys(), key=operator.itemgetter(0), reverse=True)
    
        sorted_index = [sort_dict[i] for i in sorted_key]
        sorted_population = [population[i] for i in sorted_index]
    
        return sorted_population[:length]
    
    
    # 画图
    def plot(sequence):
        global record_distance, coordinates
    
        plt.figure(figsize=(15, 6))
        plt.subplot(121)
    
        plt.plot(record_distance)
        plt.ylabel('distance')
        plt.xlabel('iteration ')
    
        plt.subplot(122)
    
        x_list = []
        y_list = []
        for i in range(len(sequence)):
            lat = coordinates[sequence[i]][0]
            lon = coordinates[sequence[i]][1]
            x_list.append(lat)
            y_list.append(lon)
            for one in name.itertuples(index=False):
                lat1 = one[0]
                if lat1 == lat:
                    lon1 = one[1]
                    if lon1 == lon:
                        # 显示城市名
                        city = one[2]
                        plt.text(lat, lon, city)
                        print(city)
                        print(' |')
                        break
        lat = coordinates[sequence[0]][0]
        lon = coordinates[sequence[0]][1]
        x_list.append(lat)
        y_list.append(lon)
        for one in name.itertuples(index=False):
            lat1 = one[0]
            if lat1 == lat:
                lon1 = one[1]
                if lon1 == lon:
                    # 显示城市名
                    city = one[2]
                    print(city)
                    break
    
        plt.plot(x_list, y_list, 'c-', label='Route')
        plt.plot(x_list, y_list, 'ro', label='Location')
    
        # 防止科学计数法
        ax = plt.gca()
        ax.xaxis.set_major_formatter(FormatStrFormatter('%.2f'))
        ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
    
        plt.xlabel("Longitude")
        plt.ylabel("Latitude")
        plt.title("Tsp Route")
        plt.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文标签
        plt.rcParams['axes.unicode_minus'] = False  # 这两行需要手动设置
        plt.grid(True)
        plt.legend()
        plt.show()
    
    
    # 获取最好的数据
    def get_elite(fitness, population):
        max_index = fitness.index(max(fitness))
        max_fitness = fitness[max_index]
        max_individual = population[max_index]
    
        return max_fitness, max_individual
    
    
    def record(f):
        global record_distance
        # 经纬度转千米的单位要乘以111
        record_distance.append(1 / f * 111)
    
    
    # 轮赌盘选择算子
    def selection(fitness, num):
        def select_one(fitness, fitness_sum):
            size = len(fitness)
            i = random.randint(0, size - 1)
            while True:
                if random.random() < fitness[i] / fitness_sum:
                    return i
                else:
                    i = (i + 1) % size
    
        res = set()
        fitness_sum = sum(fitness)
        while len(res) < num:
            t = select_one(fitness, fitness_sum)
            res.add(t)
        return res
    
    
    # 获得一个旅行路径的距离
    def get_distance(sequence):
        global distmat
    
        cost = 0
        for i in range(len(sequence)):
            cost += distmat[sequence[i - 1]][sequence[i]]
        return cost
    
    
    # 计算适应值
    def get_fitness(population):
        fitness = []
    
        for i in range(len(population)):
            fitness.append(1 / get_distance(population[i]))
    
        return fitness
    
    
    # 杂交算子
    def crossover(parent1, parent2):
        global individual_size
    
        a = random.randint(1, individual_size - 1)
        child1, child2 = parent1[:a], parent2[:a]
    
        for i in range(individual_size):
            if parent2[i] not in child1:
                child1.append(parent2[i])
    
            if parent1[i] not in child2:
                child2.append(parent1[i])
    
        return child1, child2
    
    
    # 初始化种群
    def init_population():
        global individual_size, population_size
    
        population_init = []
        for i in range(population_size):
            l = list(range(individual_size))
            population_init.append(random.sample(l, individual_size))
    
        return population_init
    
    
    # 获得城市之间的距离矩阵
    def get_distmat(M):
        length = M.shape[0]
        distmat = np.zeros((length, length))
    
        for i in range(length):
            for j in range(i + 1, length):
                distmat[i][j] = distmat[j][i] = np.linalg.norm(M[i] - M[j])
        return distmat
    
    
    if __name__ == "__main__":
        # 准备数据
        file = 'data1.csv'
        demo = 'demo.csv'
        coordinates = pd.read_csv(file, usecols=[1, 2], header=None).values
        name = pd.read_csv(demo, usecols=[0, 1, 2], encoding='gbk', header=None)
        distmat = get_distmat(coordinates)
    
        # 参数初始化
        individual_size = coordinates.shape[0]
        max_generation = 3500
        population_size = 10
        p_mutation = 0.2
        record_distance = []
    
        # 运行
        main()
    
    

    mutations.py

    # coding:utf:8
    import random
    
    
    # SMB
    def select_best_mutaion(s, distmat):
        s_res = [slide_mutation(s[:]), inversion_mutation(s[:]), irgibnnm_mutation(s[:], distmat)]
        res = [get_distance(s_res[0], distmat), get_distance(s_res[1], distmat), get_distance(s_res[2], distmat)]
    
        min_index = res.index(min(res))
    
        return s_res[min_index]
    
    
    # 滑动变异
    def slide_mutation(s):
        a, b = get_two_randint(len(s))
        t = s[a]
        for i in range(a + 1, b + 1):
            s[i - 1] = s[i]
        s[b] = t
        return s
    
    
    # 获得一个旅行路径的距离
    def get_distance(sequence, distmat):
        cost = 0
        for i in range(len(sequence)):
            cost += distmat[sequence[i - 1]][sequence[i]]
        return cost
    
    
    # 倒置变异
    def inversion_mutation(s):
        # 自己手写的2变换
        a, b = get_two_randint(len(s))
        for i in range(a, (a + b) // 2 + 1):
            s[i], s[b + a - i] = s[b + a - i], s[i]
    
        return s
    
    
    # 返回(小,大)两个随机数
    def get_two_randint(size):
        b = a = random.randint(0, size - 1)
        while a == b:
            b = random.randint(0, size - 1)
    
        if a > b:
            return b, a
        return a, b
    
    
    # irgibnnm
    def irgibnnm_mutation(s, distmat):
        a, b = get_two_randint(len(s))
        # 先inversion
        for i in range(a, (a + b) // 2 + 1):
            s[i], s[b + a - i] = s[b + a - i], s[i]
    
        # 再RGIBNNM
        b = (b + 1) % len(s)
    
        min = b - 1
        for i in range(len(s)):
            if i == b:
                continue
            if distmat[b][min] > distmat[b][i]:
                min = i
        s[b], s[min - 4] = s[min - 4], s[b]
    
        return s
    
    

    运行结果示例,只截图了一半
    在这里插入图片描述
    迭代次数与距离的关系
    在这里插入图片描述
    data1.csv

    Beijing,116.41667,39.91667
    shanghai,121.43333,34.5
    tianjin,117.2,39.13333
    guangzhou,113.23333,23.16667
    zhuhai,113.51667,22.3
    hangzhou,120.2,30.26667
    chongqing,106.45,29.5667
    qingdao,120.33333,36.03333
    fuzhou,119.3,26.08333
    lanzhou,103.73333,36.03333
    nanchang,115.9,28.68333
    nanjing,118.78333,32.05
    shenyang,123.38333,41.8
    zhengzhou,113.65,34.76667
    wuhan,114.31667,30.51667
    xian,108.95,34.26667
    changchun,125.35,43.8833
    haikou,110.35,20.01667
    guiyang,106.71667,26.56667
    changsha,113,28.2
    

    demo.csv

    116.41667,39.91667,北京
    121.43333,34.5,上海
    117.2,39.13333,天津
    113.23333,23.16667,广州
    113.51667,22.3,珠海
    120.2,30.26667,杭州
    106.45,29.5667,重庆
    120.33333,36.03333,青岛
    119.3,26.08333,福州
    103.73333,36.03333,兰州
    115.9,28.68333,南昌
    118.78333,32.05,南京
    123.38333,41.8,沈阳
    113.65,34.76667,郑州
    114.31667,30.51667,武汉
    108.95,34.26667,西安
    125.35,43.8833,长春
    110.35,20.01667,海口
    106.71667,26.56667,贵阳
    113,28.2,长沙
    

    如需更换城市找到相关城市的经纬度,修改此文件即可。

    展开全文
  • 网上有很多博客讲解遗传算法,但是大都只是“点到即止”,虽然给了一些代码实现,但也是“浅尝辄止”,没能很好地帮助大家进行扩展应用,抑或是进行深入的研究。 这是我的开篇之作~之前没有写博客的习惯,一般是将...

    加了个小目录~方便定位查看~

    前言

    正文

    一. 基础术语:

    二. 遗传算法基本算子:

    三.完整实现遗传算法:

    四.后记:


    前言

    网上有很多博客讲解遗传算法,但是大都只是“点到即止”,虽然给了一些代码实现,但也是“浅尝辄止”,没能很好地帮助大家进行扩展应用,抑或是进行深入的研究。

    这是我的开篇之作~之前没有写博客的习惯,一般是将笔记存本地,但久而久之发现回看不便,而且无法与大家交流和学习。现特此写下开篇之作,若有疏漏之处,敬请指正,谢谢!

    本文对遗传算法的原理进行梳理,相关代码是基于国内高校学生联合团队开源的高性能实用型进化算法工具箱geatpy来实现,部分教程引用了geatpy的官方文档:

    http://geatpy.com/index.php/details/

    geatpy官网:http://www.geatpy.com

    若有错误之处,恳请同志们指正和谅解,谢谢!

    因为是基于geatpy遗传和进化算法工具箱,所以下文的代码部分在执行前,需要安装geatpy:

    pip install geatpy

    安装时会自动根据系统版本匹配下载安装对应的版本。这里就有个小坑:如果最新版Geatpy没有与当前版本相匹配的包的话,会自动下载旧版的包。而旧版的包在Linux和Mac下均不可用。

    安装好后,在Python中输出版本检查是否是最新版(version 2.6.0):

    import geatpy as ea
    print(ea.__version__)
    

    下面切入主题:

    自然界生物在周而复始的繁衍中,基因的重组、变异等,使其不断具有新的性状,以适应复杂多变的环境,从而实现进化。遗传算法精简了这种复杂的遗传过程而抽象出一套数学模型,用较为简单的编码方式来表现复杂的现象,并通过简化的遗传过程来实现对复杂搜索空间的启发式搜索,最终能够在较大的概率下找到全局最优解,同时与生俱来地支持并行计算。

    下图展示了常规遗传算法(左侧) 和某种在并行计算下的遗传算法(右侧) 的流程。

    本文只研究经典的遗传算法,力求最后能够解决各种带约束单目标优化问题,并能够很轻松地进行扩展,让大家不仅学到算法理论,还能轻松地通过“复制粘贴”就能够将相关遗传算法代码结合到各类不同的现实问题的求解当中。

    从上面的遗传算法流程图可以直观地看出,遗传算法是有一套完整的“固定套路”的,我们可以把这个“套路”写成一个“算法模板”,即把:种群初始化、计算适应度、选择、重组、变异、生成新一代、记录并输出等等这些基本不需要怎么变的“套路”写在一个函数里面,而经常要变的部分:变量范围、遗传算法参数等写在这个函数外面,对于要求解的目标函数,由于在遗传进化的过程中需要进行调用目标函数进行计算,因此可以把目标函数、约束条件写在另一个函数里面。

    另外我们还可以发现,在遗传算法的“套路”里面,执行的“初始化种群”、“选择”、“重组”、“变异”等等,其实是一个一个的“算子”,在geatpy工具箱里,已经提供现行的多种多样的进化算子了,因此直接调用即可。

    Geatpy工具箱提供一个面向对象的进化算法框架,因此一个完整的遗传算法程序就可以写成这个样子:

    关于算法框架的详细介绍可参见:http://geatpy.com/index.php/geatpy%E6%95%99%E7%A8%8B/

    下面就来详细讲一下相关的理论和代码实现:

    正文

    一. 基础术语:

        先介绍一下遗传算法的几个基础的术语,分别是:”个体“、”种群“、”编码与解码“、”目标函数值“、”适应度值“。

    1.个体:“个体”其实是一个抽象的概念,与之有关的术语有:

    (1)个体染色体:即对决策变量编码后得到的行向量。

    比如:有两个决策变量x1=1,x2=2,各自用3位的二进制串去表示的话,写成染色体就是:

    (2)个体表现型:即对个体染色体进行解码后,得到的直接指代各个控制变量的值的行向量。

     比如对上面的染色体“0 0 1 0 1 0”进行解码,得到 “1 2”,它就是个体的表现型,可看出该个体存储着两个变量,值分别是1和2。

    注意:遗传算法中可以进行“实值编码”,即可以不用二进制编码,直接用变量的实际值来作为染色体。这个时候,个体的染色体数值上是等于个体的表现型的。

    (3)染色体区域描述器:用于规定染色体每一位元素范围,详细描述见下文。

    2. 种群:“种群”也是一个抽象的概念,与之有关的术语有:

    (1)种群染色体矩阵(Chrom):它每一行对应一个个体的染色体。此时会发出疑问:一个个体可以有多条染色体吗?答:可以有,但一般没有必要,一条染色体就可以存储很多决策变量的信息了,如果要用到多条染色体,可以用两个种群来表示。

    例如:

    它表示有3个个体(因为有3行),因此有3条染色体。此时需要注意:这些染色体代表决策变量的什么值,我们是不知道的,我们也不知道该种群的染色体采用的是什么编码。染色体具体代表了什么,取决于我们采用什么方式去解码。假如我们采用的是二进制的解码方式,并约定上述的种群染色体矩阵中前3列代表第一个决策变量,后3列代表第二个决策变量,那么,该种群染色体就可以解码成:

    (2)种群表现型矩阵(Phen):它每一行对应一个个体的表现型。比如上图就是根据Chrom种群染色体矩阵解码得到的种群表现型矩阵。同样地,当种群染色体采用的是“实值编码”时,种群染色体矩阵与表现型矩阵实际上是一样的。

    (3)种群个体违反约束程度矩阵(CV):它每一行对应一个个体,每一列对应一种约束条件(可以是等式约束或不等式约束)。CV矩阵中元素小于或等于0表示对应个体满足对应的约束条件,大于0则表示不满足,且越大表示违反约束条件的程度越高。比如有两个约束条件:

    如何计算CV矩阵?可以创建两个列向量CV1和CV2,然后把它们左右拼合而成一个CV矩阵。

    假设x1、x2、x3均为存储着种群所有个体的决策变量值的列向量(这里可以利用种群表现型矩阵Phen得到,比如x1=Phen[:, [0]];x2=Phen[:, [1]]);x3=Phen[:, [2]]),这样就可以得到种群所有个体对应的x1、x2和x3)。

    那么:

    比如在某一代中,种群表现型矩阵Phen为:

    则有:

    此时CV矩阵的值为:

    由此可见,第一个个体满足两个约束条件;第二个个体违反了2个约束条件;第三和第四个个体满足第一个约束条件但违反了第二个约束条件。

    下面看下如何用代码来生成一个种群染色体矩阵:

    代码1. 实整数值种群染色体矩阵的创建:

    import numpy as np
    from geatpy import crtpc
    help(crtpc) # 查看帮助
    # 定义种群规模(个体数目)
    Nind = 4
    Encoding = 'RI' # 表示采用“实整数编码”,即变量可以是连续的也可以是离散的
    # 创建“区域描述器”,表明有4个决策变量,范围分别是[-3.1, 4.2], [-2, 2],[0, 1],[3, 3],
    # FieldDR第三行[0,0,1,1]表示前两个决策变量是连续型的,后两个变量是离散型的
    FieldDR=np.array([[-3.1, -2, 0, 3],
                      [ 4.2,  2, 1, 5],
                      [ 0,    0, 1, 1]])
    # 调用crtri函数创建实数值种群
    Chrom=crtpc(Encoding, Nind, FieldDR)
    print(Chrom)
    

    代码1的运行结果:

    这里要插入讲一下“区域描述器”(见代码1中的FieldDR),它是用于描述种群染色体所表示的决策变量的一些信息,如变量范围、连续/离散性。另外还有一种区域描述器(FIeldD),用于描述二进制/格雷码的种群染色体。FieldDR和FieldD两个合称“Field”,又可以认为它们是“译码矩阵”。FieldDR具有以下的数据结构:

    代码1中的FieldDR矩阵的第三行即为这里的varTypes。它如果为0,表示对应的决策变量是连续型的变量;为1表示对应的是离散型变量。

    另一种则是用于二进制/格雷编码种群染色体解码的译码矩阵FieldD,它是具有以下结构的矩阵:

    其中,lens, lb, ub, codes, scales, lbin, ubin, varTypes均为长度等于决策变量个数的行向量。

    lens 包含染色体的每个子染色体的长度。sum(lens) 等于染色体长度。

    lb 和ub 分别代表每个决策变量的上界和下界。

    codes 指明染色体子串用的是二进制编码还是格雷编码。codes[i] = 0 表示第i 个变量使用的是标准二进制编码;codes[i] = 1 表示使用格雷编码。

    scales 指明每个子串用的是算术刻度还是对数刻度。scales[i] = 0 为算术刻度,scales[i] = 1 为对数刻度。对数刻度可以用于变量的范围较大而且不确定的情况,对于大范围的参数边界,对数刻度让搜索可用较少的位数,从而减少了遗传算法的计算量。(注意:当某个变量是对数刻度时,其取值范围中不能有0,即要么上下界都大于0,要么上下界都小于0。)

    从2.5.0版本开始,取消了对对数刻度的支持,该参数暂时保留,但不在起作用。

    lbin 和ubin 指明了变量是否包含其范围的边界。0 表示不包含边界;1 表示包含边界。

    varTypes 指明了决策变量的类型,元素为0 表示对应位置的决策变量是连续型变量;1 表示对应的是离散型变量。

    对于二进制编码,二进制种群的染色体具体代表决策变量的什么含义是不由染色体本身决定的,而是由解码方式决定的。因此在创建二进制种群染色体之初就要设定好译码矩阵(又称“区域描述器”)。

    因此,可以通过以下代码构建一个二进制种群染色体矩阵

    代码2. 二进制种群染色体矩阵的创建:

    import numpy as np
    from geatpy import crtpc
    help(crtpc) # 查看帮助
    # 定义种群规模(个体数目)
    Nind = 4
    Encoding = 'BG' # 表示采用“实整数编码”,即变量可以是连续的也可以是离散的
    # 创建“译码矩阵”
    FieldD = np.array([[3, 2], # 各决策变量编码后所占二进制位数,此时染色体长度为3+2=5
                       [0, 0], # 各决策变量的范围下界
                       [7, 3], # 各决策变量的范围上界
                       [0, 0], # 各决策变量采用什么编码方式(0为二进制编码,1为格雷编码)
                       [0, 0], # 各决策变量是否采用对数刻度(0为采用算术刻度)
                       [1, 1], # 各决策变量的范围是否包含下界(对bs2int实际无效,详见help(bs2int))
                       [1, 1], # 各决策变量的范围是否包含上界(对bs2int实际无效)
                       [0, 0]])# 表示两个决策变量都是连续型变量(0为连续1为离散)
    # 调用crtpc函数来根据编码方式和译码矩阵来创建种群染色体矩阵
    Chrom=crtpc(Encoding, Nind, FieldD)
    print(Chrom)
    

    代码2运行结果:

    3. 编码与解码

    对于实整数编码(即上面代码1所创建的实整数种群染色体),它是不需要解码,染色体直接就对应着它所代表的决策变量值。而对于代码2生成的二进制种群染色体矩阵,它需要根据译码矩阵FieldD来进行解码。在代码2后面添加以下语句即可解码:

    代码3(上接代码2):

    from geatpy import bs2ri
    help(bs2ri)
    Phen = bs2ri(Chrom, FieldD)
    print('表现型矩阵 = \n', Phen)

    运行效果如下:

    4.目标函数值

    种群的目标函数值存在一个矩阵里面(一般命名为ObjV),它每一行对应一个个体的目标函数值。对于单目标而言,这个目标函数值矩阵只有1列,而对于多目标而言,就有多列了,比如下面就是一个含两个目标的种群目标函数值矩阵:

    (这里Nind表示种群的规模,即种群含多少个个体;Nvar表示决策变量的个数)

    下面来跑一下代码,接着代码3,在对二进制染色体解码成整数值种群后,我们希望计算出f(x,y)=x+y这个目标函数值。同时设置一个等式约束:要求x + y = 3。于是完整代码如下:

    代码4:

    import numpy as np
    from geatpy import crtpc
    from geatpy import bs2ri
    
    def aim(Phen):
        x = Phen[:, [0]] # 取出种群表现型矩阵的第一列,得到所有个体的决策变量x
        y = Phen[:, [1]] # 取出种群表现型矩阵的第二列,得到所有个体的决策变量y
        CV = np.abs(x + y - 3) # 生成种群个体违反约束程度矩阵CV,以处理等式约束:x + y == 3
        f = x + y # 计算目标函数值
        return f, CV # 返回目标函数值矩阵
    
    # 定义种群规模(个体数目)
    Nind = 4
    Encoding = 'BG' # 表示采用“实整数编码”,即变量可以是连续的也可以是离散的
    # 创建“译码矩阵”
    FieldD = np.array([[3, 2], # 各决策变量编码后所占二进制位数,此时染色体长度为3+2=5
                       [0, 0], # 各决策变量的范围下界
                       [7, 3], # 各决策变量的范围上界
                       [0, 0], # 各决策变量采用什么编码方式(0为二进制编码,1为格雷编码)
                       [0, 0], # 各决策变量是否采用对数刻度(0为采用算术刻度)
                       [1, 1], # 各决策变量的范围是否包含下界(对bs2int实际无效,详见help(bs2int))
                       [1, 1], # 各决策变量的范围是否包含上界(对bs2int实际无效)
                       [0, 0]])# 表示两个决策变量都是连续型变量(0为连续1为离散)
    # 调用crtpc函数来根据编码方式和译码矩阵来创建种群染色体矩阵
    Chrom=crtpc(Encoding, Nind, FieldD)
    print('二进制染色体矩阵 = \n', Chrom)
    
    # 解码
    Phen = bs2ri(Chrom, FieldD)
    print('表现型矩阵 = \n', Phen)
    
    # 计算目标函数值矩阵
    ObjV, CV = aim(Phen)
    print('目标函数值矩阵 = \n', ObjV)
    print('CV矩阵 = \n', CV)
    

    运行结果如下:

    由上面对CV矩阵的描述可知,第三个个体的CV值为0,表示第三个个体满足x+y=3这个等式约束。其他都大于0,表示不满足该约束。

    疑问:CV矩阵有什么用呢?

    :CV矩阵既可用于标记非可行解,在含约束条件的优化问题中有用,又可用于度量种群个体违反各个约束条件的程度的高低。对于含约束条件的优化问题,我们可以采用罚函数或者是可行性法则来进行处理。罚函数法这里就不展开赘述了,最简单的罚函数可以是直接找到非可行解个体的索引,然后修改其对应的ObjV的目标函数值即可。

    对于可行性法则,它需要计算每个个体违反约束的程度,并把结果保存在种群类的CV矩阵中。CV矩阵的每一行对应一个个体、每一列对应一个约束条件(可以是等式约束也可以是不等式约束),CV矩阵中元素小于或等于0表示对应个体满足对应的约束条件,否则是违反对应的约束条件,大于0的值越大,表示违反约束的程度越高。生成CV标记之后,在后面调用适应度函数计算适应度时,只要把CV矩阵作为函数传入参数传进函数体,就会自动根据CV矩阵所描述的种群个体违反约束程度来计算出合适的种群个体适应度。

    5.适应度值

    适应度值通俗来说就是对种群个体的”生存能力的评价“。对于简单的单目标优化,我们可以简单地把目标函数值直接当作是适应度值(注意:当用geatpy遗传和进化算法工具箱时,则需要对目标函数值加个负号才能简单地把它当作适应度值,因为geatpy遵循的是”目标函数值越小,适应度值越大“的约定。)

    对于多目标优化,则需要根据“非支配排序”或是其他方法来确定种群个体的适应度值,本文不对其展开叙述。

    种群适应度(FitnV):它是一个列向量,每一行代表一个个体的适应度值:

    (这里Nind表示种群的规模,即种群含多少个个体)

    geatpy遗传和进化算法工具箱里面有几个函数可以计算种群个体的适应度 ,分别是:

    ranking、indexing、powing、scaling。具体的用法,可以用help命令查看,如help(ranking)。

    下面接着代码4,利用ranking(基于目标函数值排序的适应度分配)计算种群的适应度:

    代码5(接着代码4):

    from geatpy import ranking
    help(ranking)
    FitnV = ranking(ObjV, CV)
    print('种群适应度矩阵 = \n', FitnV)

    运行结果:

    分析这个结果我们发现,由于第1、2、4个体违反约束条件,而第三个个体满足约束条件,因此第3个个体的适应度最高。而在第1、2、4个体中,个体1的目标函数值最大,因此适应度最低。可见遵循“最小化目标”的约定,即目标函数值越小,适应度越大。

    好了,基本的术语和用法讲完后,下面讲一下遗传算法的基本算子:

     

    二. 遗传算法基本算子:

    我们不用破费时间去写底层的遗传算子,因为geatpy工具箱提供丰富的进化算子,以下所列算子不仅限于遗传算子:

     

     (如果是在iPython 控制台中调用可视化绘图函数(例如使用Spyder 开发工具),一般图像会默认显示在控制台或者是开发工具中。此时可以在iPython控制台下执行%matplotlib 来设置把图像显示在一个独立窗口中。)

    对于多目标优化,Geatpy中内置以下算子:

    可以用help(算子名)来查看对应的API文档,查看更详细的用法和例子。

    下面讲一下理论:

    1.选择:

    在进化算法中存在两个阶段的选择。第一次是参与进化操作的个体的选择。这个阶段的选择可以是基于个体适应度的、也可以是完全随机地选择交配个体。一旦个体被选中,那么它们就会参与交叉、变异等进化操作。未被选中的个体不会参与到进化操作中。

    第二次是常被称为“重插入”或“环境选择”的选择,它是指在个体经过交叉、变异等进化操作所形成的子代(或称“育种个体”)后用某种方法来保留到下一代从而形成新一代种群的过程。这个选择过程对应的是生物学中的” 自然选择”。它可以是显性地根据适应度(再次注意:适应度并不等价于目标函数值)来进行选择的,也可以是隐性地根据适应度(即不刻意去计算个体适应度)来选择。例如在多目标优化的NSGA-II 算法中,父代与子代合并后,处于帕累托分层中第一层级的个体以及处于临界层中的
    且拥挤距离最大的若干个个体被保留到下一代。这个过程就没有显性地去计算每个个体的适应度。

    经典的选择算子有:“轮盘赌选择”、“随机抽样选择”、“锦标赛选择”、“本地选择”、“截断选择”、“一对一生存者竞争选择”等等,这里不展开叙述了,可以参考:

    http://geatpy.com/index.php/2019/07/28/%E7%AC%AC%E5%9B%9B%E7%AB%A0%EF%BC%9A%E9%80%89%E6%8B%A9/

    这里要注意:遗传算法选择出的后代是可以有重复的。

    下面以低级选择函数:锦标赛选择算子(tour)为例,使用help(tour)查看其API,得到:

    实战演练如下:

    代码6:

    import numpy as np
    from geatpy import tour
    
    help(tour)
    FitnV = np.array([[1.2],[0.8],[2.1], [3.2],[0.6],[2.2],[1.7],[0.2]])
    chooseIdx = tour(FitnV, 6)
    print('个体的适应度为:\n', FitnV)
    print('选择出的个体的下标为:\n', chooseIdx)
    

    运行结果:

    光这样还不够,这里只是得出了选择个体的下标,如果我们需要得到被选中个体的染色体,同时尝试改用高级选择函数“selecting”来调用低级选择算子“tour”来进行选择,则可以如下操作:

    代码7:

    import numpy as np
    from geatpy import selecting
    
    help(selecting)
    Chrom=np.array([[1,11,21],
    [2,12,22],
    [3,13,23],
    [4,14,24],
    [5,15,25],
    [6,16,26],
    [7,17,27],
    [8,18,28]])
    
    FitnV = np.array([[1.2],[0.8],[2.1], [3.2],[0.6],[2.2],[1.7],[0.2]])
    SelCh = Chrom[selecting('tour', FitnV, 6), :] # 使用'tour'锦标赛选择算子,同时片取Chrom得到所选择个体的染色体
    print('个体的适应度为:\n', FitnV)
    print('选择后得到的种群染色矩阵为:\n', SelCh)
    

    代码7运行结果如下:

    将代码7中的'tour'换成工具箱中的其他选择算子的名称(如etour, rws, sus),就可以使用相应的选择算子进行选择。

    2.重组

    在很多的国内教材、博客文章、论文中,只提到“交叉”的概念。事实上,遗传算法的“交叉”是属于“重组”算子里面的。因为交叉算子经常使用,而对于“离散重组”、“中间重组”、“线性重组”等等,因为用的很少,所以我们常常只谈“交叉”算子了。交叉算子实际上是“值互换重组”(Values exchanged recombination)。这里就不展开叙述了,可以参考:

    http://geatpy.com/index.php/2019/07/28/%E7%AC%AC%E4%BA%94%E7%AB%A0%EF%BC%9A%E9%87%8D%E7%BB%84/

    与重组有关的遗传算法参数是“重组概率”(对于交叉而言就是“交叉概率”)(Pc),它是指两个个体的染色体之间发生重组的概率。

    下面以两点交叉(xovdp)为例,类似上面的“tour”那样我们使用help(xovdp)查看其API:

    代码实战演练如下:

    代码8:

    from geatpy import xovdp
    import numpy as np
    
    help(xovdp)
    OldChrom=np.array([[1,1,1,1,1],[1,1,1,1,1],[0,0,0,0,0],[0,0,0,0,0]]) #创建一个二进制种群染色体矩阵
    print('交叉前种群染色矩阵为:\n', OldChrom)
    NewChrom = xovdp(OldChrom, 1) # 设交叉概率为1
    print('交叉后种群染色矩阵为:\n', NewChrom)
    

    代码8运行结果如下:

    由此可见,xovdp是将前一半个体和后一半个体进行配对交叉的,有人认为应该随机选择交叉个体。事实上,在遗传算法进化过程中,各个位置的个体是什么,本身是随机的,不必要在交叉这里再增添一个随机(当然,可以在执行xovdp两点交叉之前,将种群染色体矩阵按行打乱顺序然后再交叉,从而满足自身需求)。

    3.变异

    下面以均匀突变(mutuni)为例,编写代码实现实数值编码的种群染色体的均匀突变,使用help(mutuni)查看API文档。

    代码9:

    from mutuni import mutuni
    import numpy as np
    # 自定义种群染色体矩阵,表示有3个个体,且染色体元素直接表示变量的值(即实值编码)
    OldChrom = np.array([[9,10],
                         [10,10],
                         [10,10]])
    # 创建区域描述器(又称译码矩阵)
    FieldDR = np.array([[7,8],
                        [10,10],
                        [1, 1]])
    # 此处设编码方式为实值编码中的“实整数编码”RI,表示染色体可代表实数和整数
    NewChrom = mutuni('RI', OldChrom, FieldDR, 1)
    print(NewChrom)
    

    代码9的运行结果如下:

    4.重插入

    经典遗传算法通过选择、重组和变异后,我们得到的是育种后代,此时育种后代的个体数有可能会跟父代种群的个体数不相同。这时,为了保持种群的规模,这些育种后代可以重新插入到父代中,替换父代种群的一部分个体,或者丢弃一部分育种个体,最终形成新一代种群。前面曾提到过“重插入”也是一种选择,但它是环境选择,是用于生成新一代种群的;而前面在交叉变异之前的选择是用于选择个体参与交叉变异,那个选择常被称作“抽样”。

    现考虑使用精英个体保留的遗传算法“EGA”,则重插入操作如下:

    代码10:

    from mutuni import mutuni
    import numpy as np
    # 自定义父代种群染色体(仅作为例子):
    Chrom = np.array([[1.1, 1.3],
                      [2.4, 1.2],
                      [3,   2.1],
                      [4,   3.1]])
    # 若父代个体的适应度为:
    FitnV = np.array([[1],
                     [2],
                     [3],
                     [4]])
    # 考虑采用“精英保留策略”的遗传算法,此时从父代选择出4-1=3个个体,经过交叉变异后假设子代的染色体为:
    offChrom = np.array([[2.1, 2.3],
                         [2.3, 2.2],
                         [3.4, 1.1]])
    # 假设直接把目标函数值当作适应度,且认为适应度越大越好。则通过以下代码重插入生成新一代种群:
    bestIdx = np.argmax(FitnV) # 得到父代精英个体的索引
    NewChrom = np.vstack([Chrom[bestIdx, :], offChrom]) # 得到新一代种群的染色体矩阵
    print('新一代种群的染色体矩阵为:\n', NewChrom)
    

    在“EGA”中,假设父代种群规模为N,则选择出(N-1)个个体进行交叉变异,然后找出父代的精英个体,用着个精英个体和交叉变异得到的子代个体进行“拼合”,得到新一代种群。

    除了这种重插入生成新一代种群的方法外,还有“父子两代个体合并选择”等更优秀的生成新一代种群的方法,这里就不一一赘述了。

    讲完上面的基本术语以及遗传算法基本算子后,我们就可以来利用遗传算法的“套路”编写一个遗传算法求解问题的程序了:

     

    三.完整实现遗传算法:

    上文提到遗传算法程序可以写成这个样子:

    其核心是算法模板类。在遗传算法模板里,我们根据遗传算法的“套路”,进行:初始化种群、目标函数值计算、适应度评价、选择、重组、变异、记录各代最优个体等操作。geatpy工具箱内置有开源的“套路模板”,源代码参见:

    https://github.com/geatpy-dev/geatpy/tree/master/geatpy/source-code/templets

    这些内置算法模板有详细的输入输出参数说明,以及遗传算法整个流程的完整注释,它们可以应对简单或复杂的、单目标优化的、多目标优化的、约束优化的、组合优化的等等的问题。

    但为了学习,我这里先不采用框架,直接利用工具箱提供的库函数来写一个带精英个体保留的遗传算法。这样代码量比较大,但有利于入门。

    编写代码 11、12,分别放在同一个文件夹下:

    代码11(目标函数aimfuc.py)(这里要回顾一下前面,Phen是种群表现型矩阵,存储的是种群所有个体的表现型,而不是单个个体。因而计算得到的目标函数值矩阵也是包含所有个体的目标函数值):

    # -*- coding: utf-8 -*-
    """
    aimfunc.py - 目标函数文件
    描述:
        目标:max f = 21.5 + x1 * np.sin(4 * np.pi * x1) + x2 * np.sin(20 * np.pi * x2)
        约束条件:
            x1 != 10
            x2 != 5
            x1 ∈ [-3, 12.1] # 变量范围是写在遗传算法的参数设置里面
            x2 ∈ [4.1, 5.8]
    """
    
    import numpy as np
    
    def aimfunc(Phen, CV):
        x1 = Phen[:, [0]] # 获取表现型矩阵的第一列,得到所有个体的x1的值
        x2 = Phen[:, [1]]
        f = 21.5 + x1 * np.sin(4 * np.pi * x1) + x2 * np.sin(20 * np.pi * x2)
        exIdx1 = np.where(x1 == 10)[0] # 因为约束条件之一是x1不能为10,这里把x1等于10的个体找到
        exIdx2 = np.where(x2 == 5)[0]
        CV[exIdx1] = 1
        CV[exIdx2] = 1
        return [f, CV]
    

    然后是编写算法:

    代码12:

    # -*- coding: utf-8 -*-
    """main.py"""
    import numpy as np
    import geatpy as ea # 导入geatpy库
    from aimfunc import aimfunc # 导入自定义的目标函数
    import time
    
    """============================变量设置============================"""
    x1 = [-3, 12.1] # 第一个决策变量范围
    x2 = [4.1, 5.8] # 第二个决策变量范围
    b1 = [1, 1] # 第一个决策变量边界,1表示包含范围的边界,0表示不包含
    b2 = [1, 1] # 第二个决策变量边界,1表示包含范围的边界,0表示不包含
    ranges=np.vstack([x1, x2]).T # 生成自变量的范围矩阵,使得第一行为所有决策变量的下界,第二行为上界
    borders=np.vstack([b1, b2]).T # 生成自变量的边界矩阵
    varTypes = np.array([0, 0]) # 决策变量的类型,0表示连续,1表示离散
    """==========================染色体编码设置========================="""
    Encoding = 'BG' # 'BG'表示采用二进制/格雷编码
    codes = [0, 0] # 决策变量的编码方式,设置两个0表示两个决策变量均使用二进制编码
    precisions =[4, 4] # 决策变量的编码精度,表示二进制编码串解码后能表示的决策变量的精度可达到小数点后6位
    scales = [0, 0] # 0表示采用算术刻度,1表示采用对数刻度
    FieldD = ea.crtfld(Encoding,varTypes,ranges,borders,precisions,codes,scales) # 调用函数创建译码矩阵
    """=========================遗传算法参数设置========================"""
    NIND      = 100; # 种群个体数目
    MAXGEN    = 200; # 最大遗传代数
    maxormins = [-1] # 列表元素为1则表示对应的目标函数是最小化,元素为-1则表示对应的目标函数是最大化
    selectStyle = 'rws' # 采用轮盘赌选择
    recStyle  = 'xovdp' # 采用两点交叉
    mutStyle  = 'mutbin' # 采用二进制染色体的变异算子
    pc        = 0.7 # 交叉概率
    pm        = 1 # 整条染色体的变异概率(每一位的变异概率=pm/染色体长度)
    Lind = int(np.sum(FieldD[0, :])) # 计算染色体长度
    obj_trace = np.zeros((MAXGEN, 2)) # 定义目标函数值记录器
    var_trace = np.zeros((MAXGEN, Lind)) # 染色体记录器,记录历代最优个体的染色体
    """=========================开始遗传算法进化========================"""
    start_time = time.time() # 开始计时
    Chrom = ea.crtpc(Encoding, NIND, FieldD) # 生成种群染色体矩阵
    variable = ea.bs2ri(Chrom, FieldD) # 对初始种群进行解码
    CV = np.zeros((NIND, 1)) # 初始化一个CV矩阵(此时因为未确定个体是否满足约束条件,因此初始化元素为0,暂认为所有个体是可行解个体)
    ObjV, CV = aimfunc(variable, CV) # 计算初始种群个体的目标函数值
    FitnV = ea.ranking(maxormins * ObjV, CV) # 根据目标函数大小分配适应度值
    best_ind = np.argmax(FitnV) # 计算当代最优个体的序号
    # 开始进化
    for gen in range(MAXGEN):
        SelCh = Chrom[ea.selecting(selectStyle,FitnV,NIND-1),:] # 选择
        SelCh = ea.recombin(recStyle, SelCh, pc) # 重组
        SelCh = ea.mutate(mutStyle, Encoding, SelCh, pm) # 变异
        # 把父代精英个体与子代的染色体进行合并,得到新一代种群
        Chrom = np.vstack([Chrom[best_ind, :].astype(int), SelCh])
        Phen = ea.bs2ri(Chrom, FieldD) # 对种群进行解码(二进制转十进制)
        ObjV, CV = aimfunc(Phen, CV) # 求种群个体的目标函数值
        FitnV = ea.ranking(maxormins * ObjV, CV) # 根据目标函数大小分配适应度值
        # 记录
        best_ind = np.argmax(FitnV) # 计算当代最优个体的序号
        obj_trace[gen,0]=np.sum(ObjV)/ObjV.shape[0] #记录当代种群的目标函数均值
        obj_trace[gen,1]=ObjV[best_ind] #记录当代种群最优个体目标函数值
        var_trace[gen,:]=Chrom[best_ind,:] #记录当代种群最优个体的染色体
    # 进化完成
    end_time = time.time() # 结束计时
    ea.trcplot(obj_trace, [['种群个体平均目标函数值', '种群最优个体目标函数值']]) # 绘制图像
    """============================输出结果============================"""
    best_gen = np.argmax(obj_trace[:, [1]])
    print('最优解的目标函数值:', obj_trace[best_gen, 1])
    variable = ea.bs2ri(var_trace[[best_gen], :], FieldD) # 解码得到表现型(即对应的决策变量值)
    print('最优解的决策变量值为:')
    for i in range(variable.shape[1]):
        print('x'+str(i)+'=',variable[0, i])
    print('用时:', end_time - start_time, '秒')

    执行代码12得到如下结果:

    终于,我们把遗传算法完整地实现了,但扩展性还不够高。下面学习下如何使用Geatpy提供的进化算法框架来求解上述问题:(关于使用框架来优化的介绍可详见http://geatpy.com/index.php/geatpy%E6%95%99%E7%A8%8B/

    在这里我们可以回顾以下在本文开头提到的采用遗传算法的“套路”来编程求解问题的基本流程:

    其中执行脚本和问题类是需要编写的,算法模板类我直接调用Geatpy内置的"soea_EGA_templet"(带精英个体保留的单目标遗传算法模板)。下面开始编写代码:

    代码13:问题类"MyProblem"的编写:

    # -*- coding: utf-8 -*-
    """
    MyProblem.py
    该案例展示了一个简单的连续型决策变量最大化目标的单目标优化问题。
    max f = x * np.sin(10 * np.pi * x) + 2.0
    s.t.
    -1 <= x <= 2
    """
    
    import numpy as np
    import geatpy as ea
    
    class MyProblem(ea.Problem): # 继承Problem父类
        def __init__(self):
            name = 'MyProblem' # 初始化name(函数名称,可以随意设置)
            M = 1 # 初始化M(目标维数)
            maxormins = [-1] # 初始化maxormins(目标最小最大化标记列表,1:最小化该目标;-1:最大化该目标)
            Dim = 2 # 初始化Dim(决策变量维数)
            varTypes = [0] * Dim # 初始化varTypes(决策变量的类型,元素为0表示对应的变量是连续的;1表示是离散的)
            lb = [-3, 4.1] # 决策变量下界
            ub = [12.1, 5.8] # 决策变量上界
            lbin = [1] * Dim # 决策变量下边界
            ubin = [1] * Dim # 决策变量上边界
            # 调用父类构造方法完成实例化
            ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)
        
        def aimFunc(self, pop): # 目标函数
            x1 = pop.Phen[:, [0]] # 获取表现型矩阵的第一列,得到所有个体的x1的值
            x2 = pop.Phen[:, [1]]
            f = 21.5 + x1 * np.sin(4 * np.pi * x1) + x2 * np.sin(20 * np.pi * x2)
            exIdx1 = np.where(x1 == 10)[0] # 因为约束条件之一是x1不能为10,这里把x1等于10的个体找到
            exIdx2 = np.where(x2 == 5)[0]
            pop.CV = np.zeros((pop.sizes, 2))
            pop.CV[exIdx1, 0] = 1
            pop.CV[exIdx2, 1] = 1
            pop.ObjV = f # 计算目标函数值,赋值给pop种群对象的ObjV属性
        

    第二步:编写执行脚本“main.py”

    # -*- coding: utf-8 -*-
    """main.py"""
    import geatpy as ea # import geatpy
    from MyProblem import MyProblem # 导入自定义问题接口
    """===============================实例化问题对象================================"""
    problem = MyProblem() # 生成问题对象
    """==================================种群设置=================================="""
    Encoding = 'BG'       # 编码方式
    NIND = 100            # 种群规模
    Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges, problem.borders) # 创建区域描述器
    population = ea.Population(Encoding, Field, NIND) # 实例化种群对象(此时种群还没被初始化,仅仅是完成种群对象的实例化)
    """================================算法参数设置================================="""
    myAlgorithm = ea.soea_EGA_templet(problem, population) # 实例化一个算法模板对象
    myAlgorithm.MAXGEN = 200 # 最大进化代数
    myAlgorithm.logTras = 1  # 设置每隔多少代记录日志,若设置成0则表示不记录日志
    myAlgorithm.verbose = True  # 设置是否打印输出日志信息
    myAlgorithm.drawing = 1  # 设置绘图方式(0:不绘图;1:绘制结果图;2:绘制目标空间过程动画;3:绘制决策空间过程动画)
    """===========================调用算法模板进行种群进化==============--==========="""
    [BestIndi, population] = myAlgorithm.run()  # 执行算法模板,得到最优个体以及最后一代种群
    BestIndi.save()  # 把最优个体的信息保存到文件中
    """==================================输出结果=================================="""
    print('用时:%f 秒' % myAlgorithm.passTime)
    print('评价次数:%d 次' % myAlgorithm.evalsNum)
    if BestIndi.sizes != 0:
        print('最优的目标函数值为:%s' % BestIndi.ObjV[0][0])
        print('最优的控制变量值为:')
        for i in range(BestIndi.Phen.shape[1]):
            print(BestIndi.Phen[0, i])
    else:
        print('没找到可行解。')

    运行"main.py"执行脚本即可得到以下结果:

                                                                                                         ......

    代码解析:在“main.py”执行脚本中,一开始需要实例化一个问题对象。然后是种群对象的实例化。在实例化种群对象前,需要设定种群的编码方式Encoding、种群规模NIND,并且生成区域描述器Field(或称译码矩阵),因为种群类的构造方法中需要至少用到这三个参数(详见“Population.py”中种群类的构造方法)。

    在完成了问题类对象和种群对象的实例化后,将其传入算法模板类的构造方法来实例化一个算法模板对象。这里我实例化的是单目标优化的EGA算法(即带精英个体保留的遗传算法)的模板类对象,即代码中的"soea_EGA_templet"。里面的进化算法具体是如何操作的,可详见https://github.com/geatpy-dev/geatpy/blob/master/geatpy/templates/soeas/GA/EGA/soea_EGA_templet.py

    采用Geatpy提供的进化算法框架可以既能最大程度地描述清楚所要求解的问题,而且与进化算法是高度脱耦的,即上面在编写问题类的时候完全不需要管后面采用什么算法、采用什么样编码的种群,只需把问题描述清楚即可。

    而且,遗传算法有个好处是:目标函数可以写得相当复杂,可以解决各种复杂的问题,比如神经网络。以BP神经网络为例,可以把神经网络的参数作为决策变量,神经网络的训练误差作为目标函数值,只需把上面的例子修改一下就行了。

    而且,一般而言我们不需要像我刚刚最开始那样刻意去手写进化算法,可以直接调用geatpy内置的算法模板就可以解决问题了。geatpy工具箱提供这些内置的算法模板:

    应用Geatpy求解数学建模、工业设计、资源调度等实际优化问题的的朋友们可以直接使用这些算法模板快速解决各种灵活的优化问题。

    四.后记:

    最后十分感谢由Geatpy团队提供的高性能实用型遗传和进化算法工具箱,它提供开源的进化算法框架为遗传算法求解单目标/多目标优化、约束优化、组合优化等等给出了相当准确和快捷的解决方案。据称,geatpy的运行性能相当的高,远高于matlab的遗传算法工具箱、以及采用JAVA、matlab或者Python编写的一些进化优化平台或框架,比如jMetal、platemo、pymoo、deap等,后面有时间我将进行详细的性能对比实验分析,有相关经验的读者也可以自行对比性能。而且依我的体验来看,这是我网上到处找代码、找资料学习、碰了无数次壁后所看到的比较易学易用的工具箱了。

    最后值得注意的是Geatpy的顶层是面向对象的进化算法框架,底层是面向过程的进化算子。下面放一张geatpy的UML类图、算法调用的层次结构和库函数调用关系图,以此记录方便查看:

    下面附一张一位同行朋友使用犀牛软件(Rhinoceros)结合geatpy工具箱进行产品优化设计的截图:

    很多工程软件都提供Python接口,当需要用到进化优化时,就可以编写Python代码进行优化了。

    下一篇博客将介绍如何用遗传算法求解有向图的最短路径问题:

    https://blog.csdn.net/weixin_37790882/article/details/100622338

    后续我将继续学习和挖掘该工具箱的更多深入的用法。希望这篇文章在帮助自己记录学习点滴之余,也能帮助大家!

     

    展开全文
  • 基于遗传算法的移动机器人路径规划

    万次阅读 多人点赞 2019-02-10 18:09:40
    之前在网上找基于遗传算法的移动机器人路径规划资料的时候,没有找到理想的开源代码。最后参照论文,用matlab写了代码。平时学习的时候,会看很多人写的一些博客和公众号文章,觉得很有帮助,所以也想分享一点自己学...
  • 模拟退火算法路径规划(python代码实现)一、引言二、算法伪代码三、算法流程及...  本方法与遗传算法有类似之处,对比上一篇博客可以加深理解:算法学习之遗传算法路径规划python代码实现) 二、算法伪代码   
  • 03-0003 Python遗传算法解决旅行商问题

    千次阅读 2019-04-20 21:27:57
    Python遗传算法解决旅行商问题0.前言1.介绍1.1算法介绍1.2问题介绍2.流程2.1流程图2.2选择过程2.2.1轮盘赌选择法2.2.2比例选择法2.2.3精英保留策略2.3交叉过程2.3变异3.源码4.结果5.总结 0.前言 物竞天择,适者生存...
  • Python编写的,基于遗传算法的,旅行商路径优化。通过遗传算法迭代的方式,得到旅行商路径问题的优化解。
  • 文章将基于遗传算法的原理和传统求解步骤依据具体的TSP问题做出优化改进求解51个城市最短路径规划问题,并借助python语言实现交互功能(用户可从51个城市中自行选择旅游城市,程序将为用户推荐最佳旅行方案)。
  • 基于遗传算法的多机器人栅格路径规划,完成无碰撞路径的规划!
  • 遗传算方法求解最短路径问题1.个体编码2.代码实现 一.前言 在优化问题中,网络模型是很重要的一类问题,各种物流配送计划、供应链管理、公路网络设计等等问题都可以简化为网络模型。从这里开始,我们会进入基本网络...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,412
精华内容 564
关键字:

python遗传算法路径规划

python 订阅