模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。模拟退火算法属于启发式搜索算法的一种。本文主要借鉴B站数学建模清风的视频，文章仅做笔记使用，侵权立删。

目录

爬山法寻找函数最大值

爬山法的算法搜索流程：

模拟退火的算法搜索流程：

Ct的设置：

题目1：求一个给定函数的最值问题

题目2：旅行商问题

题目3：书店买书问题

问题4：背包问题 

爬山法寻找函数最大值

爬山法的算法搜索流程：

  1. 根据初始解的位置，我们在下一步有两种走法:向左邻域走一步或者向右邻域走一步（走的步长越小越好，但会加大计算量);
  2. 比较不同走法下目标函数的大小，并决定下一步往哪个方向走。上图中向左走目标函数较大，因此我们应该往左走一小步;
  3. 不断重复这个步骤，直到找到一个极大值点(或定义域边缘处)，此时我们就结束搜索。
     

模拟退火的算法搜索流程：

  1. 随机生成一个解A，计算解A对应的目标函数值f(A)。
  2. 在A附近随机生成一个解B，计算解B对应的目标函数值f(B)
  3. 如果f(B)>f(A)，则将解B赋值给解A，然后重复上面步骤（爬山法的思想);如果f(B)≤f(A)，那么我们计算接受B的概率p_{t}=e^{​{-|f(B)-f(A)|}\times C_{t}} ,然后我们生成一个[0,1]之间的随机数r,如果r<p，我们就将解B赋值给解A,然后重复上面的步骤;否则我们返回第(2)步，在原来的A附近再重新生成一个解B，然后继续下去。·

Ct的设置：

定义初始温度T_0=100，温度下降公式为T_{t+1}=\alpha T_t，\alpha常取0.95，那么时刻t时的温度：T_{t}=\alpha^{t}T_{0}=100\times0.95^{t}。

取C_t=\frac{1}{T_t}=\frac{1}{100\times0.95^{t}}，那么P_t=e^{-\frac{f(B)-f(A)}{T_t}}=e^{-\frac{f(B)-f(A)}{100\times0.95^{t}}}。

注：这里取倒数是为了保证C_t关于t递增。

 （上图从右往左看就是退火过程，\Delta f=|f(B)-f(A)|）

  • 温度一定时，△f越小，概率越大，即目标函数相差越小接受的可能性越大。
  • Af一定时，温度越高，概率越大，即搜索前期温度较高时更有可能接受新解。 

题目1：求一个给定函数的最值问题

求函数y=11sinx+7cox(5x)在[-3,3]内的最大值

利用模拟退火算法求解新解的过程中，需要在原解A附近产生一个新解B。我们以n元函数求最值为例，来介绍新解产生的规则。

假设当前解为：(x_1,x_2,...,x_n)

首先生成一组随机数：(y_1,y_2,...,y_n),其中y_i服从N(0,1)

接下来计算(z_1,z_2,...,z_n),其中z_i=\frac{y_i}{\sqrt{y_{1}^{2}+y_{2}^{2}+...+y_{n}^{2}}}，对于每个i\in [1,n]，进行下面操作：

  1. 计算：x_{i}^{new}=x_i+T\times{}z_i，T是当前温度（也可以使用x_{i}^{new}=x_i+\sqrt{T}\times{z_i}这个公式）
  2. 判断：x_{i}^{new}是否位于上下界内，即是否满足lb_i\leqslant x_{i}^{new}\leqslant{ub_i}这个约束
  •         如果lb_i\leqslant x_{i}^{new}\leqslant{ub_i}满足，则直接令x_j=x_{j}^{new}
  •         如果x_{i}^{new}<lb_i，则x_j=r\times{lb_i}+(1-r)\times{x_i}，r\in[0,1]
  •         如果x_{i}^{new}>ub_i，则x_j=r\times{ub_i}+(1-r)\times{x_i}，r\in[0,1]

