【数学建模】《实战数学建模：例题与讲解》第十四讲-模拟退火、遗传算法（含Matlab代码）

【数学建模】《实战数学建模：例题与讲解》第十四讲-模拟退火、遗传算法（含Matlab代码）

• • 基本概念
  • • 模拟退火（Simulated Annealing）
    • 遗传算法（Genetic Algorithms）
  • 习题14.1（1）
  • • 1. 题目要求
    • 2.解题过程——模拟退火算法
    • 3.程序
    • 4.结果
  • 习题14.1（2）
  • • 1. 题目要求
    • 2.解题过程——遗传算法
    • 3.程序
    • 4.结果
  • 习题14.2（1）
  • • 1. 题目要求
    • 2.解题过程——模拟退火算法
    • 3.程序
    • 4.结果
  • 习题14.2（2）
  • • 1. 题目要求
    • 2.解题过程
    • 3.程序——遗传算法
    • 4.结果

本系列侧重于例题实战与讲解，希望能够在例题中理解相应技巧。文章开头相关基础知识只是进行简单回顾，读者可以搭配课本或其他博客了解相应章节，然后进入本文正文例题实战，效果更佳。

如果这篇文章对你有帮助，欢迎点赞与收藏~

基本概念

现代优化算法，自20世纪80年代初开始流行，主要包括了一系列基于启发式原理的算法。这些算法，如禁忌搜索（Tabu Search）、模拟退火（Simulated Annealing）、遗传算法（Genetic Algorithms）和人工神经网络（Neural Networks），被广泛应用于解决各种实际应用问题。它们在理论研究和实际应用领域均取得了显著发展。

尽管这些算法的产生背景各异，它们共同的目标是寻找NP-hard问题的全局最优解。NP-hard问题的复杂性使得这些算法通常以启发式的方式寻找问题的实际解决方案。

启发式算法的范畴非常广泛，包括了针对复杂优化问题设计的蚁群算法（Ant Colony Algorithms）。某些启发式算法是针对特定实际问题而开发的，例如那些解决解空间分解、解空间限制等问题的算法。另外，还有一类算法是集成算法，它们结合了多种启发式算法的优点。

在解决组合优化问题方面，如旅行商问题（Traveling Salesman Problem, TSP）、二次分配问题（Quadratic Assignment Problem, QAP）和作业车间调度问题（Job-shop Scheduling Problem, JSP）等，现代优化算法表现出色。这些算法不仅在理论上有其独特之处，也在实际应用中展现了强大的效能。

模拟退火（Simulated Annealing）

模拟退火是一种概率型优化算法，其灵感来源于金属加热后再缓慢冷却的退火过程。在这个过程中，金属原子会随着温度的降低逐渐趋于稳定的状态，最终达到能量最低的结构。

工作原理：

1. 算法开始时设定一个较高的“温度”，随着迭代过程逐步降低。
2. 在每个温度下，算法通过随机选取解的邻域并计算其能量（或成本）来探索解空间。
3. 如果邻域解的能量低于当前解，或者即便能量高但满足一定的概率条件（依赖于温度和能量差），算法也会接受这个新解。
4. 随着温度的逐渐降低，接受劣质解的概率降低，算法趋于稳定。

遗传算法（Genetic Algorithms）

遗传算法是一种模拟生物进化过程的搜索启发式算法。它基于自然选择的原理，通过生物遗传机制中的交叉（crossover）、变异（mutation）和选择（selection）等操作进行优化。

工作原理：

1. 算法初始化时生成一个包含多个候选解的种群。
2. 每个解（通常称为个体）都有一个与之相关的适应度值，用于评估其优劣。
3. 算法通过选择过程保留适应度高的个体，并通过交叉和变异操作生成新个体。
4. 交叉操作模拟了染色体的交换，而变异则是对个体的随机调整。
5. 经过多代的进化，种群中的个体逐渐趋于最优解。

习题14.1（1）

1. 题目要求

假设有12件物品，质量和价值如下表所示。包的最大允许质量是46公斤。

试使用模拟退火算法和遗传算法求解包中可装载物品的最大价值。

2.解题过程——模拟退火算法

解：

记12件物品的的质量和价值为 最大容量为 。

（1）解空间

解空间可以写为:

当 时，表示选择第 个物品，否则表示不选。

（2）目标函数

目标函数为最终的物品价值。我们求价值最大值，首先通过转换求价值相反数的最小值，即：

