数模（二）遗传算法（TSP、背包问题，含代码）

本文链接：https://blog.csdn.net/qq_23096319/article/details/125725448

遗传算法（简称 GA）是启发式算法的一种，模拟了“物竞天择，适者生存”的自然规律。所有的启发式算法，本质上是具有一定随机性的搜索算法，它的搜索过程不像暴力搜索那样费时，也不像贪心法一样容易陷入局部最优。具有随机性质的现象在生活、大自然界中有许多，命名为启发式，可能说的就是人在生活中看到了这些随机现象，受到启发，一拍脑袋想到了这些算法吧。

名词解释

生物遗传概念	遗传算法中的作用
适者生存	最优目标值的可行解有最大的可能被留住（算法思想，由选择算子实现）
个体	可行解
染色体	可行解的编码（很重要，因为 GA 要对可行解作解构式的操作，可以理解为编码是可行解的表征）
基因	可行解中每一分量的特征
适应性	适应度函数值（目标函数）
种群	根据适应度函数选取的一组可行解（生存下来的适者）
交配	通过交配原则产生一组新可行解的过程（由交叉算子实现）
变异	编码的某一分量发生变化的过程（由变异算子实现）

算法流程

在这里插入图片描述

（1）算法参数初始化，包括种群大小、进化代数以及变异率；
（2）对解进行编码、表征，并生成初代种群；
（3）设父代种群为 J ，首先繁衍（复制）得到种群 A，对 A 进行交叉操作后得到子代种群 A（大自然中交叉变异是在繁衍过程中发生的，在这里我们遵循自然规律，同时注意保留父代）；
（4）对 A 进行变异操作得到子代种群 B；
（5）此时我们获得父代和子代混合种群 $G = [J; A; B]$ ，计算每个个体的适应度，下一步进行“优胜劣汰”的选择操作，维持种群大小；
（6）若达到最大进化代数或者最优解满足一定条件，则转到（7），否则返回（3）；
（7）解编码，输出结果，算法结束.

编码、解码

实数编码

用实数表示基因，做函数优化时不需要解码过程，但容易过早收敛，从而陷入局部最优。
下面举两个实数编码的例子：

函数优化
平面上给定一棵树，每个节点都有一个坐标，我们可以移动分枝点（除根节点和叶子节点之外的节点），求树所有边的最小和。
模型建立
首先，设树的邻接矩阵为 G，叶子节点和根节点数量为 $n_1$ ，依次编号为 $1-n_1$ ，分枝点数量为 $n_2$ ，编号为 $n_1+1)-(n_1+n_2)$ ，它们的位置为 $X_i,Y_i),i=1,2,...,(n_1+n_2)$ .
设编号依次为 $n_1+1)-(n_1+n_2)$ 的节点移动量为 $x_{2i-1},x_{2i})$ ，我们得到一个非线性模型：
$min\sum_{e(i,j)\in{G}}{dis(i,j)}$
$\begin{cases} X_{i+n_1}=X_{i+n_1}+x_{2i-1},i=1,2,...,n_2 \\ Y_{i+n_1}=Y_{i+n_1}+x_{2i},i=1,2,...,n_2 \\ dis(i,j)=\sqrt{(X_i-X_j)^2+(Y_i-Y_j)^2} \end{cases}$

模型求解
采用实数编码，可行解本身就是编码，采用单点交叉，基本位变异的方法进行遗传迭代。

clear
global P_x
global P_y
global G_1
global n
global n1
global cost_line
global cost_switch
load('final_graph');
P_x = P_i;
P_y = P_j;
G_1 = G_2;
n = size(P_i, 2);
n1 = 50;

% GA
w = 3000; % 种群的个数
g = 10000; % 进化代数
rand('state', sum(clock));  % 初始化随机数发生器

% 编码长度
n0 = 2 * (n - n1);

% 获取初始种群
for k = 1 : w
    J(k, :) = rand(n0, 1) * 100 - 50;
end

for k = 1 : g   % 进化轮数
    A = J; % J: 父代 A: 交叉遗传产生的子代 B：从A变异而来的子代
    c = randperm(w);

    % 交叉（融合最优解）
    for i = 1 : 2 : w
        F = 1 + floor((n0 - 1) * rand);  % [1, n0-1]
        tmp = A(c(i), [F + 1 : n0]);
        A(c(i), [F + 1 : n0]) = A(c(i + 1), [F + 1 : n0]);
        A(c(i + 1), [F + 1 : n0]) = tmp;
    end
    
    % 变异
    by = [];
    while ~length(by)
        by = find(rand(1, w) < 0.95);
    end
    B = A(by, :);
    for i = 1 : length(by)
        bw = 1 + floor(n0 * rand);
        B(i, bw) = rand * 100 - 50;
    end
    
    G = [J; A; B];
    N = size(G, 1);
    
    % 优胜略汰
    Cost = zeros(1, N);
    for i = 1 : N
        Cost(i) = f1(G(i, :));
    end
    [Val, g_2] = sort(Cost);
    g_2 = g_2(1 : w);
    G_2 = G(g_2, :);
    J = G_2;    % 产生下一代
    disp(strcat('进化代数：', int2str(k), '； 当前最优为  ',int2str(Val(1))));
end

x = J(1, :);
fval = Val(1);

% 可视化
for i = n1 + 1 : n
    P_x(i) = P_x(i) + x(2 * (i - n1) - 1);
    P_y(i) = P_y(i) + x(2 * (i - n1));
