数学建模优化算法——遗传算法

遗传算法原理

具体的原理网上还是比较多，我就不赘述了。这篇文章的主要目的是为了讲述一下自己从B站UP主连大数模上面获得一个遗传算法细致思路。主要会分享给大家这个matlab代码以及每个代码的含义。
不过为了便于大家理解我的思路，大致补充一些东西：

• 首先遗传算法是以种群进行优化的，每个种群中包含了多个个体，这个个体就是数学建模中模型（优化）的任意一个决策变量组合得到的目标函数值以及决策变量取值
• 其次优化算法多数只能得到一个局部解，所以求解方式大多数就是在循环迭代中找到一个相对较优的解。

遗传算法代码

主函数

clc
clear% 清空工作区，主函数
maxgen = 500;% 设定迭代次数
sizepop = 500;% 设定初始种群大小
pcross = 0.6;% 设定每个个体交叉的概率
pmutation = 0.01;% 设定每个个体变异的概率
lenchrom = [1, 1, 1, 1]; % 这里修改变量个数，给出决策变量
bound = [0.2, 0.8; 0.1, 0.6; 0.01, 1; 0.01, 0.25];  %这里修改变量约束范围，每个决策变量的取值范围

individuals = struct('fitness', zeros(1, sizepop), 'chrom', []);% 初始化种群（其实就是生成500个目标函数值以及对应的变量大小）
avgfitness = [];% 给一个数组记录每次迭代的种群内得到的平均目标函数值
bestfitness = [];% 给定一个数组记录每次迭代中中群内的最优目标函数值
bestchrom = [];% 给定一个数组记录每次迭代中中群内的最优目标函数值的对应决策变量值
for i = 1 : sizepop% 该循环是为了随机得到一个满足约束条件的种群
    individuals.chrom(i,:) = Code(lenchrom, bound);% 得到种群中每个个体对应的变量值
    x = individuals.chrom(i, :);
    individuals.fitness(i) = fun(x);% 得到种群中每个个体对应的目标函数值
end
[bestfitness, bestindex] = min(individuals.fitness);% 记录创建种群中的最小目标函数值和对应变量的下标
bestchrom = individuals.chrom(bestindex, :);% 根据下标找到对应的最优变量值
avgfitness = sum(individuals.fitness) / sizepop;% 记录该创建种群中的平均目标函数值

trace=[];% 创建空矩阵记录每次迭代对应最优目标函数值和平均目标函数值
for i = 1 : maxgen% 迭代500次
    individuals = Select(individuals, sizepop);% 对种群中的个体进行选择操作
    avgfitness = sum(individuals.fitness) / sizepop;% 个人感觉这个代码是多余的
    individuals.chrom = Cross(pcross, lenchrom, individuals.chrom, sizepop, bound);% 对种群中的个体基因（变量）进行交叉
    individuals.chrom = Mutation(pmutation, lenchrom, individuals.chrom, sizepop, [i, maxgen], bound);% 对种群中的个体基因（变量）进行变异

    for j = 1 : sizepop% 更新对应的种群中经过选择，交叉和变异之后每个个体的目标函数值
        x = individuals.chrom(j, :);
        individuals.fitness(j) = fun(x);
    end
    
    [newbestfitness, newbestindex] = min(individuals.fitness);% 记录该次迭代种群中最好（最优值求最小）的目标函数值和对应的下标
    [worestfitness, worestindex] = max(individuals.fitness);% 记录该次迭代种群中最差（大）（最优值求最小）的目标函数值和对应的下标
    if bestfitness > newbestfitness% 判段该次迭代种群中新生成的目标值是否比原来的最优值小
        bestfitness = newbestfitness;% 进行最优值的更新
        bestchrom = individuals.chrom(newbestindex, :);% 记录最优值对应的决策变量值
    end
    individuals.chrom(worestindex, :) = bestchrom;% 这个的意思就是把这个种群中最差的个体进行替换得到对应的变量值
    individuals.fitness(worestindex) = bestfitness;% 这个的意思就是把这个种群中最差的个体进行替换得到对应的目标函数值
    
    avgfitness=sum(individuals.fitness) / sizepop;% 记录该次迭代的平均值
    
    trace = [trace; avgfitness, bestfitness];% 记录轨迹
