前言
这个例子展示了如何使用遗传算法来最小化使用自定义数据类型的函数。对遗传算法进行了定制化处理以解决旅行商问题。
一、旅行商问题（TSP）
旅行推销员问题（英语：Travelling salesman problem, TSP）是这样一个问题：给定一系列城市和每对城市之间的距离，求解访问每一座城市一次并回到起始城市的最短回路。它是组合优化中的一个NP难问题，在运筹学和理论计算机科学中非常重要。
旅行商问题(TravelingSalesmanProblem，TSP)是一个经典的组合优化问题。经典的TSP可以描述为：一个商品推销员要去若干个城市推销商品，该推销员从一个城市出发，需要经过所有城市后，回到出发地。应如何选择行进路线，以使总的行程最短。从图论的角度来看，该问题实质是在一个带权完全无向图中，找一个权值最小的Hamilton回路。由于该问题的可行解是所有顶点的全排列，随着顶点数的增加，会产生组合爆炸，它是一个NP完全问题。由于其在交通运输、电路板线路设计以及物流配送等领域内有着广泛的应用，国内外学者对其进行了大量的研究。早期的研究者使用精确算法求解该问题，常用的方法包括：分枝定界法、线性规划法、动态规划法等。但是，随着问题规模的增大，精确算法将变得无能为力，因此，在后来的研究中，国内外学者重点使用近似算法或启发式算法，主要有遗传算法、模拟退火法、蚁群算法、禁忌搜索算法、贪婪算法和神经网络等。
二、MATLAB 步骤
1.引入库
旅行商问题是一个有有限个城市且每个城市之间的旅行费用已知的优化问题。目标是找到一个有序的集合，其中包含所有的城市供销售人员访问，使成本最小化。为了解决旅行商问题，我们需要一个城市位置和每个城市之间的距离(或成本)的列表。
我们的推销员正在美国的各个城市访问。文件 usborder.mat 包含以变量 x 和 y 表示的美国地图，以及以变量 xx 和 yy 表示的相同地图的几何简化版本。
load('usborder.mat','x','y','xx','yy'); plot(x,y,'Color','red'); hold on;
我们将随机生成美国边境内城市的位置。我们可以使用 inpolygon 函数来确保所有城市都在美国边界之内或非常接近美国边界。
cities = 40; locations = zeros(cities,2); n = 1; while (n <= cities) xp = rand*1.5; yp = rand; if inpolygon(xp,yp,xx,yy) locations(n,1) = xp; locations(n,2) = yp; n = n+1; end end plot(locations(:,1),locations(:,2),'bo');
蓝色圆圈代表销售人员需要出差、送货或取货的城市位置。给定城市位置列表，我们可以计算所有城市的距离矩阵。
distances = zeros(cities); for count1=1:cities, for count2=1:count1, x1 = locations(count1,1); y1 = locations(count1,2); x2 = locations(count2,1); y2 = locations(count2,2); distances(count1,count2)=sqrt((x1-x2)^2+(y1-y2)^2); distances(count2,count1)=distances(count1,count2); end; end;
2. 为自定义数据类型定制遗传算法
默认情况下，遗传算法求解器解决基于 double 和二进制字符串数据类型的优化问题。用于创建、交叉和变异的函数假设总体是一个 double 类型的矩阵，如果是二进制字符串，则是合乎逻辑的。遗传算法求解器还可以处理涉及任意数据类型的优化问题。你可以使用任何你喜欢的数据结构。例如，可以使用 MATLAB® cell 数组指定自定义数据类型。为了将 ga 用于类型为 cell 数组的种群，您必须提供一个创建函数，交叉函数和变异函数，这些函数将对您的数据类型起作用，例如，单元数组。
3. 旅行商问题所需函数
本节将展示如何创建和注册这三个必需的函数。对于旅行商问题，种群中的个体是一个有序集合，因此可以很容易地用单元格数组表示种群。旅行商问题的自定义创建函数将创建一个单元格数组，例如 P，其中每个元素表示作为排列向量的城市的有序集合。也就是说，推销员将按照 P{i} 中指定的顺序行进。创建函数将返回一个大小为 PopulationSize 的单元格数组。
type create_permutations.m
function pop = create_permutations(NVARS,FitnessFcn,options) %CREATE_PERMUTATIONS Creates a population of permutations. % POP = CREATE_PERMUTATION(NVARS,FITNESSFCN,OPTIONS) creates a population % of permutations POP each with a length of NVARS. % % The arguments to the function are % NVARS: Number of variables % FITNESSFCN: Fitness function % OPTIONS: Options structure used by the GA % Copyright 2004-2007 The MathWorks, Inc. totalPopulationSize = sum(options.PopulationSize); n = NVARS; pop = cell(totalPopulationSize,1); for i = 1:totalPopulationSize pop{i} = randperm(n); end
自定义交叉函数接受一个单元格数组(population)，并返回一个单元格数组(交叉结果的子元素)。
type crossover_permutation.m
function xoverKids = crossover_permutation(parents,options,NVARS, ... FitnessFcn,thisScore,thisPopulation) % CROSSOVER_PERMUTATION Custom crossover function for traveling salesman. % XOVERKIDS = CROSSOVER_PERMUTATION(PARENTS,OPTIONS,NVARS, ... % FITNESSFCN,THISSCORE,THISPOPULATION) crossovers PARENTS to produce % the children XOVERKIDS. % % The arguments to the function are % PARENTS: Parents chosen by the selection function % OPTIONS: Options created from OPTIMOPTIONS % NVARS: Number of variables % FITNESSFCN: Fitness function % STATE: State structure used by the GA solver % THISSCORE: Vector of scores of the current population % THISPOPULATION: Matrix of individuals in the current population % Copyright 2004-2015 The MathWorks, Inc. nKids = length(parents)/2; xoverKids = cell(nKids,1); % Normally zeros(nKids,NVARS); index = 1; for i=1:nKids % here is where the special knowledge that the population is a cell % array is used. Normally, this would be thisPopulation(parents(index),:); parent = thisPopulation{parents(index)}; index = index + 2; % Flip a section of parent1. p1 = ceil((length(parent) -1) * rand); p2 = p1 + ceil((length(parent) - p1- 1) * rand); child = parent; child(p1:p2) = fliplr(child(p1:p2)); xoverKids{i} = child; % Normally, xoverKids(i,:); end
自定义的 mutation 函数接受一个个体，它是一个城市的有序集合，并返回一个修改后的有序集合。
type mutate_permutation.m
function mutationChildren = mutate_permutation(parents ,options,NVARS, ... FitnessFcn, state, thisScore,thisPopulation,mutationRate) % MUTATE_PERMUTATION Custom mutation function for traveling salesman. % MUTATIONCHILDREN = MUTATE_PERMUTATION(PARENTS,OPTIONS,NVARS, ... % FITNESSFCN,STATE,THISSCORE,THISPOPULATION,MUTATIONRATE) mutate the % PARENTS to produce mutated children MUTATIONCHILDREN. % % The arguments to the function are % PARENTS: Parents chosen by the selection function % OPTIONS: Options created from OPTIMOPTIONS % NVARS: Number of variables % FITNESSFCN: Fitness function % STATE: State structure used by the GA solver % THISSCORE: Vector of scores of the current population % THISPOPULATION: Matrix of individuals in the current population % MUTATIONRATE: Rate of mutation % Copyright 2004-2015 The MathWorks, Inc. % Here we swap two elements of the permutation mutationChildren = cell(length(parents),1);% Normally zeros(length(parents),NVARS); for i=1:length(parents) parent = thisPopulation{parents(i)}; % Normally thisPopulation(parents(i),:) p = ceil(length(parent) * rand(1,2)); child = parent; child(p(1)) = parent(p(2)); child(p(2)) = parent(p(1)); mutationChildren{i} = child; % Normally mutationChildren(i,:) end
我们还需要一个适合旅行商问题的适应度函数。个体的适应度是为一组有序城市所走的总距离。适应度函数也需要距离矩阵来计算总距离。
type traveling_salesman_fitness.m
function scores = traveling_salesman_fitness(x,distances) %TRAVELING_SALESMAN_FITNESS Custom fitness function for TSP. % SCORES = TRAVELING_SALESMAN_FITNESS(X,DISTANCES) Calculate the fitness % of an individual. The fitness is the total distance traveled for an % ordered set of cities in X. DISTANCE(A,B) is the distance from the city % A to the city B. % Copyright 2004-2007 The MathWorks, Inc. scores = zeros(size(x,1),1); for j = 1:size(x,1) % here is where the special knowledge that the population is a cell % array is used. Normally, this would be pop(j,:); p = x{j}; f = distances(p(end),p(1)); for i = 2:length(p) f = f + distances(p(i-1),p(i)); end scores(j) = f; end
ga 将只用一个参数 x 调用我们的适应度函数，但我们的适应度函数有两个参数：x，distances。我们可以使用一个匿名函数来获取额外的参数——距离矩阵的值。我们创建一个函数，将 FitnessFcn 处理为一个匿名函数，该函数接受一个输入 x，但使用 x 和 distances 调用 traveling_salesman_fitness。在创建函数 FitnessFcn 时，变量 distance 有一个值，因此这些值将被匿名函数捕获。
%distances defined earlier FitnessFcn = @(x) traveling_salesman_fitness(x,distances);
我们可以添加一个自定义绘图函数来绘制城市的位置和当前的最佳路线。红色圆圈表示城市，蓝色线表示两个城市之间的有效路径。
type traveling_salesman_plot.m
function state = traveling_salesman_plot(options,state,flag,locations) % TRAVELING_SALESMAN_PLOT Custom plot function for traveling salesman. % STATE = TRAVELING_SALESMAN_PLOT(OPTIONS,STATE,FLAG,LOCATIONS) Plot city % LOCATIONS and connecting route between them. This function is specific % to the traveling salesman problem. % Copyright 2004-2006 The MathWorks, Inc. persistent x y xx yy if strcmpi(flag,'init') load('usborder.mat','x','y','xx','yy'); end plot(x,y,'Color','red'); axis([-0.1 1.5 -0.2 1.2]); hold on; [unused,i] = min(state.Score); genotype = state.Population{i}; plot(locations(:,1),locations(:,2),'bo'); plot(locations(genotype,1),locations(genotype,2)); hold off
我们将再次使用一个匿名函数来创建一个匿名函数的函数句柄，该函数调用带有额外参数位置的 traveling_salesman_plot。
%locations defined earlier my_plot = @(options,state,flag) traveling_salesman_plot(options, ... state,flag,locations);
4.设置遗传算法选项
首先，我们将创建一个选项容器来指示自定义数据类型和人口范围。
options = optimoptions(@ga, 'PopulationType', 'custom','InitialPopulationRange', ... [1;cities]);
我们选择了已经创建的自定义创建、交叉、变异和绘图函数，并设置了一些停止条件。
options = optimoptions(options,'CreationFcn',@create_permutations, ... 'CrossoverFcn',@crossover_permutation, ... 'MutationFcn',@mutate_permutation, ... 'PlotFcn', my_plot, ... 'MaxGenerations',500,'PopulationSize',60, ... 'MaxStallGenerations',200,'UseVectorized',true);
最后，调用遗传算法求解问题。
numberOfVariables = cities; [x,fval,reason,output] = ... ga(FitnessFcn,numberOfVariables,[],[],[],[],[],[],[],options)
ga stopped because it exceeded options.MaxGenerations.
x = {[14 12 36 3 5 11 40 25 38 37 7 30 28 10 23 21 27 4 1 29 26 31 9 24 16 19 13 8 18 2 22 34 6 35 20 17 15 39 33 32]}
fval = 5.3846
reason = 0
output = problemtype: 'unconstrained' rngstate: [1×1 struct] generations: 500 funccount: 28563 message: 'ga stopped because it exceeded options.MaxGenerations.' maxconstraint: [] hybridflag: []
该图用蓝色圆圈显示了城市的位置，以及遗传算法找到的推销员将要走过的路径。推销员可以从路线的任何一端出发，也可以在终点回到出发城市再回到家。
