遗传算法原理与实战（python、matlab）

遗传算法

---

1.什么是遗传算法

遗传算法（Genetic Algorithm，简称GA）是一种基于生物进化论和遗传学原理的全局优化搜索算法。它通过模拟自然界中生物种群的遗传机制和进化过程来解决复杂问题，如函数优化、组合优化、机器学习等。遗传算法具有并行性、全局搜索能力和对问题描述的简单性，在很多领域有着广泛应用。

2.遗传算法的基本原理

遗传算法的基本思想是将问题的解表示为“染色体”（Chromosome），每个染色体由一系列“基因”（Gene）组成，这些基因可以是二进制编码、实数或离散值。算法通过选择、交叉（也称为杂交）、变异等遗传操作，模仿自然选择和遗传变异的过程，使种群中的个体逐渐适应环境，即找到问题的最优解或近似最优解。

下面我试着用通俗一点的例子来解释什么是“遗传算法”：

想象一下，在一座起伏的山脉中，住着一群人。这些人散居在山脉的不同高度上，有的住在山顶，有的住在半山腰，还有的住在山谷里。突然有一天，洪水来袭，水位不断上升，那些住在低处的人们开始面临危险。但住在高处的人因为地势较高，所以能够幸免于难，我们可以说，这些人对这种“洪水”环境的适应能力更强。

为了确保人类的延续，这些适应能力强的人开始繁衍后代。他们的目标是希望通过繁殖，能够产生出生活在山脉最高点的后代，这样就能完全避免洪水的威胁，相当于找到了最完美的生存策略。在这个过程中，父母双方的特征（我们称之为“基因”）会通过一种类似“混合”的方式传递给子女，有时候，子女身上还会出现一些意料之外的新特征，这就是“基因变异”。

遗传算法就像是指导这群人在面对洪水时如何繁衍出能在最高山峰生存的后代的一种智慧。它通过模拟自然界中的选择、交配和变异过程，帮助寻找出最适应环境的解决方案。在这个例子中，“最高山峰”的横坐标就是我们要找的答案，也就是那组能让人们在任何情况下都能生存下来的最优“DNA”。

即变成了求解该函数在区间[0 , 5]上的最大值点的问题。

简而言之，遗传算法就像是一种搜索策略，它通过模仿自然进化的过程，逐步改进并找到解决问题的最佳方案。在这个过程中，那些表现更好（适应度更高）的解决方案会被优先保留和进一步优化，直到找到最优解。

3.遗传算法的名词解释

种群（Population）：遗传算法中的一组可能解的集合。每个解都是一个染色体。

染色体（Chromosome）：代表一个可能的解，由多个基因组成。

基因（Gene）：染色体上的基本单元，通常是一个二进制位或数值，表示解的一个部分。

适应度函数（Fitness Function）：用于评估染色体的好坏，即解的质量。适应度高的染色体有更大的机会被选中进行遗传操作。

选择（Selection）：根据适应度函数，从当前种群中选择优秀的个体进入下一代种群的过程。常用的选择方法有轮盘赌选择、锦标赛选择等。

交叉（Crossover）：两个染色体交换部分基因，生成新的染色体。这是遗传算法中实现全局搜索的主要方式。

变异（Mutation）：以很小的概率改变染色体中的某个基因，增加种群的多样性，防止算法过早收敛。

收敛（Convergence）：当种群中的个体适应度不再显著提高时，认为算法已经找到最优解或近似最优解，此时算法结束。

遗传算法通过迭代上述过程，不断优化种群，最终找到满足要求的解。

4.遗传算法程序实现过程

遗传算法实现的一般流程如下：

（1）初始化种群
生成一个初始种群，种群中的每一个个体（染色体）都代表了一个可能的解。
每个个体由一组基因构成，基因的编码形式取决于问题的特性，可以是二进制、实数或符号编码等。
（2）评估适应度
使用适应度函数（Fitness Function）来评估每个个体的适应度，这反映了个体在解决特定问题上的性能。
适应度函数应该设计得能够量化个体解的质量，以便于比较不同个体的优劣。
（3）选择（Selection）
根据个体的适应度值进行选择，适应度较高的个体有更高的概率被选中参与后续的遗传操作。
常见的选择策略包括轮盘赌选择、锦标赛选择、排名选择等。
（4）遗传操作
交叉（Crossover）：选定的个体之间进行基因信息的交换，以产生新的后代个体。
变异（Mutation）：以小概率随机改变后代个体的某些基因，以增加种群的多样性。
这些操作通常遵循一定的概率规则，例如交叉概率和变异概率。
（5）替换（Replacement）
新产生的后代个体可能会替代旧的种群成员，这个过程可以通过多种策略进行，如世代替换、稳态遗传算法等。
有时会使用精英策略（Elitism），确保种群中最好的个体被保留到下一代。
（6）终止条件检查
检查是否达到预定的终止条件，这可能是固定的迭代次数、种群的平均适应度或最佳适应度达到某个阈值。
如果满足终止条件，则算法停止；否则，继续回到步骤3。
（7）输出结果
算法结束后，输出种群中适应度最高的个体作为最优解。

