智能优化算法及其Matlab实现——遗传算法

本文参考：

• https://search.bilibili.com/all?vt=49410436&keyword=%E9%81%97%E4%BC%A0%E7%AE%97%E6%B3%95&from_source=webtop_search&spm_id_from=333.1007&search_source=3
• 《智能优化算法及其MATLAB实例》-包子阳

1遗传算法概念

• 遗传算法的定义
    遗传算法是一类模拟自然选择和遗传机制的搜索优化方法，通过模仿生物进化过程中的选择、交叉和变异等操作，来解决复杂的优化问题。
• 遗传算法的原理
    遗传算法基于达尔文的生物进化论，通过模拟自然选择的过程，不断迭代生成新的解，直到找到满足条件的最优解或近似最优解
• 遗传算法的应用
    遗传算法广泛应用于机器学习、人工智能、工程设计等领域，可以解决各种复杂的优化问题，如路径规划、参数优化等；遗传算法能有效地求解NP问题（Non-deterministic Polynomial problem）以及非线性、多峰函数优化和多目标优化问题。

2.遗传算法理论基础

2.1遗传算法理论——生物学基础

  自然选择——“适者生存，不适者淘汰”

• 遗传和变异是决定生物进化的内在因素。
• 遗传是指父代与子代之间，在性状上存在的相似现象；
  
• 变异是指父代与子代之间，以及子代的个体之间，在性状上存在的差异现象。

生物遗传和进化的规律有：
（1）生物的所有遗传信息都包含在其染色体中，染色体决定了生物的性状。染色体是由基因及其有规律的排列所构成的。
（2）生物的繁殖过程是由其基因的复制过程来完成的。同源染色体的交叉或变异会产生新的物种，使生物呈现新的性状。
（3）对环境适应能力强的基因或染色体，比适应能力差的基因或染色体有更多的机会遗传到下一代。

  总结：通过对环境的选择、基因的交叉和变异 这一生物演化的迭代过程的模仿，提出了能够用于求解最优化问题的强鲁棒性和自适应性的遗传算法。

2.2遗传算法——理论基础

2.2.1 模式定理

模式定理的核心思想

  模式定理指出：适应度高于种群平均、定义距短、低阶的模式在遗传算法中会以指数速度增长。其数学表达式为：
$$m(H,t+1)\ge m(H,t)\ast \frac{f(H)}{\bar{f}}[1-p_c\ast \delta (H)/L-1-o(H)\ast p_m]$$
其中：
$m(H,t)$：第t代中模式H的个体数。
$f(H)$：模式H的平均适应度。
$\bar{f}$ ：种群的平均适应度。
$p_c$：交叉概率，p_m：变异概率。
$\delta (H)$：模式的定义距。
$o(H)$：模式的阶。
$L$：染色体长度。

关键机制分析

• 选择操作：
• 适应度高的个体被选中的概率更高，因此适应度高于平均的模式H在下一代中期望数量增长，比例为$f(H)/\bar{f}$。

2. 交叉操作：

• 单点交叉可能破坏模式。破坏概率与定义距δ(H)成正比，近似为$p_c\ast \frac{\delta (H)}{L-1}$。
• 短定义距的模式更不易被破坏。例如，定义距为1的模式仅在交叉点位于其确定位之间时被破坏，概率极低。

3. 变异操作：

• 变异可能改变模式中的确定位。破坏概率近似为$o(H)\ast p_m$。
• 低阶模式（确定位少）更易存活。例如，阶为1的模式仅需一个确定位不变异即可保留。

积木块假设
  积木块定义：具有低阶、短定义距以及高平均适应度的模式称作积木块。(之所以称之为积木块，是由于遗传算法的求解过程并不是在搜索空间中逐一地测试各个基因的枚举组合，而是通过一些较好的模式，像搭积木一样，将它们拼接在一起,从而逐渐地构造出适应度越来越高的个体编码串。)

  积木块假设：个体的积木块通过选择、交叉、变异等遗传算子的作用，能够相互结合在一起，形成高阶、长距、高平均适应度的个体编码串 。(模式定理说明了积木块的样本数呈指数级增长，亦即说明了用遗传算法寻求最优样本的可能性，但它并未指明遗传算法一定能够寻求到最优样本；而积木块假设却说明了遗传算法的这种能力。)

  总结：从遗传算法的模式定理得到：具有高适应度、低阶、短定义矩的模式的数量会在种群的进化中呈指数级增长，从而保证了算法获得最优解的一个必要条件。而另一方面，积木块假设则指出：遗传算法有能力使优秀的模式向着更优的方向进化，即遗传算法有能力搜索到全局最优解。

2.2.2 基本概念

遗传编码

二进制编码 EX：求实数区间[0，4]上函数f （x ）的最大值?

• 可以由长度为6的位串表示变量 x ，即从“000000”到“111111”，并将中间的取值映射到实数区间[0，4]内。

• 由于从整数上来看，6位长度的二进制编码位串可以表示0～63，所以对应[0，4]的区间，每个相邻值之间的阶跃值为4/63≈0.0635，这个就是编码精度。