一次迭代由下列三步产生。

（3）新解的产生

任选序号 ，将取反，此时新的选取方法为：

计算此时的重量：

若 则新解有效，否则重新生成。

（4）代价函数差

路径差可以表示为：

（5）接受准则

（6）降温

选定降温系数 降温，取新温度 ，初始温度 。

（7）结束条件

选定终止温度 判断退火是否结束，当 时，结束模拟，输出当前状态。

3.程序

求解的MATLAB程序如下：

clc, clear

mass = [2, 5, 18, 3, 2, 5, 10, 4, 11, 7, 14, 6]; % 物品质量
value = [5, 10, 13, 4, 3, 11, 13, 10, 8, 16, 7, 4]; % 物品价值
solution = zeros(1, length(mass)); % 初始化解
max_mass = 46; % 最大允许质量
min_temperature = 0.1^30;
iterations = 20000;
alpha = 0.999;
temperature = 1;

for k = 1:iterations
    old_value = -sum(solution.*value);
    temp_solution = solution;
    while 1
        item = 1 + floor(length(mass)*rand(1, 2));
        temp_solution(item) = ~temp_solution(item); % 改变选取物品的状态
        if sum(temp_solution.*mass) > max_mass % 如果超过背包最大允许质量，重新选取
            temp_solution = solution;
            continue
        else
            break
        end
    end
    new_value = -sum(temp_solution.*value);
    df = new_value - old_value;
    if df < 0 || exp(-df/temperature) >= rand % 接受新解
        solution = temp_solution;
    end
    temperature = temperature * alpha; % 降温
    if temperature < min_temperature
        break
    end
end

solution, total_mass = sum(solution.*mass), best_value = sum(solution.*value)

4.结果

得到结果为：

即选取第1，2，5，6，7，8，9，10件物品，总重量为46，价值为76。

习题14.1（2）

1. 题目要求

同上。

2.解题过程——遗传算法

解：

设定种群大小 ，最大迭代次数 ，交叉率 ，变异率 。

基因采取编码：

使用改良圈算法产生 个可行解，转化为初始基因编码，目标函数为最终的物品价值即

接下来每一次的迭代都以概率 进行基因重组和基因变异，最后选择目标函数值最大的 个个体进化到下一代。

3.程序

求解的MATLAB程序如下：

clc, clear

% 输入数据
weights = [2, 5, 18, 3, 2, 5, 10, 4, 11, 7, 14, 6];
values = [5, 10, 13, 4, 3, 11, 13, 10, 8, 16, 7, 4];

num_items = length(weights);
cross_rate = .95; % 交叉概率
mutation_rate = .15; % 变异概率
max_iter = 30; % 最大迭代次数
pop_size = 100; % 种群大小
max_capacity = 46; % 背包最大承重

% 构建初始种群
population = init_population(pop_size, num_items, weights, max_capacity);

% 初始化统计指标
avg_values = zeros(1, max_iter);
max_values = zeros(1, max_iter);