5.遗传算法代码实现

5.1python代码实现

种群类设计：

import numpy as np
# 定义种群类
class Population:
 
    def __init__(self, pop_size, chromosome_size):
        # 初始化种群大小和染色体大小
        self.pop_size = pop_size
        self.chromosome_size = chromosome_size
        # 随机生成初始种群，每个基因是0或1
        self.population = np.round(np.random.rand(pop_size, chromosome_size)).astype(np.int)
        # 初始化适应度值为0
        self.fit_value = np.zeros((pop_size, 1))
 
    # 选择染色体方法
    def select_chromosome(self):
        # 计算总适应度值
        total_fitness_value = self.fit_value.sum()
        # 计算每个个体的选择概率
        p_fit_value = self.fit_value / total_fitness_value
        # 计算累积概率
        p_fit_value = np.cumsum(p_fit_value)
        # 生成随机点，用于轮盘赌选择
        point = np.sort(np.random.rand(self.pop_size, 1), 0)
        # 初始化选择和新的种群索引
        fit_in = 0
        new_in = 0
        new_population = np.zeros_like(self.population)
        # 轮盘赌选择新种群
        while new_in < self.pop_size:
            if point[new_in] < p_fit_value[fit_in]:
                new_population[new_in, :] = self.population[fit_in, :]
                new_in += 1
            else:
                fit_in += 1
        # 更新种群
        self.population = new_population
 
    # 交叉染色体方法
    def cross_chromosome(self, cross_rate):
        # 获取种群大小和染色体大小
        x = self.pop_size
        y = self.chromosome_size
        new_population = np.zeros_like(self.population)
        # 对种群中的染色体进行交叉操作
        for i in range(0, x-1, 2):
            # 根据交叉率决定是否交叉
            if np.random.rand(1) < cross_rate:
                # 生成交叉点
                insert_point = int(np.round(np.random.rand(1) * y).item())
                # 交叉染色体
                new_population[i, :] = np.concatenate([self.population[i, 0:insert_point], self.population[i+1, insert_point:y]], 0)
                new_population[i+1, :] = np.concatenate([self.population[i+1, 0:insert_point], self.population[i, insert_point:y]], 0)
            else:
                # 不交叉则直接复制
                new_population[i, :] = self.population[i, :]
                new_population[i + 1, :] = self.population[i + 1, :]
        # 更新种群
        self.population = new_population
 
    # 寻找最佳个体方法
    def best(self):
        # 初始化最佳个体和适应度值
        best_individual = self.population[0, :]
        best_fit = self.fit_value[0]
        # 遍历种群寻找最佳个体
        for i in range(1, self.pop_size):
            if self.fit_value[i] > best_fit:
                best_individual = self.population[i, :]
                best_fit = self.fit_value[i]
        return best_individual, best_fit
 
    # 突变染色体方法
    def mutation_chromosome(self, mutation_rate):
        # 获取种群大小
        x = self.pop_size
        # 对种群中的染色体进行突变操作
        for i in range(x):
            if np.random.rand(1) < mutation_rate:
                # 生成突变点
                m_point = int(np.round(np.random.rand(1) * self.chromosome_size).item())
                # 突变操作
                if self.population[i, m_point] == 1:
                    self.population[i, m_point] = 0
                else:
                    self.population[i, m_point] = 1
 
    # 二进制转十进制方法
    def binary2decimal(self, population):
        pop1 = np.zeros_like(population)
        y = self.chromosome_size
        # 转换二进制到十进制
        for i in range(y):
            pop1[:, i] = 2 ** (y - i - 1) * population[:, i]
        pop = np.sum(pop1, 1)
        pop2 = pop * 10 / (1 << y)
        return pop2
 
    # 计算目标函数值方法
    def cal_obj_value(self):
        # 将二进制表示转换为十进制表示
        x = self.binary2decimal(self.population)
        # 计算并更新适应度值
        self.fit_value = 10 * np.sin(5 * x) + 7 * np.abs(x - 5) + 