• 一般来说，编码精度越高，所得到的解的质量也越高，意味着解更为优良；但同时，由于遗传操作所需的计算量也更大，因此算法的耗时将更长。因而在解决实际问题时，编码位数需要适当选择。
  实数编码：

• 基于二进制编码的个体尽管操作方便，计算简单，但也存在着一些难以克服的困难而无法满足所有问题的要求。
    例如，对于高维、连续优化问题，由于从一个连续量离散化为一个二进制量本身存在误差，使得算法很难求得精确解。而要提高解的精度又必须加长编码串的长度，造成解空间扩大，算法效率下降。
    同时，二进制编码也不利于反映所求问题的特定知识，对问题信息和知识利用得不充分也会阻碍算法效率的进一步提高。为了解决二进制编码产生的问题，人们在解决一些数值优化问题（尤其是高维、连续优化问题）时，经常采用实数编码方式。
    实数编码的优点是计算精确度高，便于和经典连续优化算法结合，适用于数值优化问题；但其缺点是适用范围有限，只能用于连续变量问题。

二进制编码 一个长度为n的串，能表示多少个数呢？$2^n$。
如图：在区间[0,10]上，区间长度为L：10-1=9；二进制长度为$n：2$。对应精度为$E：\frac{10-1}{2^2-1}=3$。

三者关系为：$\frac{L}{2^n-1}=E$。

二进制解码 以[1,0]为例，其十进制转换结果为1∗2⁰ +1∗2¹ = 2
对应表示的数值为 7= 1+2*3

一般的，区间为[a,b]，区间长度为L，即L=b-a，串长为n，当前串对应十进制为T，则该串对应实值解为：
$$X=a+T\ast \frac{b-a}{2^n-1}$$

遗传操作

  遗传操作是优选强势个体的“选择”、个体间交换基因产生新个体 的“交叉”、个体基因信息突变而产生新个体的“变异”这三种变换的统 称。在生物进化过程中，一个群体中生物特性的保持是通过遗传来继 承的。生物的遗传主要是通过选择、交叉、变异三个过程把当前父代群体的遗传信息遗传到下一代（子代）成员。与此对应，遗传算法中 最优解的搜索过程也模仿生物的这个进化过程，使用所谓的遗传算子 来实现，即选择算子、交叉算子、变异算子。

• 选择算子：根据个体的适应度，按照一定的规则或方法，从第t 代群体P （t ）中选择出一些优良的个体遗传到下一代群体P （t+1）中。
    其中，“轮盘赌”*（轮盘赌参考代码理解更容易或者参考文章开头的链接视频）*选择法是遗传算法中最早提出的一种选择方法，由Holland提出，因为它简单实用，所以被广泛采用。它是一种基于比例的选择，利用各个个体适应度所占比例的大小来决定其子孙保留可能性。若某个个体i 的适应度为fi ，种群大小为NP ，则它被选取的概率表示为:
  $$p_i=\frac{f_i}{\sum _{i=1}^{NP}{f}_i}$$
    个体适应度越大，则其被选择的机会也越大；反之亦然。为了选择交叉个体，需要进行多轮选择。每一轮产生一个[0，1]内的均匀随机数，将该随机数作为选择指针来确定被选个体。
• 交叉算子：将群体P (t )中选中的各个个体随机搭配，对每一对个体，以某一概率（交叉概率P c ）交换它们之间的部分染色体。通过交叉，遗传算法的搜索能力得以飞跃提高。
    交叉操作一般分为以下几个步骤：首先，从交配池中随机取出交配的一对个体；然后，根据位串长度 L ，对要交配的一对个体，随机选取[1，L -1]中的一个或多个整数k 作为交叉位置；最后，根据交叉概率P c 实施交叉操作，配对个体在交叉位置处，相互交换各自的部分基因，从而形成新的一对个体。
  
• 变异算子：对群体中的每个个体，以某一概率（变异概率 Pm ）将某一个或某一些基因座上的基因值改变为其他的等位基因值。根据个体编码方式的不同，变异方式有：实值变异、二进制变异。对于二进制的变异，对相应的基因值取反；对于实值的变异，对相应的基因值用取值范围内的其他随机值替代。
    变异操作的一般步骤为：首先，对种群中所有个体按事先设定的变异概率判断是否进行变异；然后，对进行变异的个体随机选择变异位进行变异。
  

关键参数

• 群体规模：NP
    群体规模将影响遗传优化的最终结果以及遗传算法的执行效率。当群体规模NP 太小时，遗传优化性能一般不会太好。采用较大的群体规模可以减小遗传算法陷入局部最优解的机会，但较大的群体规模意味着计算复杂度较高。一般 NP 取10～200。
• 交叉概率：$p_c$
    交叉概率$p_c$ 控制着交叉操作被使用的频度。较大的交叉概率可以增强遗传算法开辟新的搜索区域的能力，但高性能的模式遭到破坏的可能性增大；若交叉概率太低，遗传算法搜索可能陷入迟钝状态。一般$p_c$ 取0.25～1.00。