%% SA 模拟退火: 求解函数y = 11*sin(x) + 7*cos(5*x)在[-3,3]内的最大值(动画演示）
tic
clear; clc

%% 绘制函数的图形
x = -3:0.1:3;
y = 11*sin(x) + 7*cos(5*x);
figure
plot(x,y,'b-')
hold on  % 不关闭图形，继续在上面画图

%% 参数初始化
narvs = 1; % 变量个数
T0 = 100;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 200;  % 最大迭代次数
Lk = 100;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数
x_lb = -3; % x的下界
x_ub = 3; % x的上界

%%  随机生成一个初始解
x0 = zeros(1,narvs);
for i = 1: narvs
    x0(i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(1);    
end
y0 = Obj_fun1(x0); % 计算当前解的函数值
h = scatter(x0,y0,'*r');  % scatter是绘制二维散点图的函数（这里返回h是为了得到图形的句柄，未来我们对其位置进行更新）

%% 定义一些保存中间过程的量，方便输出结果和画图
max_y = y0;     % 初始化找到的最佳的解对应的函数值为y0
MAXY = zeros(maxgen,1); % 记录每一次外层循环结束后找到的max_y (方便画图）

%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  % 内循环，在每个温度下开始迭代
        y = randn(1,narvs);  % 生成1行narvs列的N(0,1)随机数
        z = y / sqrt(sum(y.^2)); % 根据新解的产生规则计算z
        x_new = x0 + z*T; % 根据新解的产生规则计算x_new的值
        % 如果这个新解的位置超出了定义域，就对其进行调整
        for j = 1: narvs
            if x_new(j) < x_lb(j)
                r = rand(1);
                x_new(j) = r*x_lb(j)+(1-r)*x0(j);
            elseif x_new(j) > x_ub(j)
                r = rand(1);
                x_new(j) = r*x_ub(j)+(1-r)*x0(j);
            end
        end
        x1 = x_new;    % 将调整后的x_new赋值给新解x1
        y1 = Obj_fun1(x1);  % 计算新解的函数值
        if y1 > y0    % 如果新解函数值大于当前解的函数值
            x0 = x1; % 更新当前解为新解
            y0 = y1;
        else
            p = exp(-(y0 - y1)/T); % 根据Metropolis准则计算一个概率
            if rand(1) < p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                x0 = x1; % 更新当前解为新解
                y0 = y1;
            end
        end
        % 判断是否要更新找到的最佳的解
        if y0 > max_y  % 如果当前解更好，则对其进行更新
            max_y = y0;  % 更新最大的y
            best_x = x0;  % 更新找到的最好的x
        end
    end
    MAXY(iter) = max_y; % 保存本轮外循环结束后找到的最大的y
    T = alfa*T;   % 温度下降
    pause(0.01)  % 暂停一段时间(单位：秒)后再接着画图
    h.XData = x0;  % 更新散点图句柄的x轴的数据（此时解的位置在图上发生了变化）
    h.YData = Obj_fun1(x0); % 更新散点图句柄的y轴的数据（此时解的位置在图上发生了变化）
end

disp('最佳的位置是：'); disp(best_x)
disp('此时最优值是：'); disp(max_y)

pause(0.5)
h.XData = [];  h.YData = [];  % 将原来的散点删除
scatter(best_x,max_y,'*r');  % 在最大值处重新标上散点
title(['模拟退火找到的最大值为', num2str(max_y)])  % 加上图的标题

%% 画出每次迭代后找到的最大y的图形
figure
plot(1:maxgen,MAXY,'b-');
xlabel('迭代次数');
ylabel('y的值');
toc

 

 

题目2：旅行商问题

给定一系列城市和每对城市之间的距离，求解访问每一座城市一次并回到起始城市的最短回路。

function  result =  calculate_tsp_d(path,d)
% 输入：path:路径（1至n的一个序列），d：距离矩阵
    n = length(path);
    result = 0; % 初始化该路径走的距离为0
    for i = 1:n-1  
        result = d(path(i),path(i+1)) + result;  % 按照这个序列不断的更新走过的路程这个值
    end   
    result = d(path(1),path(n)) + result;  % 别忘了加上从最后一个城市返回到最开始那个城市的距离
end
function path1 = gen_new_path(path0)
    % path0: 原来的路径
    n = length(path0);
    % 随机选择两种产生新路径的方法
    p1 = 0.33;  % 使用交换法产生新路径的概率
    p2 = 0.33;  % 使用移位法产生新路径的概率
    r = rand(1); % 随机生成一个[0 1]间均匀分布的随机数
    if  r< p1    % 使用交换法产生新路径 
        c1 = randi(n);   % 生成1-n中的一个随机整数
        c2 = randi(n);   % 生成1-n中的一个随机整数
        path1 = path0;  % 将path0的值赋给path1
        path1(c1) = path0(c2);  %改变path1第c1个位置的元素为path0第c2个位置的元素
        path1(c2) = path0(c1);  %改变path1第c2个位置的元素为path0第c1个位置的元素
    elseif r < p1+p2 % 使用移位法产生新路径
        c1 = randi(n);   % 生成1-n中的一个随机整数
        c2 = randi(n);   % 生成1-n中的一个随机整数
        c3 = randi(n);   % 生成1-n中的一个随机整数
        sort_c = sort([c1 c2 c3]);  % 对c1 c2 c3从小到大排序
        c1 = sort_c(1);  c2 = sort_c(2);  c3 = sort_c(3);  % c1 < = c2 <=  c3
        tem1 = path0(1:c1-1);
        tem2 = path0(c1:c2);
        tem3 = path0(c2+1:c3);
        tem4 = path0(c3+1:end);
        path1 = [tem1 tem3 tem2 tem4];
    else  % 使用倒置法产生新路径
        c1 = randi(n);   % 生成1-n中的一个随机整数
        c2 = randi(n);   % 生成1-n中的一个随机整数
        if c1>c2  % 如果c1比c2大，就交换c1和c2的值
            tem = c2;
            c2 = c1;
            c1 = tem;
        end
        tem1 = path0(1:c1-1);
        tem2 = path0(c1:c2);
        tem3 = path0(c2+1:end);
        path1 = [tem1 fliplr(tem2) tem3];   %矩阵的左右对称翻转 fliplr，上下对称翻转 flipud
    end
end
%% 模拟退火解决TSP问题
tic
rng('shuffle')  % 控制随机数的生成，否则每次打开matlab得到的结果都一样
clear;clc
% TSP数据集网站(http://www.tsp.gatech.edu/world/djtour.html) 上公测的最优结果6656。
coord = [11003.611100,42102.500000;11108.611100,42373.888900;11133.333300,42885.833300;
    11155.833300,42712.500000;11183.333300,42933.333300;11297.500000,42853.333300;
    11310.277800,42929.444400;11416.666700,42983.333300;11423.888900,43000.277800;
    11438.333300,42057.222200;11461.111100,43252.777800;11485.555600,43187.222200;
    11503.055600,42855.277800;11511.388900,42106.388900;11522.222200,42841.944400;
    11569.444400,43136.666700;11583.333300,43150.000000;11595.000000,43148.055600;
    11600.000000,43150.000000;11690.555600,42686.666700;11715.833300,41836.111100;
    11751.111100,42814.444400;11770.277800,42651.944400;11785.277800,42884.444400;
    11822.777800,42673.611100;11846.944400,42660.555600;11963.055600,43290.555600;
    11973.055600,43026.111100;12058.333300,42195.555600;12149.444400,42477.500000;
    12286.944400,43355.555600;12300.000000,42433.333300;12355.833300,43156.388900;
    12363.333300,43189.166700;12372.777800,42711.388900;12386.666700,43334.722200;
    12421.666700,42895.555600;12645.000000,42973.333300];
n = size(coord,1);  % 城市的数目

figure  % 新建一个图形窗口
plot(coord(:,1),coord(:,2),'o');   % 画出城市的分布散点图
hold on % 等一下要接着在这个图形上画图的

d = zeros(n);   % 初始化两个城市的距离矩阵全为0
for i = 2:n  
    for j = 1:i  
        coord_i = coord(i,:);   x_i = coord_i(1);     y_i = coord_i(2);  % 城市i的横坐标为x_i，纵坐标为y_i
        coord_j = coord(j,:);   x_j = coord_j(1);     y_j = coord_j(2);  % 城市j的横坐标为x_j，纵坐标为y_j
        d(i,j) = sqrt((x_i-x_j)^2 + (y_i-y_j)^2);   % 计算城市i和j的距离
    end
end
d = d+d';   % 生成距离矩阵的对称的一面

%% 参数初始化
T0 = 1000;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 1000;  % 最大迭代次数
Lk = 500;  % 每个温度下的迭代次数
alpfa = 0.95;  % 温度衰减系数

%%  随机生成一个初始解
path0 = randperm(n);  % 生成一个1-n的随机打乱的序列作为初始的路径
result0 = calculate_tsp_d(path0,d); % 调用我们自己写的calculate_tsp_d函数计算当前路径的距离

%% 定义一些保存中间过程的量，方便输出结果和画图
min_result = result0;     % 初始化找到的最佳的解对应的距离为result0
RESULT = zeros(maxgen,1); % 记录每一次外层循环结束后找到的min_result (方便画图）

%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  %  内循环，在每个温度下开始迭代
        path1 = gen_new_path(path0);  % 调用我们自己写的gen_new_path函数生成新的路径
        result1 = calculate_tsp_d(path1,d); % 计算新路径的距离
        %如果新解距离短，则直接把新解的值赋值给原来的解
        if result1 < result0    
            path0 = path1; % 更新当前路径为新路径
            result0 = result1; 
        else
            p = exp(-(result1 - result0)/T); % 根据Metropolis准则计算一个概率
            if rand(1) < p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                path0 = path1;  % 更新当前路径为新路径
                result0 = result1; 
            end
        end
        % 判断是否要更新找到的最佳的解
        if result0 < min_result  % 如果当前解更好，则对其进行更新
            min_result = result0;  % 更新最小的距离
            best_path = path0;  % 更新找到的最短路径
        end
    end
    RESULT(iter) = min_result; % 保存本轮外循环结束后找到的最小距离
    T = alpfa*T;   % 温度下降       
end

disp('最佳的方案是：'); disp(mat2str(best_path))
disp('此时最优值是：'); disp(min_result)

best_path = [best_path,best_path(1)];   % 在最短路径的最后面加上一个元素，即第一个点（我们要生成一个封闭的图形）
n = n+1;  % 城市的个数加一个（紧随着上一步）
for i = 1:n-1 
    j = i+1;
    coord_i = coord(best_path(i),:);   x_i = coord_i(1);     y_i = coord_i(2); 
    coord_j = coord(best_path(j),:);   x_j = coord_j(1);     y_j = coord_j(2);
    plot([x_i,x_j],[y_i,y_j],'-b')    % 每两个点就作出一条线段，直到所有的城市都走完
%     pause(0.02)  % 暂停0.02s再画下一条线段
    hold on
end

%% 画出每次迭代后找到的最短路径的图形
figure
plot(1:maxgen,RESULT,'b-');
xlabel('迭代次数');
ylabel('最短路径');
toc

 

 

题目3：书店买书问题

（假设有m个书店，n本书，那么买书的方案为m^n种）

某同学要从六家线上商城选购五本书籍B1,B2,B3,B4,B5，每本书籍在不同商家的售价以及每个商家的单次运费如下表所示，请给该同学制定最省钱的选购方案。
(注:在同一个店买多本书也只会收取一次运费)

M-每本书在每家店的价格
B1B2B3B4B5B6B7B8B9B10B11B12B13B14B15B16B17B18B19B20
A13131412125282334382938333224232023262132
A24027382623292422372932343127312226272527
A33525412026213724342242313726282323212628
A43326222938253432342427252631393421214134
A53329362421243328252924262629372425253227
A62532202120324222332435283826342139254023
A73522352929263830272125303332303225232623
A83622392634253223352920323431252438252925
A93223222127222030272441273327292231262524
A102728363238272933292529333425242237274230
A113928262737282331352730282032312132314321
A122228383340234330353423263623342440244130
A133025273227304027362230292132413333293131
A143421272925213633212821303522222440272523
A153127242539234030222838312129212540223135
freight——每家店的运费
A1A2A3A4A5A6A7A8A9A10A11A12A13A14A15
10101471251081491261159
function  money =  calculate_money(way,freight,M,b)
% 输入：way: 购买方案； fright:运费;  M: 每本书在每家店的价格; b：一共要购买几本书
   index = unique(way);  % 在哪些商店购买了商品，因为我们等下要计算运费
   money = sum(freight(index)); % 计算买书花费的运费
    % 计算总花费：刚刚计算出来的运费 + 五本书的售价
    for i = 1:b 
        money = money + M(way(i),i);  
    end
end
function way1 = gen_new_way(way0, s, b)
% way0：原来的买书方案，是一个1*b的向量，每一个元素都位于1-s之间
        index =  randi([1, b],1) ;  % 看哪一本书要更换书店购买
        way1 = way0;  % 将原来的方案赋值给way1
        way1(index) = randi([1, s],1);  % 将way1中的第index本书换一个书店购买   
end
%% 模拟退火解决书店买书问题  
tic
rng('shuffle')  % 控制随机数的生成，否则每次打开matlab得到的结果都一样

load book_data  % 这个文件一定要在当前文件夹下面
% 这个数据文件里面保存了两个矩阵：M是每本书在每家店的价格; freight表示每家店的运费
[s, b] = size(M);  % s是书店的数量，b是要购买的书的数量

%% 参数初始化
T0 = 1000;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 500;  % 最大迭代次数
Lk = 200;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数

%%  随机生成一个初始解
way0 = randi([1, s],1,b); % 在1-s这些整数中随机抽取一个1*b的向量，表示这b本书分别在哪家书店购买
money0 = calculate_money(way0,freight,M,b); % 调用我们自己写的calculate_money函数计算这个方案的花费

%% 定义一些保存中间过程的量，方便输出结果和画图
min_money = money0;     % 初始化找到的最佳的解对应的花费为money0
MONEY = zeros(maxgen,1); % 记录每一次外层循环结束后找到的min_money (方便画图）

%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  %  内循环，在每个温度下开始迭代
        way1 = gen_new_way(way0,s,b);  % 调用我们自己写的gen_new_way函数生成新的方案
        money1 = calculate_money(way1,freight,M,b); % 计算新方案的花费
        if money1 < money0    % 如果新方案的花费小于当前方案的花费
            way0 = way1; % 更新当前方案为新方案
            money0 = money1;
        else
            p = exp(-(money1 - money0)/T); % 根据Metropolis准则计算一个概率
            if rand(1) < p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                way0 = way1;
                money0 = money1;
            end
        end
        % 判断是否要更新找到的最佳的解
        if money0 < min_money  % 如果当前解更好，则对其进行更新
            min_money = money0;  % 更新最小的花费
            best_way = way0;  % 更新找到的最佳方案
        end
    end
    MONEY(iter) = min_money; % 保存本轮外循环结束后找到的最小花费
    T = alfa*T;   % 温度下降
end

disp('最佳的方案是：'); disp(mat2str(best_way))
disp('此时最优值是：'); disp(min_money)

%% 画出每次迭代后找到的最佳方案的图形
figure
plot(1:maxgen,MONEY,'b-');
xlabel('迭代次数');
ylabel('最小花费');
toc

问题4：背包问题 

 （如果有n件货物，那么可能性有2^n种）

Logo
开放原子开发者工作坊

开放原子开发者工作坊旨在鼓励更多人参与开源活动，与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动，如meetup、训练营等，主打技术交流，干货满满，真诚地邀请各位开发者共同参与！

更多推荐

按键精灵初步，使用winIO实现驱动级模拟按键

winIO是目前非常古老的技术，是用底层驱动来实现程序控制键盘按键，要使用winIO,需要用到两个文件WinIo64.dll，与WinIo64.sys。 也有32为的 WinIO32.dll, WinIO32.sys。基于现在的操作系统，都是用的32位的，且都是win7 以上的了，所以，本文用的是winIO64.dll，WinIO64.sys第一步:命令行运行： cmd-> bcd...

avatar 开放原子开发者工作坊

记一次Spark引擎执行Sql超时优化

前几天用spark引擎执行了一个较大的sql，涉及的表和数据量都不少，不同时间段执行了几次都超时，经过上网及分析，尝试解决了此问题，使用spark引擎测试几次大概都在半个小时左右完成，不再出现超时现象一、问题现象摘抄部分现场日志如下：2022-02-01 13:02:12 INFO 22/02/01 13:02:12 [dag-scheduler-event-loop] INFO DAGSched

avatar 开放原子开发者工作坊

HDMI设计1--HDMI 1.4b SPEC的阅读个人总结

工作需要设计HDMI1.4/2.0 IP,支持4K视频的需求。因此，近期对HDMI的相关SPEC进行了阅读，在这里进行总结，记录个人的知识成长的过程。本文将从一下7个方面,讲述自己对HDMI 1.4b spec相关内容的理解。因为，工作需求是只需要支持4K视频的接受和传输，不涉及音频，所以略去了HDMI中Audio部分的内容。1, Overview2, Signaling and Encoding

avatar 开放原子开发者工作坊