*测试：*

from GA遗传算法.Population import Population
import numpy as np
 
cross_rate = 0.6
mutation_rate = 0.001
pop_size = 100
chromosome_size = 10
population = Population(100, 10)
 
for i in range(100):
    population.cal_obj_value()
    population.select_chromosome()
    population.cross_chromosome(cross_rate)
    population.mutation_chromosome(mutation_rate)
    best_individual, best_fit = population.best()
    best_individual = np.expand_dims(best_individual, 0)
    x = population.binary2decimal(best_individual)
    print("X：", x, "\t Y: ", best_fit)

5.3结果展示

在某些情况下遗传算法可能陷入局部最优。

1. 早熟收敛（Premature Convergence）
   原因：种群中的个体过早地变得相似，导致遗传多样性下降，使得算法无法探索解空间的其他区域。
   表现：种群的适应度值在算法早期迅速上升，但随后停滞不前。
2. 选择压力过大
   原因：选择操作过于偏向适应度高的个体，可能导致适应度较低的个体被过早淘汰，减少了种群的多样性。
   表现：算法快速收敛到某个局部最优解，而忽略了可能存在更好解的区域。
3. 交叉和突变概率设置不当
   原因：交叉率和突变率太低，无法有效地探索新的解空间；交叉率和突变率太高，可能导致算法变成随机搜索。
   表现：交叉和突变不足时，算法容易在局部最优解附近徘徊；过度交叉和突变则可能导致算法无法收敛。
4. 适应度函数设计不当
   原因：适应度函数未能准确反映解的质量，或者存在多个局部最优解，而适应度函数无法有效区分。
   表现：算法可能收敛到非目标解的局部最优解。
5. 解的表示方法
   原因：解的编码方式可能不适合特定问题，导致算法难以找到全局最优解。
   表现：即使算法能够探索解空间，但由于编码的限制，仍可能无法表示全局最优解。
6. 种群大小
   原因：种群大小太小，无法维持足够的遗传多样性；种群大小太大，计算成本过高。
   表现：种群太小容易陷入局部最优，种群太大则可能导致算法效率低下。
7. 算法停止条件
   原因：算法停止条件设置过早，可能导致算法未能充分探索解空间。
   表现：算法在达到停止条件时可能刚好处于局部最优解。

遗传算法的MATLAB代码

---

%ga_1
clc
clear
close all

%% 画出函数图
figure(1);
hold on;

lb = 1;  ub= 2;                                                   %函数自变量范围【1.2】
ezplot('sin(10*pi*X)/X',[lb,ub]);                 %画出函数曲线
xlabel('自变量/x')
ylabel('自变量/y')

%% 定义遗传算法参数
NIND = 40;                                  %种群大小
MAXGEN = 20;                            %最大遗传代数
PRECI = 20;                                 %个体长度
GGAP = 0.95;                                %代沟
px = 0.7;                                          %交叉概率
pm = 0.01;                                      %变异概率
trace = zeros(2,MAXGEN);       %寻优结果的初始值
FieldD = [PRECI;lb;ub;1;0;1;1];  %区域描述器
Chrom = crtbp(NIND,PRECI);  %创建任意离散随机种群

%% 优化
gen = 0;                                            %代计数器
X = bs2rv(Chrom,FieldD);            %初始种群二进制到十进制转换
ObjV = sin(10*pi*X)./X;                 %计算目标函数值
while gen<MAXGEN
    FitnV = ranking(ObjV);              %分配适应度值
    SelCh = select('sus',Chrom,FitnV,GGAP);               %选择
    SelCh = recombin('xovsp',SelCh,px);                     %重组
    SelCh = mut(SelCh,pm);                                      %变异
    X = bs2rv(SelCh,FieldD);                                    %自带个体的十进制转换
    ObjVSel = sin(10*pi*X)./X;                              %计算自带的目标函数值
    [Chrom,ObjV] = reins(Chrom,SelCh,1,1,ObjV,ObjVSel);         %重插入子代到父代，得到新种群
    X = bs2rv(Chrom,FieldD);
    gen = gen+1;                                         % 代计数器增加
    
    %获取每代的最优解及其序号，Y为最优解，I为个体的序号
    [Y,I] = min(ObjV);
    trace(1,gen) = X(I);                         %记下每代的最优值
    trace(2,gen) = Y;                               %记下每代的最优值