• 变异概率： $p_m$
    变异在遗传算法中属于辅助性的搜索操作，它的主要目的是保持群体的多样性。一般低频度的变异可防止群体中重要基因的可能丢失，高频度的变异将使遗传算法趋于纯粹的随机搜索。通常$p_m$ 取0.001～0.1。

遗传运算的终止进化代数 G
终止进化代数G 是表示遗传算法运行结束条件的一个参数，它表示遗传算法运行到指定的进化代数之后就停止运行，并将当前群体中的最佳个体作为所求问题的最优解输出。一般视具体问题而定，G 的取值可在100～1000之间。

2.2.3 遗传算法流程

（1）初始化。设置进化代数计数器g =0，设置最大进化代数G ， 随机生成NP 个个体作为初始群体P （0）。
（2）个体评价。计算群体P （t ）中各个个体的适应度。
（3）选择运算。将选择算子作用于群体，根据个体的适应度，按照一定的规则或方法，选择一些优良个体遗传到下一代群体。
（4）交叉运算。将交叉算子作用于群体，对选中的成对个体，以某一概率交换它们之间的部分染色体，产生新的个体。
（5）变异运算。将变异算子作用于群体，对选中的个体，以某一 概率改变某一个或某一些基因值为其他的等位基因。群体P （t ）经过选择、交叉和变异运算之后得到下一代群体P （t +1）。计算其适应度值，并根据适应度值进行排序，准备进行下一次遗传操作。
（6）终止条件判断：若g ≤G ，则g = g +1，转到步骤（2）；若g ＞ G ，则此进化过程中所得到的具有最大适应度的个体作为最优解输出，终止计算。

2.3遗传算法——代码实现（以应用案例介绍）

2.3.1 EX1：

使用遗传算法求解函数f(x) = x + 10sin(5x) +7cos(4x) 的最大值，其中x的取值范围为[0,10]。
这是一个有多个局部极值的函数，其函数值如图所示：

clear; close all; clc;
%% 函数f(x)=x+10sin(5x)+7cos(4x)
x = 0:0.01:10; 
y = x + 10*sin(5*x) + 7*cos(4*x); 
plot(x,y);
xlabel('x');
ylabel('f(x)');
title('f(x)=x+10sin(5x)+7cos(4x)');

代码实现：
（1）初始化种群数目为NP =50，染色体二进制编码长度为L=20，最大进化代数为G =100，交叉概率为Pc =0.8，变异概率为Pm=0.1。

%% 遗传算法求解
clear; % 清除所有变量
close all; %清图
clc; % 清屏
NP = 50; % 种群数量
L = 20; % 二进制位长度
Pc = 0.8; % 交叉概率
Pm = 0.1; % 变异概率
G =100; % 最大遗传代数
Xs = 10; % 上限
Xx = 0; % 下限
xlim = [0 10] ; %x取值范围
pop = round(rand(NP,L));  % 随机获得初始种群(即NP个L长度的二进制串，50个染色体，每个染色体有20个基因)

（2）产生初始种群，将二进制编码转换成十进制，计算个体适应度值，并进行归一化；采用基于轮盘赌的选择操作、基于概率的交叉和变异操作，产生新的种群，并把历代的最优个体保留在新种群中，进行下一步遗传操作。第一次迭代如下：

decpop = bintodec(pop, NP, L, xlim); % 计算初代解对应十进制
fx = calobjvalue(decpop); %计算初代解的函数值
[y(1),l] = max(fx); x(1) = decpop(l); %第一次迭代最大值和，对应位置
plotfig(decpop, fx, xlim, 1) ; %绘制图形

（3）第2次到第G次迭代。判断是否满足终止条件：若满足，则结束搜索过程，输出优化值；若不满足，则继续进行迭代优化。

%% 2-G次迭代
%% 2-G次迭代
for i =2:G
    decpop = bintodec(pop, NP, L, xlim); %计算上一代解对应十进制
    fx = calobjvalue(decpop); % 计算上一代解对应的十进制 
    fitvalue = fx; %适应度映射
    newpop = copyx(pop, fitvalue, NP); %复制
    newpop = crossover(newpop , Pc, NP, L); %交叉
    newpop = mutation(newpop , Pm, NP, L); %变异
    
    % 对新一代群体进行择优保留（即实现保底机制）
    newdecpopo = bintodec(newpop , NP, L, xlim); %计算这一代解对应十进制
    new_fx = calobjvalue(newdecpopo ); % 计算这一代解对应的十进制
    new_fitvalue = new_fx; %适应度映射
    index = find(new_fitvalue > fitvalue);
    
    pop(index, : ) = newpop(index, : ); % 更新得到最新解
    decpop = bintodec(pop, NP, L, xlim); %计算十进制
    fx = calobjvalue(decpop); % 计算上一代解对应的十进制
    plotfig(decpop, fx, xlim, i) % 绘制新解的图

    [bestindividual,bestindex] = max(fx);% 找出更新后的个体最优函数
    y(i) = bestindividual ; % 记录每一代的最优函数值
    x(i) = decpop(bestindex); %十进制解
    subplot(1,2,2); plot(1:i,y); title('适应度进化曲线');
end
[ymax, max_index] = max(y); disp(['找的最优解位置为: ', num2str(x(max_index)) ]); disp(['对应最优解为: ', num2str(ymax)]);