end
figure
G_n = graph(G_1);
col = zeros(n, 3);
for i = 1 : n1
    col(i, :) = [0 0.4470 0.7410];
end
for i = n1 + 1 : n
    col(i, :) = [1 0 0];
end
h = plot(G_n, 'NodeColor', col);
h.XData = P_x;
h.YData = P_y;
h.LineWidth = 1.5 * G_n.Edges.Weight;


function cost = f1(x)
    global n
    global n1
    delta_x = zeros(n - n1);
    delta_y = zeros(n - n1);
    for i = 1 : n - n1
        delta_x(i) = x(2 * i - 1);
        delta_y(i) = x(2 * i);
    end
    cost = f(delta_x, delta_y);
end

% 评价函数
function cost = f(delta_x, delta_y)
	global P_x
    global P_y
    global G_1
    global n
    global n1
    global cost_line
    cost = 0;
    for i = n1 + 1 : n
        P_x(i) = P_x(i) + delta_x(i - n1);
        P_y(i) = P_y(i) + delta_y(i - n1);
    end
    for i = 1 : n
        for j = i + 1 : n
            if G_1(i, j) >= 1
                cost = cost + dis(i, j);
            end
        end
    end
    for i = n1 + 1 : n
        P_x(i) = P_x(i) - delta_x(i - n1);
        P_y(i) = P_y(i) - delta_y(i - n1);
    end
end

% 平面距离函数
function d = dis(i, j)
	global P_x
    global P_y
    d = sqrt((P_x(i) - P_x(j))^2 + (P_y(i) - P_y(j))^2);
end

结果可视化

请添加图片描述

TSP 问题
TSP问题中，可行解采用实数编码。我们知道，TSP 的可行解为一个顶点序列，除首尾外不重复。在变异和交叉操作中很难保证这种唯一性。实际上，采用实数编码，我们编码的是数之间的大小关系，顶点序列的第 k 位，对应的是排第 k 大数的索引。
举个例子，在一个有 5 个节点的图中，设起点为 0，我们考察其最小哈密顿圈。可行解的编码长度为 4（因为起点终点都已经固定，我们只需要确定中间的 4 个节点），我们用 4 个随机实数表示每一个解，如 $[0.2, 0.4, 0.8, 0.4]$ ，
顶点序列的第 1 位是排第 1 大的数（0.8）的索引为 3；
顶点序列的第 2 位是排第 2 大的数（0.4）的索引为 2（或者 4）；
顶点序列的第 3 位是排第 3 大的数（0.4）的索引为 4（或者 2）；
顶点序列的第 4 位是排第 4 大的数（0.2）的索引为 1；
因此，它表征的是顶点序列 $[3, 2, 4, 1]$ .
以上也是解码的过程，在 matlab 中 sort 方法的第二返回值就是该索引序列，因此可以非常简便地实现。

二进制编码

一串二进制代码对应一个解，稳定性高，种群多样性丰富，但需要存储空间大。最容易想到的例子是背包问题，它的解是某个集合的子集。值得注意的是交叉变异的后代并不都满足背包的体积约束，因此根据适应度值（装入背包中物体的总价值）选择个体之前要根据体积约束淘汰一些个体。

 
clear

load("b1.mat");
V = 1405406868;

w = 50; % 种群的个数
g = 50; % 进化代数
rand('state', sum(clock));  % 初始化随机数发生器


% 获取初始种群
for k = 1 : w
    V_sum = inf;
    while V_sum > V
        V_sum = 0;
        for i = 1 : n
            if rand > 0.6
                J(k, i) = 1;
                V_sum = V_sum + W(i);
            else
                J(k, i) = 0;
            end
        end
    end
end

for k = 1 : g   % 进化轮数
    A = J; % J: 父代 A: 交叉遗传产生的子代 B：从A变异而来的子代
    c = randperm(w);

    % 交叉（融合最优解）
    for i = 1 : 2 : w
        F = 1 + floor((n - 1) * rand);  % [1, n-1]
        tmp = A(c(i), [F + 1 : n]);
        A(c(i), [F + 1 : n]) = A(c(i + 1), [F + 1 : n]);
        A(c(i + 1), [F + 1 : n]) = tmp;
    end
    
    % 变异
    by = [];
    while ~length(by)
        by = find(rand(1, w) < 0.1); % rand(a, b): a*b矩阵，(0, 1)
    end
    B = A(by, :);
    for i = 1 : length(by)
        bw = 1 + floor(n * rand);
        B(i, bw) = 1 - B(i, bw);
    end
    
    G = [J; A; B];
    N = size(G, 1);
    
    % 优胜略汰
    % 淘汰不满足体积约束的个体
    V_sum = zeros(1, N);    % err: 必须重新初始化，否则沿用上次部分数据
    for i = 1 : N
        for j = 1 : n
            V_sum(i) = V_sum(i) + W(j) * G(i, j);
        end
    end
    g_1 = find(V_sum <= V);
    G_1 = G(g_1, :);
    N = size(G_1, 1);
    
    % 留下价值更大的个体
    C_sum = zeros(1, N);
    for i = 1 : N
        for j = 1 : n
            C_sum(i) = C_sum(i) + C(j) * G_1(i, j);
        end
    end
    [Val, g_2] = sort(C_sum, 'descend');
    g_2 = g_2(1 : 50);
    G_2 = G_1(g_2, :);
    J = G_2;    % 产生下一代
end

C_max = 0;
for j = 1 : n
    C_max = C_max + C(j) * J(1, j);
end
format long
C_max