end

% 画图
figure
plot((1 : maxgen)', trace(:, 1), 'r-', (1: maxgen)', trace(:, 2), 'b--');
title(['函数值曲线  ' '终止代数＝' num2str(maxgen)], 'fontsize', 12);
xlabel('进化代数', 'fontsize', 12);
ylabel('函数值', 'fontsize', 12);
legend('各代平均值', '各代最佳值', 'fontsize', 12);
ylim([-0.5, 5])
disp('函数值                   变量');
disp([1 / bestfitness, x]);

对于主函数的讲解我直接放进代码块的注释了，这个地方主要讲解一下对应的模型如何修改代码。
首先比较关键的就是：

1. 迭代次数maxgen和**sizepop **是迭代次数和种群中的个体个数，后续调整的话，主要是为了减少算法复杂度和增加求解准确度。
2. 决策变量个数lenchrom这个变量（我个人是不太懂Matlab的）在我看来这个地方就是有多少个决策变量这个矩阵就创建多大（矩阵值应该都是1）。

生成个体中的每个决策变量值：Code函数

% 对于种群中的每个个体生成对应的变量（基因）
function ret = Code(lenchrom, bound)% 给出需要生成的变量个数和变量取值范围
flag = 0;
while flag == 0
    pick = rand(1, length(lenchrom));% 随机生成一个变量数组（0，1）范围内
    ret = bound(:, 1)' + (bound(:, 2) - bound(:, 1))' .* pick;%将其对应到变量的取值范围内
    flag = test(lenchrom, bound, ret);% 这是一个约束条件，检验生成的变量是否满足要求（避免极端情况的发生）
end

这个代码中主要就是用一个死循环和一个随机函数rand生成一个变量组合，这个值是一个在(0,1)范围内的值（再次声明不是很懂Matlab，不知道是不是能等于0，如果可以不等于0的话，后面很多代码的一个判断部分其实可以去除），通过每个变量的上下界及其距离可以得到个体的变量组合，下面会举一个例子，方便大家理解。其中test函数是后续会用到的检验函数，我觉得也可以后续约束条件函数，判断生成的变量组合是否满足约束条件，若满足则跳出循环。
$$\begin{gathered} \text{若}x\text{变量的取值范围为}\left[0,3\right]\text{，}pick\text{是}\left(0,1\right)\text{的值} \\ x\text{可以取值为}0+\left[\max \left(\text{取值范围}\right)-\min \left(\text{取值范围}\right)\right]pick \end{gathered}$$

检验函数（满足约束条件）：test函数

% 检验函数，我个人感觉如果需要约束条件的话，应该全部写在这个代码中是最好的选择
function flag = test(lenchrom, bound, code)% 给出变量个数和上下界以及生成的变量值
flag = 1;% 为了帮助跳出Code函数中的死循环
[n, ~] = size(code);% 给出生成变量的个数（个人理解代码可以简化一下）
for i = 1 : n% 对每一个个体的变量进行检验
    if code(i) < bound(i, 1) || code(i) > bound(i, 2)% 判断是否满足取值范围，后续也可以增加约束条件，若不满足，则为0
        flag = 0;
    end
end

这个代码对于**[n,~]**我的理解就是变量的长度是一个值，通过对个体中每个变量进行判断是否满足约束条件（在这里面是变量的取值范围），结合Code函数的话，最终会得到满足约束条件的决策变量组合也即种群中的一个个体。

选择函数，对于种群中的个体进行选择：Select函数

% 选择函数
function ret = Select(individuals, sizepop)% 给定种群的个体和种群个体大小