• 二进制转十进制函数。

一般的，区间为[a,b]，区间长度为L，即L=b-a，串长为n，当前串对应十进制为T，则该串对应实值解为：
$$X=a+T\ast \frac{b-a}{2^n-1}$$

%% 二进制解码操作
function dec = bintodec(pop, NP, L, xlim)
    dec = zeros(1,NP);
    index = L-1 : -1 : 0;
    for i = 1: NP
        dec(i) = sum(pop(i,:).*(2.^index));
    end
    dec = xlim(1) +dec*( xlim(2)-xlim(1) ) / (2^L-1);
end

• 适应度函数（目标函数）
• f(x) = x + 10sin(5x) +7cos(4x)

%% 适应度函数
function fx= calobjvalue(decpop) % 参数为十进制
    f = @(x)x + 10*sin(5*x) + 7*cos(4*x);
    fx = f(decpop);
end

• 图形绘制函数

%% 绘图函数
function plotfig(decpop, fx, xlim, k)
    f = @(x)x + 10*sin(5*x) + 7*cos(4*x);
    x = xlim(1) : 0.05 : xlim(2);
    y = f(x);
    subplot(1,2,1);
    plot(x,y,decpop,fx,'o');
    title(['第',num2str(k),'次迭代进化']);
    pause(0.2);
end

• 轮盘赌复制操作函数(此处t图片只是轮盘赌的一个举例和被题无关)

%% 复制操作
function newx = copyx(pop, fitvalue, NP) % 输入二进制串和对应的适应度
	newx = pop; %开辟一个空间
	i=1;j=1;
	% 函数值域包含正负解，所以使用轮盘赌是，需要对适应度进行归一化处理
	maxFit = max(fitvalue); 
	minFit = min(fitvalue);
	Fit=(fitvalue - minFit)/(maxFit - minFit);
	p = Fit./ sum(Fit);
	Cs = cumsum(p);
	R = sort(rand(NP, 1)); % 每个个体的复制概率
	while j <= NP
		if R(j) < Cs(i)
			newx(j, :) = pop(i, :);
			j = j + 1;
		else
			i = i + 1;
			if i>NP
				return;
			end
		end
	end
end

• 交叉操作函数
  

%% 交叉操作
function newx = crossover(pop, Pc, NP, L) % 输入二进制串、交叉概率、种群数量、串长度
	newx = pop;% 开辟一个空间
	i = 2;
 	while i + 2 <= NP
		% 将第i个染色体和第i-1个染色体交叉
		if rand < Pc
			x1 = pop(i - 1, :);
			x2 = pop(i, :);
			r = randperm(L, 2);% 返回范围内两个整数作为交叉点
			r1 = min(r); r2 = max(r);% 交叉复制的位点
			newx(i-1, :) = [x1(1 : r1-1), x2(r1:r2), x1(r2+1 : end)];
			newx(i, :) = [x2(1 : r1-1), x1(r1:r2), x2(r2+1 : end)];
		end
		i = i + 2;
	end
end

• 变异操作函数
  
  本案例中的选择是用子代、父代两代中的优秀染色体；

%% 变异操作
function newx=mutation(pop, Pm, NP, L) % 输入二进制串、交叉概率、种群数量、串长度
	i = 1;
	while i + 2 <= NP
		if rand < Pm
			r = randperm(L, 1);
			pop(i, r) = ~pop(i, r);
		end
		i = i + 1;
		newx = pop;
	end
end

EX1 代码合并

clear; close all; clc;
%% 函数f(x)=x+10sin(5x)+7cos(4x)
x = 0:0.01:10; 
y = x + 10*sin(5*x) + 7*cos(4*x); 
plot(x,y);
xlabel('x');
ylabel('f(x)');
title('f(x)=x+10sin(5x)+7cos(4x)');
%% 遗传算法求解
clear; % 清除所有变量
close all; %清图
clc; % 清屏
NP = 50; % 种群数量
L = 20; % 二进制位长度
Pc = 0.8; % 交叉概率
Pm = 0.1; % 变异概率
G =100; % 最大遗传代数
Xs = 10; % 上限
Xx = 0; % 下限
xlim = [0 10] ; %x取值范围
pop = round(rand(NP,L));  % 随机获得初始种群(即NP个L长度的二进制串，50个染色体，每个染色体有20个基因)
decpop = bintodec(pop, NP, L, xlim); % 计算初代解对应十进制
fx = calobjvalue(decpop); %计算初代解的函数值
[y(1),l] = max(fx); x(1) = decpop(l); %第一次迭代最大值和，对应位置
plotfig(decpop, fx, xlim, 1) ; %绘制图形
%% 2-G次迭代
for i =2:G
    decpop = bintodec(pop, NP, L, xlim); %计算上一代解对应十进制
    fx = calobjvalue(decpop); % 计算上一代解对应的十进制 
    fitvalue = fx; %适应度映射
    newpop = copyx(pop, fitvalue, NP); %复制
    newpop = crossover(newpop , Pc, NP, L); %交叉
    newpop = mutation(newpop , Pm, NP, L); %变异
    
    % 对新一代群体进行择优保留（即实现保底机制）
    newdecpopo = bintodec(newpop , NP, L, xlim); %计算这一代解对应十进制
    new_fx = calobjvalue(newdecpopo ); % 计算这一代解对应的十进制
    new_fitvalue = new_fx; %适应度映射
    index = find(new_fitvalue > fitvalue);
    
    pop(index, : ) = newpop(index, : ); % 更新得到最新解
    decpop = bintodec(pop, NP, L, xlim); %计算十进制
    fx = calobjvalue(decpop); % 计算上一代解对应的十进制
    plotfig(decpop, fx, xlim, i) % 绘制新解的图

    [bestindividual,bestindex] = max(fx);% 找出更新后的个体最优函数
    y(i) = bestindividual ; % 记录每一代的最优函数值
    x(i) = decpop(bestindex); %十进制解
    subplot(1,2,2); plot(1:i,y); title('适应度进化曲线');