% 迭代进化
for i = 1:max_iter
    % 交叉
    offspring_cross = crossover(population, cross_rate, weights, max_capacity);
    % 变异
    offspring_mutate = mutation(population, mutation_rate, weights, max_capacity);
    % 选择
    population = selection([population; offspring_cross; offspring_mutate], pop_size, values);
    % 统计当前迭代的平均价值和最大价值
    avg_values(i) = mean(sum(population*values', 2));
    max_values(i) = max(sum(population*values', 2));
end

% 输出结果
best_solution = population(1, :)
best_value = sum(best_solution*values')
best_mass = sum(best_solution*weights')

% ************************* 以下为函数实现 *************************
% 初始化种群函数
function population = init_population(pop_size, num_items, weights, max_capacity)
population = zeros(pop_size, num_items);
for i = 1:pop_size
    chromosome = round(rand(1, num_items));
    while sum(chromosome.*weights) > max_capacity
        chromosome = round(rand(1, num_items));
    end
    population(i, :) = chromosome;
end
end

% 交叉函数
function offspring = crossover(population, cross_rate, weights, max_capacity)
offspring = population;
num_individuals = size(population, 1);
for i = 1:2:num_individuals
    if rand < cross_rate
        cross_point = ceil(rand * size(population, 2));
        temp1 = offspring(i, :);
        temp2 = offspring(i+1, :);
        temp1(cross_point) = temp2(cross_point);
        temp2(cross_point) = offspring(i, cross_point);
        if (temp1 * weights' <= max_capacity) && (temp2 * weights' <= max_capacity)
            offspring(i, :) = temp1;
            offspring(i+1, :) = temp2;
        end
    end
end
end

% 变异函数
function offspring = mutation(population, mutation_rate, weights, max_capacity)
offspring = population;
for i = 1:size(population, 1)
    if rand < mutation_rate
        mutate_point = ceil(rand * size(population, 2));
        temp = offspring(i, :);
        temp(mutate_point) = ~temp(mutate_point);
        if sum(temp.*weights) <= max_capacity
            offspring(i, :) = temp;
        end
    end
end
end

% 选择函数
function new_population = selection(population, pop_size, values)
fitness_values = sum(population*values', 2);
[~, sorted_indices] = sort(fitness_values, 'descend');
new_population = population(sorted_indices(1:pop_size), :);
end

4.结果

得到结果为：

即选取第1，2，5，6，7，8，9，10件物品，总重量为46，价值为76。

习题14.2（1）

1. 题目要求

假设有一个旅行商人要拜访全国 31 个省会城市，他需要选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。

全国 31 个省会城市的坐标如下表所示。

试使用模拟退火算法和遗传算寻找最短路径。

2.解题过程——模拟退火算法

解：

设城市 之间的距离 。

（1）解空间

解空间可以表示为：
KaTeX parse error: Undefined control sequence: \mbox at position 68: …2,…,\pi_{31})\̲m̲b̲o̲x̲{为}\{1,2,…,31…
特别的，规定 。

初始解可以选择为 。

（2）目标函数

目标函数为路径长度：

一次迭代由下列三步产生.

（3）新解的产生

任选序 ，交换 之间的顺序交换，此时新路径为：

（4）代价函数差

（5）接受准则

（6）降温

选定降温系数 降温，取新温度 ，初始温度 .

（7）结束条件

选定终止温度 ，当时，结束模拟。

3.程序

求解的MATLAB程序如下：

clc, clear

% 城市坐标
city_coordinates = [1304, 2312; 3639, 1315; 4177, 2244; 3712, 1399; 3488, 1535; 3326, 1556; ...
    3238, 1229; 4196, 1004; 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, 2367; ...
    3394, 2643; 3439, 3201; 2935, 3240; 3140, 3550; 2545, 2357; 2778, 2826; ...
    2370, 2975];

num_cities = size(city_coordinates, 1);

% 初始化路径
path = 1:num_cities;

% 计算距离矩阵
distances = pdist2(city_coordinates, city_coordinates);

% 初始化总距离
total_distance = sum(distances(sub2ind(size(distances), path, [path(2:end), path(1)])));

% 设置退火参数
initial_temperature = 100;
alpha = 0.9999;
final_temperature = 0.1^30;

temperature = initial_temperature;

while temperature > final_temperature
    % 随机选择两个城市
    cities_to_swap = sort(randperm(num_cities, 2));
    
    % 生成新路径
    new_path = path;
    new_path(cities_to_swap(1):cities_to_swap(2)) = new_path(cities_to_swap(2):-1:cities_to_swap(1));
    
    % 计算新距离
    new_distance = sum(distances(sub2ind(size(distances), new_path, [new_path(2:end), new_path(1)])));
    
    % 判断是否接受新解
    if new_distance < total_distance || rand < exp((total_distance - new_distance) / temperature)
        path = new_path;
        total_distance = new_distance;
    end
    
    % 降温
    temperature = temperature * alpha;
end

% 输出结果
disp('Optimal path:');
disp(path);
disp('Total distance:');
disp(total_distance);

% 绘制图像
plot(city_coordinates([path, path(1)], 1), city_coordinates([path, path(1)], 2), '-o');

4.结果

最终得到的结果如图，路径长度 。

习题14.2（2）

1. 题目要求

同上。

2.解题过程

解：

设定种群大小 ，最大迭代次数 。

交叉率 。

变异率 。

基因编码用随机数列 其中 在整个序列中的升序排序位置为城市 所在的位置。

使用改良圈算法产生个可行解，转化为初始基因编码，目标函数为最终的物品价值即：

接下来每一次的迭代都以概率 进行基因重组和基因变异，最后选择目标函数值最大的 个个体进化到下一代。

3.程序——遗传算法

求解的MATLAB程序如下：

clc, clear

distanceMatrix = zeros(31);
optimalPath = zeros(1, 31);
coordinates = [1304, 2312; 3639, 1315; 4177, 2244; 3712, 1399; 3488, 1535; 3326, 1556; ...
    3238, 1229; 4196, 1004; 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, 2367; ...
    3394, 2643; 3439, 3201; 2935, 3240; 3140, 3550; 2545, 2357; 2778, 2826; ...
    2370, 2975];

for i = 1:31
    optimalPath(i) = i;
    for j = 1:31
        distanceMatrix(i, j) = sqrt((coordinates(i, 1) - coordinates(j, 1))^2+(coordinates(i, 2) - coordinates(j, 2))^2);
    end
end

distance = 0;
for i = 1:30
    distance = distance + distanceMatrix(optimalPath(i), optimalPath(i+1));
end

distance = distance + distanceMatrix(optimalPath(1), optimalPath(31));
w = 100;
g = 100; 
rand('state', sum(clock)); % 初始化随机数发生器

for k = 1:w  % 通过改良圈算法选取初始种群
    c1 = randperm(31);  % 产生1,..,31的一个全排列
    for t = 1:100  % 该层循环是修改圈
        flag = 0;  % 修改圈退出标志
        for m = 1:31
            for n = m + 1:31
                m1 = m - 1;
                n1 = n + 1;
                if m == 1 && n == 31
                    continue
                end
                if m1 == 0
                    m1 = 31;
                end
                if n1 == 32
                    n1 = 1;
                end
                if distanceMatrix(c1(m1), c1(n)) + distanceMatrix(c1(m), c1(n1)) < ...
                        distanceMatrix(c1(m), c1(m1)) + distanceMatrix(c1(n), c1(n1))
                    c1 = [c1(1:m-1), c1(n:-1:m), c1(n+1:31)];
                    flag = 1; % 修改圈
                end
            end
        end
        if flag == 0
            chromosomeMatrix(k, c1) = 1:31;
            break % 记录下较好的解并退出当前层循环
        end
    end
end

chromosomeMatrix(:, 1) = 0;
chromosomeMatrix = chromosomeMatrix / 31; % 把整数序列转换成［0,1］区间上的实数，即转换成染色体编码
for k = 1:g % 该层循环进行遗传算法的操作
    childPopulation = chromosomeMatrix; % 交配产生子代 A 的初始染色体
    c = randperm(w); % 产生下面交叉操作的染色体对
    for i = 1:2:w
        F = 1 + floor(31*rand(1)); % 产生交叉操作的地址
        temp = childPopulation(c(i), [F:31]); % 中间变量的保存值
        childPopulation(c(i), [F:31]) = childPopulation(c(i+1), [F:31]); % 交叉操作
        childPopulation(c(i+1), F:31) = temp;
    end
    by = []; % 为了防止下面产生空地址，这里先初始化
    while ~length(by)
        by = find(rand(1, w) < 0.15); % 产生变异操作的地址
    end
    mutationList = childPopulation(by, :); % 产生变异操作的初始染色体
    for j = 1:length(by)
        bw = sort(1+floor(31*rand(1, 3))); % 产生变异操作的3个地址
        mutationList(j, :) = mutationList(j, [1:bw(1) - 1, bw(2) + 1:bw(3), bw(1):bw(2), bw(3) + 1:31]); % 交换位置
    end
    G = [chromosomeMatrix; childPopulation; mutationList];  % 父代和子代种群合在一起
    [SG, indl] = sort(G, 2);
    populationSize = size(G, 1);
    pathLengths = zeros(1, populationSize);  % 路径长度的初始值
    for j = 1:populationSize
        for i = 1:31
            if i == 31
                pathLengths(j) = pathLengths(j) + distanceMatrix(indl(j, i), indl(j, 1));
            else
                pathLengths(j) = pathLengths(j) + distanceMatrix(indl(j, i), indl(j, i+1)); % 计算每条路径长度
            end
        end
    end
    [slong, ind2] = sort(pathLengths); % 对路径长度从小到大排序
    chromosomeMatrix = G(ind2(1:w), :); % 精选前 w 个较短的路径对应的染色体
end
optimalPath = indl(ind2(1), :), optimalPathLength = slong(1)
plot([coordinates(optimalPath, 1); coordinates(optimalPath(1), 1)], [coordinates(optimalPath, 2); coordinates(optimalPath(1), 2)], '- o')

4.结果

最终得到的结果如图，路径长度 。

如果这篇文章对你有帮助，欢迎点赞与收藏~