individuals.fitness = 1 ./ (individuals.fitness);% 由于需要求最小值，所以将目标函数值小的个体转换一下，方便以概率大的进行保存下来
sumfitness = sum(individuals.fitness);% 求全部的最优函数值
sumf = individuals.fitness ./ sumfitness;% 得到对应个体的概率
index=[];% 创建一个数组记录
for i = 1 : sizepop% 对种群中每一个个体进行选择
    pick = rand;% 生成随机数——概率
    while pick==0% 判断概率是否等于0
        pick = rand;% 重新生存
    end
    for j = 1 : sizepop% 在这个循环我个人感觉是比较难以理解，首先生成一个随机数，这个随机数概率选择这
        % 个创建种群中个体目标值较小，记录下来，意味着可能会有相同值，但是若这个概率比较小的话，
        % 按照统计学概率来说留下来的基本都是目标值相对较小的
        pick = pick - sumf(j);
        if pick < 0
            index = [index, j];% 取得一个值就跳出循环
            break;
        end
    end
end
individuals.chrom = individuals.chrom(index, :);% 找到选择的个体将其记录下来
individuals.fitness = 1 ./ individuals.fitness(index);
ret = individuals;% 输出选择后的种群

这个函数的目的主要是对一个种群中选择初始创建种群中的相对较优的个体，按照我的理解这个种群中个体会重复。

交叉函数对于种群中的任意两个个体基因进行交叉：Cross函数

% 进行个体中的变量交叉
function ret = Cross(pcross, lenchrom, chrom, sizepop, bound)% 给定交叉概率，变量个数和种群
% 中的全部个体，给定种群大小和变量范围
for i = 1 : sizepop% 对于种群中每个个体
    pick = rand(1, 2);% 生成两个随机值为了找到两个个体的下标位置
    while prod(pick) == 0
        pick = rand(1, 2);% 防止极端情况进行重新生成
    end
    index = ceil(pick .* sizepop);% 找到两个个体对应的下标
    pick = rand;% 生成一个概率值判断这两个个体是否可以进行交叉
    while pick == 0
        pick = rand;
    end
    if pick > pcross
        continue;
    end
    flag = 0;% 跳出code代码中的循环
    while flag == 0
        pick = rand;% 生成随机数为了找到交叉的变量位置
        while pick == 0
            pick = rand;
        end
        pos = ceil(pick .* sum(lenchrom));% 找到交叉的变量位置
        pick = rand;% 生成一个随机数，后续按照等比例对个体的变量进行交叉
        v1 = chrom(index(1), pos);% 招待两个个体的对应变量
        v2 = chrom(index(2), pos);
        chrom(index(1), pos) = pick * v2 + (1 - pick) * v1;% 对两个变量个体进行等比例交叉
        chrom(index(2), pos) = pick * v1 + (1 - pick) * v2;
        flag1 = test(lenchrom, bound, chrom(index(1), :));% 判断交叉之后的变量是否满足约束条件
        flag2 = test(lenchrom, bound, chrom(index(2), :));
        if flag1 * flag2 == 0
            flag = 0;
        else
            flag = 1;
        end% 满足则跳出循环
    end
end
ret = chrom;% 最终输出种群内进行交叉之后的个体

这个代码的目的就是选择种群中的任意两个个体找到对应的基因按照如下方式进行交叉。但是代码我感觉有一点问题，就是他是随机选择了两个个体，并不是对每个个体进行交叉变异，后续的话大家自己调整看一看。

进行交叉部分的代码

		v1 = chrom(index(1), pos);% 招待两个个体的对应变量
        v2 = chrom(index(2), pos);
        chrom(index(1), pos) = pick * v2 + (1 - pick) * v1;% 对两个变量个体进行等比例交叉
        chrom(index(2), pos) = pick * v1 + (1 - pick) * v2;

交叉代码的数学公式

$$\begin{gathered} pick\text{是一个随机}\left(0,1\right)\text{的数据，}{x}_{ij}\text{和}{x}_{mj}\text{是找到的第}i,m\text{个} \\ \text{个体上的第}j\text{个变量，通过交叉得到新的变量} \\ x{\prime }_{ij}=pick\times x_{mj}+\left(1-pick\right)\times x_{ij} \\ x{\prime }_{mj}=pick\times x_{ij}+\left(1-pick\right)\times x_{mj} \end{gathered}$$

变异函数对于种群中个体的基因进行变异：Mutation函数