end
[ymax, max_index] = max(y); disp(['找的最优解位置为: ', num2str(x(max_index)) ]); disp(['对应最优解为: ', num2str(ymax)]);
%% 二进制解码操作
function dec = bintodec(pop, NP, L, xlim)
    dec = zeros(1,NP);
    index = L-1 : -1 : 0;
    for i = 1: NP
        dec(i) = sum(pop(i,:).*(2.^index));
    end
    dec = xlim(1) +dec*( xlim(2)-xlim(1) ) / (2^L-1);
end
%% 适应度函数
function fx= calobjvalue(decpop) % 参数为十进制
    f = @(x)x + 10*sin(5*x) + 7*cos(4*x);
    fx = f(decpop);
end
%% 绘图函数
function plotfig(decpop, fx, xlim, k)
    f = @(x)x + 10*sin(5*x) + 7*cos(4*x);
    x = xlim(1) : 0.05 : xlim(2);
    y = f(x);
    subplot(1,2,1);
    plot(x,y,decpop,fx,'o');
    title(['第',num2str(k),'次迭代进化']);
    pause(0.2);
end
%% 复制操作
function newx = copyx(pop, fitvalue, NP) % 输入二进制串和对应的适应度
	newx = pop; %开辟一个空间
	i=1;j=1;
	% 函数值域包含正负解，所以使用轮盘赌是，需要对适应度进行归一化处理
	maxFit = max(fitvalue); 
	minFit = min(fitvalue);
	Fit=(fitvalue - minFit)/(maxFit - minFit);
	p = Fit./ sum(Fit);
	Cs = cumsum(p);
	R = sort(rand(NP, 1)); % 每个个体的复制概率
	while j <= NP
		if R(j) < Cs(i)
			newx(j, :) = pop(i, :);
			j = j + 1;
		else
			i = i + 1;
			if i>NP
				return;
			end
		end
	end
end
%% 交叉操作
function newx = crossover(pop, Pc, NP, L) % 输入二进制串、交叉概率、种群数量、串长度
	newx = pop;% 开辟一个空间
	i = 2;
 	while i + 2 <= NP
		% 将第i个染色体和第i-1个染色体交叉
		if rand < Pc
			x1 = pop(i - 1, :);
			x2 = pop(i, :);
			r = randperm(L, 2);% 返回范围内两个整数作为交叉点
			r1 = min(r); r2 = max(r);% 交叉复制的位点
			newx(i-1, :) = [x1(1 : r1-1), x2(r1:r2), x1(r2+1 : end)];
			newx(i, :) = [x2(1 : r1-1), x1(r1:r2), x2(r2+1 : end)];
		end
		i = i + 2;
	end
end
%% 变异操作
function newx=mutation(pop, Pm, NP, L) % 输入二进制串、交叉概率、种群数量、串长度
	i = 1;
	while i + 2 <= NP
		if rand < Pm
			r = randperm(L, 1);
			pop(i, r) = ~pop(i, r);
		end
		i = i + 1;
		newx = pop;
	end
end

2.3.2 EX2：

计算函数f(x)=∑_i=1^n▒x_i2(−20≤x_i≤20)的最小值，其中个体x 的维数n =10。这是一个简单的平方和函数，只有一个极小点x=（0，0，…，0），理论最小值f （0，0，…，0）=0。
解： 仿真过程如下：
（1）初始化种群数目为NP =100，染色体基因维数为D =10，最大进化代数为G =1000，交叉概率为P c =0.8，变异概率为P m =0.1。

（2）产生初始种群，计算个体适应度值；进行实数编码的选择以及交叉和变异操作。选择和交叉操作采用“君主方案”，即在对群体根据适应度值高低进行排序的基础上，用最优个体与其他偶数位的所有个体进行交叉，每次交叉产生两个新的个体。在交叉过后，对新产生的群体进行多点变异产生子群体，再计算其适应度值，然后和父群体合并，并且根据适应度值进行排序，取前NP 个个体为新群体，进行下一次遗传操作。

（3）判断是否满足终止条件：若满足，则结束搜索过程，输出优化值；若不满足，则继续进行迭代优化。