end

plot(trace(1,:),trace(2,:),'bo');            %画出每代的最优点
grid on;
plot(X,ObjV,'b*');                                       %画出最后一代的种群hold off
hold off
%% 画进化图
figure(2);
plot(1:MAXGEN,trace(2,:));
grid on 
xlabel('遗传代数')
ylabel('解的变化')
title('进化过程')
bestY = trace(2,end);
bestX = trace(1,end);
fprintf(['最优解：\nX = ',num2str(bestX),'\nY = ',num2str(bestY),'\n'])

% ga_2
clc
clear all
close all

%% 画出函数图
figure(1);
lbx = -2;ubx = 2;
lby = -2;uby = 2;
ezmesh('y*sin(2*pi*x)+x*cos(2*pi*y)',[lbx,ubx,lby,uby],50);    %画出函数曲线
hold on;

%% 定义遗传算法参数
NIND = 40;                  %种群大小
MAXGEN = 50;            %最大遗传代数
PRECI = 20;                 %个体长度
GGAP = 0.95;                %代沟
px = 0.7;                       %交叉变异
pm = 0.01;                      %变异概率
trace = zeros(3,MAXGEN);     %寻优结果的初始值
FieldD = [PRECI PRECI;lbx lby;ubx uby;1 1;0 0;1 1;1 1];         %区域扫描器
Chrom = crtbp(NIND,PRECI*2);                                      %创建任意离散随即种群

%% 优化
gen = 0;                            %代计数器
XY = bs2rv(Chrom,FieldD);                       %初代种群的十进制转换
X = XY(:,1);   Y = XY(:,2);
ObjV = Y.*sin(2*pi*X)+X.*cos(2*pi*Y);                   %计算目标函数值
while gen<MAXGEN
    FitnV = ranking(-ObjV);                                 %分配适应度值
    SelCh = select('sus',Chrom,FitnV,GGAP);             %选择
    SelCh = recombin('xovsp',SelCh,px);                         %重组
    SelCh = mut(SelCh,pm);                                      %变异
    XY = bs2rv(SelCh,FieldD);                                   %子代个体的十进制转换
    X = XY(:,1);   Y = XY(:,2);
    ObjVSel = Y.*sin(2*pi*X)+X.*cos(2*pi*Y);         %计算子代的目标函数
    [Chrom,ObjV] = reins(Chrom,SelCh,1,1,ObjV,ObjVSel);              %重插入子代到父代，得到新种群
    XY= bs2rv(Chrom,FieldD);
    gen = gen + 1 ;                               %代计数器增加
    %获取每代的最优解及其序号，Y为最优解，I为个体的序号
    [Y,I] = max(ObjV);
    trace(1:2,gen) = XY(I,:);                       %记下每代的最优值
    trace(3,gen) = Y;                                                              %记下每代的最优解
end
plot3(trace(1,:),trace(2,:),trace(3,:),'bo');          %画出每代的最优点
grid on;
plot3(XY(:,1),XY(:,2),ObjV,'bo');                               %画出最后一代的种群
hold off

%% 画进化图
figure(2);
plot(1:MAXGEN,trace(3,:));
grid on
xlabel('遗传代数')
ylabel('解的变化')
title('进化进程')
bestZ = trace(3,end);
bestX = trace(1,end);
bestY = trace(2,end);
fprintf(['最优解:\nX =',num2str(bestX),'\nY = ',num2str(bestY),'\nz = ',num2str(bestZ),'\n'])

% ga_tsp
clear;clc;
% 输入数据
vertexs=importdata('city.txt');          %城市坐标
n=length(vertexs);                        %城市数目
dist=zeros(n);                            %城市距离矩阵
for i = 1:n
    for j = 1:n
        dist(i,j)=distance(vertexs(i,:),vertexs(j,:));
    end