% 变异操作
% 给定参数变异概率，变量个数，种群个体，种群大小和一个二维矩阵用于变异操作，以及变量对应范围
function ret = Mutation(pmutation, lenchrom, chrom, sizepop, pop, bound)
for i = 1 : sizepop% 对于种群中的每个个体
    pick = rand;% 创建一个随机数
    while pick == 0% 防止极端情况
        pick=rand;
    end
    index = ceil(pick * sizepop);% 这个代码应该也没有用到
    pick = rand;
    if pick > pmutation% 判断该个体是否有概率变异
        continue;
    end
    flag = 0;% 判断code函数是否可以跳出循环
    while flag == 0
        pick = rand;% 生成一个随机数为了找到该个体的变异基因位置
        while pick == 0
            pick = rand;
        end
        pos = ceil(pick * sum(lenchrom));% 找到个体变异基因位置
        v = chrom(i, pos);% 找到该个体的变异基因
        v1 = v - bound(pos, 1);% 弄出两个数值对基因进行变异
        v2 = bound(pos, 2) - v;
        pick = rand;% 随机对基因进行变异
        if pick > 0.5
            delta = v2 * (1 - pick ^ ((1 - pop(1) / pop(2)) ^ 2));
            chrom(i, pos) = v + delta;
        else
            delta=v1 * (1 - pick ^ ((1 - pop(1) / pop(2)) ^ 2));
            chrom(i, pos) = v - delta;
        end
        % 判断生成的基因是否满足约束条件
        flag = test(lenchrom, bound, chrom(i, :));
    end
end
ret = chrom;% 得到该种群中的每个个体

这个代码的目的就是对一个种群的个体上的基因进行变异，变异方式如下。

变异主要代码块

		v = chrom(i, pos);% 找到该个体的变异基因
        v1 = v - bound(pos, 1);% 弄出两个数值对基因进行变异
        v2 = bound(pos, 2) - v;
        pick = rand;% 随机对基因进行变异
        if pick > 0.5
            delta = v2 * (1 - pick ^ ((1 - pop(1) / pop(2)) ^ 2));
            chrom(i, pos) = v + delta;
        else
            delta=v1 * (1 - pick ^ ((1 - pop(1) / pop(2)) ^ 2));
            chrom(i, pos) = v - delta;
        end

$$\begin{gathered} v\text{是你找到的个体变异基因，}{v}_1\text{和}{v}_2\text{是变异需要用到的值} \\ {v}_1=v-\min \left(v\text{的取值范围}\right)\text{，}{v}_2=\max \left(v\text{的取值范围}\right)-v \\ \text{生成一个值}pick\left(0,1\right)\text{，判断}pick\text{是否}<0.5\text{进行不同的变异} \\ \{\begin{array}{l}\Delta v=v_2\left(1-pic{k}^{{\left(\frac{1-pop\left(1\right)}{pop\left(2\right)}\right)}^2}\right),v=v+\Delta vifpick<0.5\\ \Delta v=v_1\left(1-pic{k}^{{\left(\frac{1-pop\left(1\right)}{pop\left(2\right)}\right)}^2}\right),v=v-\Delta velse\end{array} \\ \text{其中}pop\left(1\right)\text{和}pop\left(2\right)\text{也是两个数值} \end{gathered}$$
我个人是不太懂这个变异方式的，但是按照我的理解就是随机在这个决策变量的范围内进行跳动。

目标函数fun函数

% 目标函数值，修改目标函数和对应参数即可
function y = fun(x)

y = 1 / (62.17 * x(2) * sqrt(2 * x(2) * x(1) * sqrt(x(3))) * exp(x(4)));

Sdown = (sqrt(3) / 4) * x(1) ^ 2 * 6;
Sup = x(3) * Sdown;
V = (1 / 3) * x(2) * (Sup + Sdown + sqrt(Sup * Sdown));

if V <= 1 && V >= 0.01
    y = y;
else
    y = y + 10000;
end

通过对这个函数传入决策变量组合，然后计算该目标函数的值，最终大家可以在这个函数中写出自己建立模型的对应目标函数。