clear all; clc; close all;
D = 10; % 基因数目
NP = 100;% 染色体数目
Xs = 20;% 上限
Xx = -20;% 下限
G = 1000;% 最大遗传迭代次数
f = zeros(D, NP); % 初始种群赋空间
nf = zeros(D, NP);% 子种群赋空间
Pc = 0.8; % 交叉概率
Pm = 0.1; % 变异概率
f = rand(D, NP) * (Xs - Xx) + Xx;% 随机获得初始种群

%% 按适应度升序排列
for np = 1:NP
    MSLL(np) = func2(f(:, np));
end
[SortMSLL, Index] = sort(MSLL);
Sortf = f(:, Index);
%% 遗传算法循环
for gen = 1:G
    % %采用君主方案进行选择交叉操作
    Emper = Sortf(:, 1);% 君主染色体
    NoPoint = round(D * Pc); % 每次交叉点的个数
    PoPoint = randi([1 D], NoPoint, NP / 2 );% 交叉基因的位置
    nf = Sortf;
    for i = 1:NP / 2
        nf(:, 2 * i - 1) = Emper;
        nf(:, 2 * i) = Sortf(:, 2 * i);
        for k = 1:NoPoint
            nf(PoPoint(k, i), 2 * i - 1) = nf(PoPoint(k, i), 2 * i);
            nf(PoPoint(k, i), 2 * i) = Emper(PoPoint(k, i));
        end
    end

    %% 变异操作
    for m = 1:NP
         for n = 1 : D
             r = rand(1, 1);
             if r<Pm
                 nf(n, m) = rand(1, 1) * (Xs - Xx) + Xx;
             end
         end
    end

    %% 子种群按适应度升序排列
    for np = 1:NP
        NMSLL(np) = func2(nf(:, np));
    end
    [NSortMSLL, Index] = sort(NMSLL);
    NSortf = nf(:, Index);

    %% 产生新种群
    f1 = [Sortf, NSortf]; % 子代和父代合并
    MSLL1 = [SortMSLL, NSortMSLL];% 子代和父代的适应度值合并
    [SortMSLL1, Index] = sort(MSLL1);% 适应度按升序排列
    Sortf1 = f1(:,Index); % 按适应度排列个体
    SortMSLL = SortMSLL1(1:NP); % 取前NP个适应度值
    Sortf = Sortf1(:, 1:NP);% 取前NP个个体
    trace(gen) = SortMSLL(1);% 历代最优适应度值
end

Bestf = Sortf(:, 1); % 最优个体值
trace(end);% 最优值
figure
plot(trace);
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度计划曲线')

%% 适应度函数
function result = func2(x)
    summ = sum(x.^2);
    result = summ;
end

2.3.3 EX3：

旅行商问题（TSP问题）。假设有一个旅行商人要拜访全国31个省会城市，他需要选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。对路径选择的要求是：所选路径的路程为所有路径之中的最小值。全国31个省会城市的坐标为 [1304 2312；3639 1315；4177 2244；3712 1399；3488 1535；3326 1556；3238 1229；4196 1004；4312790；4386 570；3007 1970；2562 1756；2788 1491；2381 1676；1332695；3715 1678；3918 2179；4061 2370；3780 2212；3676 2578；40292838；4263 2931；3429 1908；3507 2367；3394 2643；3439 3201；2935 3240；3140 3550；2545 2357；2778 2826；2370 2975]。

解： 仿真过程如下：
（1）初始化种群数目为NP =200，染色体基因维数为N =31，最大进化代数为G =1000。

（2）产生初始种群，计算个体适应度值，即路径长度；采用基于概率的方式选择进行操作的个体；对选中的成对个体，随机交叉所选中的成对城市坐标，以确保交叉后路径每个城市只到访一次；对选中的单个个体，随机交换其一对城市坐标作为变异操作，产生新的种群，进行下一次遗传操作。

（3）判断是否满足终止条件：若满足，则结束搜索过程，输出优化值；若不满足，则继续进行迭代优化。

%% 遗传算法解决TSP问题
% 清除所有变量
clear all;close all;clc;
C=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
3238 1229;4196 1044;4312 790;4386 570;3007 1970;2562 1756;...
2788 1491;2381 1676;1332 695;3715 1678;3918 2179;4061 2370;...
3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;...
3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;...
2370 2975];
%31个省会城市坐标TSP问题的规模，即城市数目
N=size(C,1);
%任意两个城市距离间隔矩阵
D=zeros(N);
%%求任意两个城市距离间隔矩阵
for i=1:N
	for j=1:N
		D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;
	end
end
NP=200; % 种群规模
G=1000; %最大遗传代数
f=zeros(NP,N); %用于存储种群
F=[]; %种群更新中间存储
for i=1:NP
    f(i,:)=randperm(N); %随机生成初始种群