end
% 遗传算法参数设置
NIND=50;                        %种群大小
MAXGEN=150;                     %迭代次数
Pm=0.9;                         %交叉概率
Pc=0.1;                         %变异概率
pSwap=0.2;                      %选择交换结构的概率
pReversion=0.5;                 %选择逆转结构的概率
pInsertion=0.3;                 %选择插入结构的概率
N=n;                            %染色体长度=城市数目
% 种群初始化
Chrom=InitPop(NIND,N);
% 优化
gen=1;                          %计数器
bestChrom=Chrom(1,:);           %初始全局最优个体
bestL=RouteLength(bestChrom,dist);%初始全局最优个体的总距离
BestChrom=zeros(MAXGEN,N);      %记录每次迭代过程中全局最优个体
BestL=zeros(MAXGEN,1);          %记录每次迭代过程中全局最优个体的总距离
while gen<=MAXGEN
    SelCh=BinaryTourment_Select(Chrom,dist);            %二元锦标赛选择          
    SelCh=Recombin(SelCh,Pm);                           %OX交叉
    SelCh=Mutate(SelCh,Pc,pSwap,pReversion,pInsertion); %变异
    Chrom=SelCh;                                        %将Chrom更新为SelCh
    Obj=ObjFunction(Chrom,dist);                        %计算当前代所有个体总距离
    [minObj,minIndex]=min(Obj);                         %找出当前代中最优个体
    if minObj<=bestL                                    %将当前代中最优个体与全局最优个体进行比较，如果当前代最优个体更好，则将全局最优个体进行替换
        bestChrom=Chrom(minIndex,:); 
        bestL=minObj;
    end
    BestChrom(gen,:)=bestChrom;                         %记录每一代全局最优个体，及其总距离
    BestL(gen,:)=bestL;
    disp(['第' num2str(gen) '次迭代：全局最优路线总距离 = ' num2str(bestL)]); %显示外层循环每次迭代的信全局最优路线的总距离
    figure(1);                                          %画出每次迭代的全局最优路线图
    PlotRoute(bestChrom,vertexs(:,1),vertexs(:,2))
    pause(0.01);
    gen=gen+1;    %计数器加1
end
figure;                                                 % 打印每次迭代的全局最优个体的总距离变化趋势图
plot(BestL,'LineWidth',1);
title('优化过程')
xlabel('迭代次数');
ylabel('总距离');

function dist = distance(a,b)
%a          第一个城市坐标
%b          第二个城市坐标
%dist       两个城市之间距离
    x = (a(1)-b(1))^2;
    y = (a(2)-b(2))^2;
    dist = (x+y)^(1/2);
end

function Dist=ObjFunction(Chrom,dist)
%Chrom           种群
%dist            距离矩阵
%Dist            每个个体的目标函数值，即每个个体的总距离
    NIND=size(Chrom,1);                     %种群大小
    Dist=zeros(NIND,1);                      %目标函数初始化为0
    for i=1:NIND
        route=Chrom(i,:);                   %当前个体
        Dist(i,1)=RouteLength(route,dist);   %计算当前个体的总距离
    end
end

function L=RouteLength(route,dist)
%route               路线
%dist                距离矩阵
%L                   该路线总距离
    n=length(route);
    route=[route route(1)];
    L=0;
    for k=1:n 
        i=route(k);
        j=route(k+1); 
        L=L+dist(i,j); 
    end
end

function Chrom=InitPop(NIND,N)
%种群初始化
%NIND        种群大小
%N           染色体长度
%Chrom       随机生成的初始种群
    Chrom=zeros(NIND,N);                %种群初始化为NIND行N列的零矩阵
    for i=1:NIND
        Chrom(i,:)=randperm(N);         %每个个体为1~N的随机排列
    end
end

function FitnV=Fitness(Obj)
%适应度函数，总距离的倒数    
%输入Obj：     每个个体的总距离
%输出FitnV：   每个个体的适应度值，即总距离的倒数
    FitnV=1./Obj;
end

function Selch=BinaryTourment_Select(Chrom,dist)
%Chrom           种群
%dist            距离矩阵
%Selch           二元锦标赛选择出的个体
    Obj=ObjFunction(Chrom,dist);            %计算种群目标函数值，即每个个体的总距离
    FitnV=Fitness(Obj);                     %计算每个个体的适应度值，即总距离的倒数
    [NIND,N]=size(Chrom);                   %NIND-种群个数、N-种群长度
    Selch=zeros(NIND,N);                    %初始化二元锦标赛选择出的个体
    for i=1:NIND
        R=randperm(NIND);                   %生成一个1~NIND的随机排列
        index1=R(1);                        %第一个比较的个体序号
        index2=R(2);                        %第二个比较的个体序号
        fit1=FitnV(index1,:);               %第一个比较的个体的适应度值（适应度值越大，说明个体质量越高）
        fit2=FitnV(index2,:);               %第二个比较的个体的适应度值
        %如果个体1的适应度值 大于等于 个体2的适应度值，则将个体1作为第i选择出的个体
        if fit1>=fit2
            Selch(i,:)=Chrom(index1,:);
        else
            %如果个体1的适应度值 小于 个体2的适应度值，则将个体2作为第i选择出的个体
            Selch(i,:)=Chrom(index2,:);
        end
    end