end
R=f(1,:); % 存储最优种群
len =zeros(NP,1); %存储路径长度
fitness =zeros(NP,1); %存储归一化适应值
gen =0;
gen =0;
%% 遗传算法循环
while gen <G
	%% 计算路径长度
	for i=1:NP
		len(i,1)=D(f(i,N),f(i,1));
		for j=1:(N-1)
			len(i,1)=len(i,1)+D(f(i,j),f(i,j+1));
		end
	end
	% 最长路径
	maxlen = max(len);
	% 最短路径
	minlen = min(len);
	%更新最短路径
	rr = find(len==minlen);
	R=f(rr(1,1),:);
	%计算归一化适应值
	for i= 1:length(len)
		fitness(i,1)=(1-((len(i,1)-minlen)/(maxlen-minlen+0.001)));
	end
	%选择操作
	nn=0;
	for i=1:NP
		if fitness(i,1)>= rand
			nn= nn+1;
			F(nn,:)=f(i,:);
		end
	end
	[aa,bb]= size(F);
	while aa< NP
		nnper = randperm(nn);
		A=F(nnper(1),:);
		B=F(nnper(2),:);
		%交叉操作
		W=ceil(N/10); % 交叉点个数
		p= unidrnd(N-W+1);% 随机选择交叉范围,从p到p+w
		for i= 1:W
			x=find(A==B(p+i-1));
			y=find(B==A(p+i-1));
			temp=A(p+i-1);
			A(p+i-1)=B(p+i-1);
			B(p+i-1)=temp;
			temp=A(x);
			A(x)=B(y);
			B(y)=temp;
		end
		%变异操作
		p1=floor(1+N*rand());
		p2=floor(1+N*rand());
		while p1==p2
			p1= floor(1+N*rand());
			p2=floor(1+N*rand());
		end
		tmp=A(p1);
		A(p1)=A(p2);
		A(p2)=tmp;
		tmp=B(p1);
		B(p1)=B(p2);
		B(p2)=tmp;
		F=[F;A;B];
		[aa,bb]= size(F);
	end
	if aa >NP
		F=F(1:NP,:);% 保持种群规模为NP
	end
	f=F; % 更新种群
	f(1,:)=R; % 保留每代最优个体
	clear F;
	gen = gen+1;
	Rlength(gen)= minlen;
end
figure
for i=1:N-1
	plot([C(R(i),1),C(R(i+1),1)],[C(R(i),2),C(R(i+1),2)],'bo-');
	hold on;
end
plot([C(R(N),1),C(R(1),1)],[C(R(N),2),C(R(1),2)],'ro-');
title(['优化最短距离:',num2str(minlen)]);
figure
plot(Rlength)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')

2.3.4 EX4：

旅背包问题。有N 件物品和一个容量为V的背包。第i 件物品的体积是 c （i ），价值是 w （i ）。求解将哪些物品放入背包可使物品的体积总和不超过背包的容量，且价值总和最大。假设物品数量为10，背包的容量为300。每件物品的体积为[95，75，23，73，50，22，6，57，89，98]，价值为[89，59，19，43，100，72，44，16，7，64]。

解 ：仿真过程如下：
（1）初始化种群数目为NP =50，染色体基因维数为L =10，最大进化代数为G =100。

（2）产生二进制初始种群，其中1表示选择该物品，0表示不选择该物品。取适应度值为选择物品的价值总和，计算个体适应度值，当物品体积总和大于背包容量时，对适应度值进行惩罚计算。

（3）对适应度进行归一化，采用基于轮盘赌的选择操作、基于概率的交叉和变异操作，产生新的种群，并把历代的最优个体保留在新种群中，进行下一步遗传操作。

（4）判断是否满足终止条件：若满足，则结束搜索过程，输出优化值；若不满足，则继续进行迭代优化。
（1）先确定染色体、基因。
   a.背包问题的基因是确定的即每个基因对应每个物品，这里有10个物品，则基因数（基因串长度）L= 10；
   b.染色体数目（种群数目）可自己定义，这里取NP = 50；

（2）交叉概率Pc = 0.8、变异概率Pm = 0.05；遗代代数G = 100；

（3）背包容量 V = 300 ；物品体积C = [95，75，23，73，50，22，6，57，89，98]，物品价值为W=[89，59，19，43，100，72，44，16，7，64]；

（4）惩罚函数系数 afa = 2 ；
说明：随机从十个物品中选择物品放入背包，存在总容量超出限定范围的情况。此时采取评价策略，即总价值-超出容量部分*afa；这样便可以成功选择合适的染色体

clear all;close all;clc;
% 种群规模
NP= 50;
% 物品件数
L=10;
% 交叉率
Pc=0.8;
% 变异率
Pm=0.05;
% 最大遗传代数
G=100;
% 背包容量
V=300;
% 物品体积
C=[95,75,23,73,50,22,6,57,89,98];
% 物品价值
W=[89,59,19,43,100,72,44,16,7,64];
% 惩罚函数系数
afa= 2;

（5）由题目分析来看，物品只有取用和被取用两种情况。用基因来表示就是：0、1的串表示。
说明：f为50*10的数组表示，50个染色体；

% 随机获得初始种群
f = randi([0,1],NP,L);

（6）适应度函数。-matlab调用函数放置到末尾
$$y=\{\begin{array}{llc}\sum (容量\ast 价值),& \sum (容量）<背包容量& \\ \sum (容量\ast 价值)-[惩罚系数\ast （背包容量-\sum (容量））],& \sum (容量）<背包容量,& x>0\end{array}$$

%% 适应度函数
function result = func4(f,c,w,v,afa)
	fit = sum(f.*w);
	TotalSize = sum(f.*c);
	if TotalSize <= v
		fit = fit;
	else
		fit = fit - afa *(TotalSize -v);
	end
	result = fit;
end

（7）基于轮盘赌复制。
a.第一步求取每个染色体对应的适应度函数；即每个染色体对应的价值总量。
b.对每个概率进行归一化操作，即将其映射到[minFit,maxFit]上，并记录每代最优个体；
c.求取每个适应度对应的概率将其映射到[0,1]的轮盘上；
d.使用轮盘赌方法进行复制；

%% 遗传算法循环
for k=1:G
% 适应度计算
	for i=1:NP
		Fit(i)= func4(f(i,:),C,W,V,afa);
	end
	% 最大值
	maxFit = max(Fit);
	% 最小值
	minFit = min(Fit);
	rr = find(Fit==maxFit);
	% 历代最优个体
	fBest= f(rr(1,1),:);
	% 归一化适应度值
	Fit =(Fit - minFit)/(maxFit - minFit);
	%% 基于轮盘赌的复制操作 
	sum_Fit = sum(Fit);
	fitvalue = Fit./sum_Fit;
	fitvalue = cumsum(fitvalue);
	ms = sort(rand(NP,1));
	fiti= 1;
	newi = 1;
	while newi <= NP
		if(ms(newi))<fitvalue(fiti)
			nf(newi,:)=f(fiti,:);
			newi = newi + 1;
		else
			fiti= fiti+ 1;
		end
	end

（8）基于概率的交叉
注：随机生成的概率小于交叉概率，并进行交叉操作

%% 基于概率的交叉操作
	for i= 1:2:NP
		p = rand;
		if p<Pc
			q= randi([0,1],1,L);
			for j=1:L
				if q(j)==1;
					temp = nf(i+1,j);
					nf(i+1,j)=nf(i,j);
					nf(i,j)=temp;
				end
			end
		end
	end

（8）基于概率的变异
注：随机生成的概率小于变异概率，并进行变异操作，基因取反

	%% 基于概率的变异操作
	for m= 1:NP
		for n=1:L
			r = rand(1,1);
			if r< Pm
				nf(m,n)=~nf(m,n);
			end
		end
	end

（9）保留最优个体

	f = nf;
	%% 保留最优个体在新种群中
	f(1,:)=fBest;
	% 历代最优适应度
	trace(k)= maxFit;
end
% 最优个体
fBest;
figure
plot(trace)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')

EX4整合代码

clear all;close all;clc;
% 种群规模
NP= 50;
% 物品件数
L=10;
% 交叉率
Pc=0.8;
% 变异率
Pm=0.05;
% 最大遗传代数
G=100;
% 背包容量
V=300;
% 物品体积
C=[95,75,23,73,50,22,6,57,89,98];
% 物品价值
W=[89,59,19,43,100,72,44,16,7,64];
% 惩罚函数系数
afa= 2;
% 随机获得初始种群
f = randi([0,1],NP,L);
%% 遗传算法循环
for k=1:G
% 适应度计算
	for i=1:NP
		Fit(i)= func4(f(i,:),C,W,V,afa);
	end
	% 最大值
	maxFit = max(Fit);
	% 最小值
	minFit = min(Fit);
	rr = find(Fit==maxFit);
	% 历代最优个体
	fBest= f(rr(1,1),:);
	% 归一化适应度值
	Fit =(Fit - minFit)/(maxFit - minFit);
	%% 基于轮盘赌的复制操作 
	sum_Fit = sum(Fit);
	fitvalue = Fit./sum_Fit;
	fitvalue = cumsum(fitvalue);
	ms = sort(rand(NP,1));
	fiti= 1;
	newi = 1;
	while newi <= NP
		if(ms(newi))<fitvalue(fiti)
			nf(newi,:)=f(fiti,:);
			newi = newi + 1;
		else
			fiti= fiti+ 1;
		end
	end
	%% 基于概率的交叉操作
	for i= 1:2:NP
		p = rand;
		if p<Pc
			q= randi([0,1],1,L);
			for j=1:L
				if q(j)==1;
					temp = nf(i+1,j);
					nf(i+1,j)=nf(i,j);
					nf(i,j)=temp;
				end
			end
		end
	end
	%% 基于概率的变异操作
	for m= 1:NP
		for n=1:L
			r = rand(1,1);
			if r< Pm
				nf(m,n)=~nf(m,n);
			end
		end
	end
	f = nf;
	%% 保留最优个体在新种群中
	f(1,:)=fBest;
	% 历代最优适应度
	trace(k)= maxFit;
end
% 最优个体
fBest;
figure
plot(trace)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
%% 适应度函数
function result = func4(f,c,w,v,afa)
	fit = sum(f.*w);
	TotalSize = sum(f.*c);
	if TotalSize <= v
		fit = fit;
	else
		fit = fit - afa *(TotalSize -v);
	end
	result = fit;
end