end

function SelCh=Recombin(SelCh,Pc)
% 交叉操作
%SelCh    被选择的个体
%Pc       交叉概率
% SelCh   交叉后的个体
    NSel=size(SelCh,1);
    for i=1:2:NSel-mod(NSel,2)
        if Pc>=rand %交叉概率Pc
            [SelCh(i,:),SelCh(i+1,:)]=OX(SelCh(i,:),SelCh(i+1,:));
        end
    end
end
function [a,b]=OX(a,b)
%输入：a和b为两个待交叉的个体
%输出：a和b为交叉后得到的两个个体
    L=length(a);
    while 1
        r1=randsrc(1,1,[1:L]);
        r2=randsrc(1,1,[1:L]);
        if r1~=r2
            s=min([r1,r2]);
            e=max([r1,r2]);
            a0=[b(s:e),a];
            b0=[a(s:e),b];
            for i=1:length(a0)
                aindex=find(a0==a0(i));
                bindex=find(b0==b0(i));
                if length(aindex)>1
                    a0(aindex(2))=[];
                end
                if length(bindex)>1
                    b0(bindex(2))=[];
                end
                if i==length(a)
                    break
                end
            end
            a=a0;
            b=b0;
            break
        end
    end
end
 

function SelCh=Mutate(SelCh,Pm,pSwap,pReversion,pInsertion)
% 变异操作
%SelCh           被选择的个体
%Pm              变异概率
%pSwap           选择交换结构的概率
%pReversion      选择逆转结构的概率
%pInsertion      选择插入结构的概率
%SelCh           变异后的个体
    NSel=size(SelCh,1);
    for i=1:NSel
        if Pm>=rand
            index=Roulette(pSwap,pReversion,pInsertion);
            route1=SelCh(i,:);
            if index==1                %交换结构
                route2=Swap(route1);
            elseif index==2            %逆转结构
                route2=Reversion(route1);
            else                       %插入结构
                route2=Insertion(route1);
            end
            SelCh(i,:)=route2;
        end
    end
end

function index=Roulette(pSwap,pReversion,pInsertion)
%pSwap           选择交换结构的概率
%pReversion      选择逆转结构的概率
%pInsertion      选择插入结构的概率
%index           最终选择的邻域结构
    p=[pSwap pReversion pInsertion];
    r=rand;
    c=cumsum(p);
    index=find(r<=c,1,'first');
end

function route2=Swap(route1)
%交换操作
%route1          原路线1
%route2          经过交换结构变换后的路线2
    n=length(route1);
    seq=randperm(n);
    I=seq(1:2);
    i1=I(1);
    i2=I(2);
    route2=route1;
    route2([i1 i2])=route1([i2 i1]);
end

function route2=Reversion(route1)
%逆转变换
%route1          路线1
%route2          经过逆转结构变换后的路线2
    n=length(route1);
    seq=randperm(n);
    I=seq(1:2);
    i1=min(I);
    i2=max(I);
    route2=route1;
    route2(i1:i2)=route1(i2:-1:i1);
end

function route2=Insertion(route1)
%插入变换
%route1          路线1
%route2          经过插入结构变换后的路线2
    n=length(route1);
    seq=randperm(n);
    I=seq(1:2);
    i1=I(1);
    i2=I(2);
    if i1<i2
        route2=route1([1:i1-1 i1+1:i2 i1 i2+1:end]);
    else
        route2=route1([1:i2 i1 i2+1:i1-1 i1+1:end]);
    end
end

function PlotRoute(route,x,y)
%route           路线
%x,y             x,y坐标
    route=[route route(1)];
    plot(x(route),y(route),'k-o','MarkerSize',10,'MarkerFaceColor','w','LineWidth',1.5);
    xlabel('x');
    ylabel('y');
end