模拟退火算法

引入

  模拟退火算法(SimulatedAnnealing Algorithm,简记为SAA或SA),它是基于金属退火机理而建立起来的一种最优化方法,它能够以随机搜索技术从概率意义上找出目标函数的全局最值点。

粒子群算法相比模拟退火有如下特点:

  1. 不仅能处理连续优化问题,还能很方便的处理组合优化问题,且编程简单易于实现,目标函数的收敛速度较快。
  2. 参数的选择至关重要,初始参数的合理选取是保证算法的全局收敛性和效率的关键,选择不当得到的结果可能会很差。
  3. 面对不同的问题,可以设计不同的退火过程,不同过程得到的准确度和效率可能有较大差别,因此需要一定的经验(我们可以学习相关领域的学者提出的论文)

爬山法找函数最大值

步骤

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

爬山法求解最大值

函数:

$y=11\sin \left( x \right) +7\cos \left( 5x \right)$

fun.m

function y = fun(x)
    y = 11 * sin(x) + 7 * cos(5 * x);
end

main.m

%% 
clear, close, clc
tic

%% 修改此处初始值即可
N = 20; % 一共运行N次
xMax = zeros(N, 1); % 每次的最大值所对应的x
upBound = 3; % 上界
downBound = -3; % 下界
n = 1e5; % 将可行域分为n份数
stepLength = (upBound - downBound) / n; % 步长
k = power(10, length(num2str(n)) - 3); % 间隔k次画一次图

%% 核心代码
for i = 1 : N
    tempX = unifrnd(downBound, upBound); % 在可行域内随机取一个值
    leftX = tempX; rightX = tempX; % 为左右两个变量赋初值
    x = linspace(downBound, upBound, n);
    y = fun(x);
    plot(x, y, 'b', 'LineWidth', 2)
%     text(1, 0.2, 'Interpreter', 'y=11sin\left( x \right) +7cos\left( 5x \right)')
    xlabel('x'), ylabel('y')
    axis([downBound, upBound -inf inf])
    grid on, hold on
    h = scatter(tempX, fun(tempX));
    h.LineWidth = 0.6;
    h.MarkerEdgeColor = 'r';
    h.MarkerFaceColor = 'r';
    for j = 1 : n
        leftX = leftX - stepLength;
        rightX = rightX + stepLength;
        if fun(leftX) > fun(rightX) && fun(leftX) > fun(tempX)
            tempX = leftX; % 向左移动
        elseif fun(leftX) < fun(rightX) && fun(tempX) < fun(rightX)
            tempX = rightX; % 向右移动
        else
            break;
        end
        if mod(j, k) == 0
            title(['当前最大值',num2str(fun(tempX))])
            pause(0.01)
            h.XData = tempX;
            h.YData = fun(tempX);
        end
    end
    hold off
    xMax(i) = tempX;
end
[maxY, a] = max(fun(xMax));
maxX = xMax(a);
h.XData = maxX;
h.YData = maxY;
title(['最大值为',num2str(maxY)])
fprintf('横坐标为:%.4f, 纵轴标为:%.4f\n', maxX, maxY)
toc

缺点
特别容易找到局部最优解。

模拟退火的介绍

  模拟退火算法的思想最初是由美国物理学家Metropolis(梅特罗波利斯)在1953年提出的,1983年,Kirkpatrick等将其应用于求解组合优化问题。其出发点是基于物理中固体物质的退火过程与一般的优化问题之间的相似性。

步骤

  1. 给定一个初始温度$T_0$,并随机生成一个初始解$x_0$,并计算相应的目标函数值$f\left( x_0 \right)$
  2. 令当前温度$?$等于冷却进度表中的下一个值$?_i$;(第一次迭代时$T=T_0$)
  3. 在当前解$x_i$的附近随机产生一个新解$x_j$,计算新解的目标函数值$f\left( x_j \right)$;(第一次迭代时$x_i=x_0$ )
  4. 如果$f\left( x_j \right) <f\left( x_i \right)$, 则接受新解$x_j$; 如果$f\left( x_j \right) >f\left( x_i \right)$, 则计算$\varDelta f=f\left( x_j \right) -f\left( x_i \right)$, 并计算$p=e^{-\varDelta f/T_i}$, 然后随机生成一个在区间$[0,1]$上服从均匀分布的随机数$?$,如果$r<p$, 则接受新解$x_j$;
  5. 在温度$T_i$下,将步骤(3)和(4)重复$?_i$次;
  6. 判断是否满足退出条件,如果满足则退出迭代,否则回到步骤(2)继续迭代。

接受新解的概率

  1. 温度一定时, $\varDelta f$越小,概率越大,即目标函数相差越小接受的可能性越大。
  2. $\varDelta f$一定时,温度越高,概率越大,即搜索前期温度较高时更有可能接受新解。

冷却进度表

注意
参数的设置很关键,模拟退火难就难在参数的选取上。
  1. 初始温度$T_0$最常取100
  2. 衰减函数最常用的是$T_{k+1}=\alpha \times T_k$,$\alpha$一般是0.5-0.99之间的常数。
  3. 止准则是指当温度下降到一个特别小的数时就退出搜索,即退出外循环。在后面我们会介绍其他退出外循环的规则。
  4. $Markov$链的长度$?_k$,这个是随机过程中马尔科夫链里面的概念,和上方第5点一个意思即在当前温度下需要内循环的次数。从经验来说,对于简单的情况可以令$?_k=100?$, 其中$?$为自变量的个数。

设置模拟退火的温度

设置模拟退火温度的初值$T_0$

  求解全局优化问题的随机搜索算法一般都采用大范围的粗略搜索与局部的精细搜索相结合的搜索策略。只有在初始的大范围搜索阶段找到全局最优解所在的区域,才能逐渐缩小搜索的范围,最终求出全局最优解。模拟退火算法是通过控制退火温度的初值$T_0$和其衰减变化过程来实现大范围的粗略搜索与局部的精细搜索。一般来说,只有足够大的$T_0$才能满足算法要求(但对不同的问题“足够大”的含义也不同,有的可能$T_0=100$ 就可以,有的则要1000。在问题规模较大时,过小的$T_0$往往导致算法难以跳出局部陷阱而达不到全局最优。但为了减少计算量, $T_0$不宜取得过大,而应与其他参数折中选取。

控制温度的衰减函数

  衰减函数可以有多种形式,一个常用的衰减函数是
$$
T_{k+1}=\alpha \times T_k\,\,k=0,1,2,\cdots
$$

  其中,?是一个常数,可以取为0.5-0.99,它的取值决定了降温的过程。$\alpha$较大会放慢温度衰减的速度,可能导致算法进程迭代次数的增加,从而使算法进程接受更多的变换,访问更多的邻域,搜索更大范围的解空间,返回更好的最终解,但求解的时间会有很大的增加,目前使用较多的?一般为0.95。可以多次尝试不同的参数,观察求解时间和求解结果。

新解的产生规则

模拟退火求解一元函数最大值

求解函数

fun.m

function y = fun(x)
    y = 11 * sin(x) + 7 * cos(5 * x);
end

main.m

%% 
clear, close, clc
tic

%% 修改此处初始值即可
var = 1; % 变量个数
T0 = 100; % 初始温度
T = T0; % 第一次迭代的温度T0
maxIter = 200; % 最大迭代次数
tmaxIter = 100; % 每个温度下的最大迭代次数
alpha = 0.95; % 温度衰减系数
upBound = [3]; % 上界
lowBound = [-3]; % 下界
n = 1e5; % 将可行域分为n份数

%% 绘制函数图形
x = linspace(lowBound, upBound, n);
y = fun(x);
subplot(1,2,1)
plot(x, y, 'b', 'LineWidth', 2)
xlabel('x'), ylabel('y')
hold on
grid on

%% 随机在可行域内生成一个初始解
x0 = zeros(var, 1);
for i = 1 : var
    x0(i) = unifrnd(lowBound(i), upBound(i));
end
y0 = fun(x0);
h = scatter(x0, y0);
h.LineWidth = 0.6;
h.MarkerEdgeColor = 'r';
h.MarkerFaceColor = 'r';

%%   中间变量&迭代后找到的最大y的图形
tempY = y0;
tempX = x0;
plotMaxY = tempY;
plotMaxX = 1;
subplot(1,2,2)
plotp = plot(plotMaxX, plotMaxY, 'b', 'LineWidth', 2);
xlabel('迭代次数'), ylabel('y的值')
axis([1 200 -inf 17.494])
grid on

%% 模拟退火(核心代码)
for i = 1 : maxIter
    for j = 1 : tmaxIter
        yi = randn(var, 1); % 生成1列var行的N(0,1)随机数
        zi = yi / sqrt(sum(yi .^ 2)); % 新解的产生规则计算z
        xNew = x0 + T * zi; % 新解的产生规则计算xNew的值
        for k = 1 : var
            if xNew(k) < lowBound(k)
                r = rand;
                xNew(k) = r * lowBound(k)+ (1 - r) * x0(k);
            elseif xNew(k) > upBound(k)
                r = rand;
                xNew(k) = r * upBound(k) + (1 - r) * x0(k);
            end
        end
        x1 = xNew;
        y1 = fun(x1);
        if y1 > y0
            x0 = x1;
            y0 = y1;
        else
            p = exp(-(y0 - y1)/ T);
            if rand < p
                x0 = x1;
                y0 = y1;
            end
        end
        if y0 > tempY
           tempY = y0;
           tempX = x0;
        end
    end
    plotMaxY(i) = tempY;
    plotMaxX(i) = i;
    T = alpha * T;
    title(['当前最大值为',num2str(tempY)])
    plotp.XData = plotMaxX;
    plotp.YData = plotMaxY;
    pause(0.01)
    h.XData = x0;
    h.YData = fun(x0);
end
title(['最大值为',num2str(tempY)])
fprintf('横坐标为:%.4f, 纵轴标为:%.4f\n', tempX, tempY)

toc

求解二元函数

$y=x{1}^{2}+x{2}^{2}-x_1x_2-10x_1-4x_2+60$
在$x_1,x_2\in \left[ -15,15 \right]$

%% 
clear, close, clc
tic

%% 修改此处初始值即可
var = 2; % 变量个数
T0 = 100; % 初始温度
T = T0; % 第一次迭代的温度T0
minIter = 200; % 最大迭代次数
tminIter = 100; % 每个温度下的最大迭代次数
alpha = 0.95; % 温度衰减系数
upBound = [15 15]; % 上界
lowBound = [-15 -15]; % 下界
n = 30; % 将可行域分为n份数

%% 绘制函数图形
x = linspace(lowBound(1), upBound(1), n);
[x1, x2] = meshgrid(x, x);
y = x1 .^ 2 + x2 .^ 2 - x1 .* x2 - 10 * x1 - 4 * x2 + 60;
subplot(1,2,1)
mesh(x1, x2, y)
xlabel('x1'), ylabel('x2'), zlabel('y')
axis vis3d
hold on

%% 随机在可行域内生成一个初始解
x0 = zeros(var, 1);
for i = 1 : var
    x0(i) = unifrnd(lowBound(i), upBound(i));
end
y0 = fun(x0);
h = scatter3(x0(1), x0(2), y0);
h.LineWidth = 0.6;
h.MarkerEdgeColor = 'r';
h.MarkerFaceColor = 'r';

%%   中间变量&迭代后找到的最大y的图形
tempY = y0;
tempX = x0;
plotMinY = tempY;
plotMinX = 1;
subplot(1,2,2)
plotp = plot(plotMinX, plotMinY, 'b', 'LineWidth', 2);
xlabel('迭代次数'), ylabel('y的值')
axis([0 200 7.95 inf])
grid on

%% 模拟退火(核心代码)
for i = 1 : minIter
    for j = 1 : tminIter
        yi = randn(var, 1); % 生成1列var行的N(0,1)随机数
        zi = yi / sqrt(sum(yi .^ 2)); % 新解的产生规则计算z
        xNew = x0 + T * zi; % 新解的产生规则计算xNew的值
        for k = 1 : var
            if xNew(k) < lowBound(k)
                r = rand;
                xNew(k) = r * lowBound(k)+ (1 - r) * x0(k);
            elseif xNew(k) > upBound(k)
                r = rand;
                xNew(k) = r * upBound(k) + (1 - r) * x0(k);
            end
        end
        x1 = xNew;
        y1 = fun(x1);
        if y1 < y0
            x0 = x1;
            y0 = y1;
        else
            p = exp(-(-y0 + y1)/ T);
            if rand < p
                x0 = x1;
                y0 = y1;
            end
        end
        if y0 < tempY
           tempY = y0;
           tempX = x0;
        end
    end
    plotMinY(i) = tempY;
    plotMinX(i) = i;
    T = alpha * T;
    title(['当前最小值为',num2str(tempY)])
    plotp.XData = plotMinX;
    plotp.YData = plotMinY;
    pause(0.01)
    h.XData = x0(1);
    h.YData = x0(2);
    h.ZData = fun(x0);
end
title(['最小值为',num2str(tempY)])
fprintf('x1为:%.13f\nx2为:%.13f\ny 为:%.13f\n', tempX(1), tempX(2), tempY)

toc

模拟退火解决旅行商问题

问题描述

问题描述

产生新解的方法

交换法

随机选择两个点,交换这两个点的位置

移位法

随机选择三个点,将前两个点之间的点移位到第三个点后

倒置法

随机选择两个点,将这两个点之间的顺序完全颠倒

swap.m

function [x, y] = swap(a, b)
    x = b;
    y = a;
end

newPath.m

function path = newPath(pathTmp)
    n = length(pathTmp);
    path = pathTmp;
    p1 = 1 / 3; % 使用交换法产生新路径的概率
    p2 = 2 / 3; % 使用移位法产生新路径的概率
    r = rand;  % 随机生成一个[0 1]间均匀分布的随机数
    indTmpOne = randi(n); % 生成1-n中的一个随机整数
    indTmpTwo = randi(n);
    indTmpThree = randi(n);
    if r < p1
        [path(indTmpOne), path(indTmpTwo)] = swap(pathTmp(indTmpOne), pathTmp(indTmpTwo));
    elseif r >= p1 && r < p2
        sortTmp = sort([indTmpOne, indTmpTwo, indTmpThree]);
        path = [pathTmp(1 : sortTmp(1) - 1), pathTmp(sortTmp(2) + 1 : sortTmp(3))...
            pathTmp(sortTmp(1) : sortTmp(2)), pathTmp(sortTmp(3) + 1 : end)];
    else
        if indTmpOne > indTmpTwo
            [indTmpOne, indTmpTwo] = swap(indTmpOne, indTmpTwo);
        end
        path = [pathTmp(1 : indTmpOne - 1), fliplr(pathTmp(indTmpOne : indTmpTwo))...
            pathTmp(indTmpTwo + 1 : end)];  % 矩阵的左右对称翻转 fliplr,上下对称翻转 flipud
    end
end

calcTspDist.m

function res = calcTspDist(path, d)
    n = length(path);
    res = 0;
    for i = 1 : n - 1
        res = res + d(path(i), path(i + 1)); % 按照这个序列不断的更新走过的路程这个值
    end
    res = res + d(path(1), path(n)); % 加上从最后一个城市返回到最开始那个城市的距离
end

main.m

%% 
clear, close, clc
rng('shuffle')  % 控制随机数的生成,否则每次打开matlab得到的结果都一样
tic

%% 修改此处初始值即可
cityCoord = [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]; % 38个城市坐标
n = length(cityCoord); % 城市的数目
T0 = 1000; % 初始温度
T = T0; % 第一次迭代的温度T0
maxIter = 1000; % 最大迭代次数
tmaxIter = 500; % 每个温度下的最大迭代次数
alpha = 0.95; % 温度衰减系数

%% 两个城市的距离矩阵
dist = zeros(n);
for i = 2 : n  
    for j = 1 : i - 1
        dist(i,j) = sqrt((cityCoord(i, 1) - cityCoord(j, 1))^2 + (cityCoord(i, 2) - cityCoord(j, 2))^2);   % 计算城市i和j的距离
    end
end
dist = dist + dist'; % 生成距离矩阵的对称的一面

%% 随机在可行域内生成一个初始解
pathTmp = randperm(n); % 生成一个1到n的随机打乱的序列作为初始的路径
resTmp = calcTspDist(pathTmp, dist);  % 调用我们自己写的calcTspDist函数计算当前路径的距离

%%   中间变量&迭代后找到的最小y的图形
plotMaxY = resTmp;
plotMaxX = 1;
minRes = resTmp;
% 画图可以去掉注释
% subplot(2,1,2)
% plotp = plot(plotMaxX, plotMaxY, 'b', 'LineWidth', 2);
% xlabel('迭代次数'), ylabel('y的值')
% axis([0 maxIter 6000 inf])
% grid on

%% 模拟退火(核心代码)
count = 0;
for i = 1 : maxIter
    for j = 1 : tmaxIter
        path = newPath(pathTmp); % 调用我们自己写的newPath函数生成新的路径
        res = calcTspDist(path, dist); % 计算新路径的距离
        % 如果新解距离短,则直接把新解的值赋值给原来的解
        if res < resTmp
            resTmp = res; % 更新当前路径为新路径
            pathTmp = path;
        else
            p = exp(-(res - resTmp) / T); % 根据Metropolis准则计算一个概率
            if rand < p
                resTmp = res;
                pathTmp = path; % 更新当前路径为新路径
            end
        end
        if resTmp < minRes
            minRes = resTmp;
            bestPath = pathTmp;
        end
    end
    % 如果count次最优值还没有发生变化提前退出循环
    if resTmp == minRes
        count = count + 1;
        if count == 100
           break; 
        end
    end
    plotMinY(i) = minRes;
    plotMinX(i) = i;
    T = alpha * T;
    % 画图可以去掉注释
%     plotp.XData = plotMinX;
%     plotp.YData = plotMinY;
%     subplot(2,1,1)
%     scatter(cityCoord(:, 1), cityCoord(:, 2));
%     grid on
%     hold on
%     for iteri = 1 : n - 1
%         plot([cityCoord(bestPath(iteri), 1), cityCoord(bestPath(iteri + 1), 1)]...
%             ,[cityCoord(bestPath(iteri), 2), cityCoord(bestPath(iteri + 1), 2)], '-b')
%         hold on
%     end
%     plot([cityCoord(bestPath(end), 1), cityCoord(bestPath(1), 1)]...
%             ,[cityCoord(bestPath(end), 2), cityCoord(bestPath(1), 2)], '-r')
%     hold off
%     pause(1e-8)
end

disp('最佳的方案是:'); disp(mat2str(pathTmp))
disp('此时最优值是:'); disp(resTmp)

toc

模拟退火解决书店买书问题(0-1规划)

bookData.mat

calcMoney.m

function moneyTmp = calcMoney(wayTmp, freight, M, buyBookNum)
    ind = unique(wayTmp); % 相同的书店路费只算一次
    moneyTmp = sum(freight(ind)); % 买书花的路费
    % 计算总费用
    for i = 1 : buyBookNum
        moneyTmp = moneyTmp + M(wayTmp(i), i);
    end
end

newWay.m

function way = newWay(wayTmp, bookStoreNum, buyBookNum)
    ind = randi([1, buyBookNum], 1); % 看看哪本书要更换书店购买
    way = wayTmp;
    way(ind) = randi([1, bookStoreNum], 1);
end

main.m

%% 
clear, close, clc
rng('shuffle')  % 控制随机数的生成,否则每次打开matlab得到的结果都一样
load book_Data % 导入当前目录中的mat数据
tic

%% 修改此处初始值即可
[bookStoreNum, buyBookNum] = size(M); % bookStoreNum是书店的数量,buyBookNum是要购买的书的数量
T0 = 1000; % 初始温度
T = T0; % 第一次迭代的温度T0
maxIter = 500; % 最大迭代次数
tmaxIter = 300; % 每个温度下的最大迭代次数
alpha = 0.95; % 温度衰减系数

%% 随机在可行域内生成一个初始解
wayTmp = randi([1, bookStoreNum], 1, buyBookNum); % 在1到bookStoreNum这些整数中随机抽取一个buyBookNum*1的向量,表示这buyBookNum本书分别在哪家书店购买
moneyTmp = calcMoney(wayTmp, freight, M, buyBookNum); % 调用我们自己写的ccalcMoney函数计算这个方案的花费

%%   中间变量&迭代后找到的最小y的图形
plotMaxY = moneyTmp;
plotMaxX = 1;
minMoney = moneyTmp;
% 画图可以去掉注释
plotp = plot(plotMaxX, plotMaxY, 'b', 'LineWidth', 2);
xlabel('迭代次数'), ylabel('y的值')
axis([0 200 465 inf])
grid on

%% 模拟退火(核心代码)
count = 0;
for i = 1 : maxIter
    for j = 1 : tmaxIter
        way = newWay(wayTmp, bookStoreNum, buyBookNum); % 调用我们自己写的newWay函数生成新的路径
        money = calcMoney(way, freight, M, buyBookNum); % 计算新方案的花费
        % 如果新解距离短,则直接把新解的值赋值给原来的解
        if money < moneyTmp
            moneyTmp = money; % 更新当前路径为新路径
            wayTmp = way;
        else
            p = exp(-(money - moneyTmp) / T); % 根据Metropolis准则计算一个概率
            if rand < p
                moneyTmp = money;
                wayTmp = way; % 更新当前方案
            end
        end
        if moneyTmp < minMoney
           minMoney = moneyTmp;
           bestWay = wayTmp;
        end
    end
    % 如果count次最优值还没有发生变化提前退出循环
    if moneyTmp == minMoney
       count = count + 1;
       if count == 200
           break;
       end
    end
    plotMinY(i) = minMoney;
    plotMinX(i) = i;
    T = alpha * T;
    title(['当前最小值为',num2str(minMoney)])
    % 画图可以去掉注释
    plotp.XData = plotMinX;
    plotp.YData = plotMinY;
    pause(0.001)
end
title(['最优值为',num2str(minMoney)])
disp('最佳的方案是:'); disp(num2str(bestWay))
disp('此时最优值是:'); disp(minMoney)

toc

最佳的方案是:
6 14 6 6 6 14 11 6 14 6 14 4 11 14 14 6 4 4 14 11
此时最优值是:
466


练习

练习一

求解$y=21.5+x_1\sin \left( 4\pi x_1 \right) +x_2\sin \left( 20\pi x_2 \right),x_1\in \left[ -3,12.1 \right] ,x_2\in \left[ 4.1,5.8 \right]$

fun.m

function y = fun(x)
    y = 21.5 + x(1)*sin(4*pi*x(1)) + x(2)*sin(20*pi*x(2));
end

main.m

%%
clear, close, clc
tic

%% 修改此处初始值即可
var = 2; % 变量个数
T0 = 100; % 初始温度
T = T0; % 第一次迭代的温度T0
minIter = 1000; % 最大迭代次数
tminIter = 200; % 每个温度下的最大迭代次数
alpha = 0.95; % 温度衰减系数
upBound = [12.1 5.8]; % 上界
lowBound = [-3 4.1]; % 下界
% n = 30; % 将可行域分为n份数

%% 绘制函数图形
x1 = -3:0.01:12.1;
x2 = 4.1:0.01:5.8;
[x1, x2] = meshgrid(x1, x2);
y = 21.5 * x1 .* sin(4 * pi * x1) + x2 .* sin(20 * pi * x2);
subplot(1,2,1)
mesh(x1, x2, y)
xlabel('x1'), ylabel('x2'), zlabel('y')
axis vis3d
hold on

%% 随机在可行域内生成一个初始解
x0 = zeros(var, 1);
for i = 1 : var
    x0(i) = unifrnd(lowBound(i), upBound(i));
end
y0 = fun(x0);
h = scatter3(x0(1), x0(2), y0);
h.LineWidth = 0.6;
h.MarkerEdgeColor = 'r';
h.MarkerFaceColor = 'r';

%%   中间变量&迭代后找到的最大y的图形
tempY = y0;
tempX = x0;
plotMinY = tempY;
plotMinX = 1;
subplot(1,2,2)
plotp = plot(plotMinX, plotMinY, 'b', 'LineWidth', 2);
xlabel('迭代次数'), ylabel('y的值')
% axis([0 200 7.95 inf])
grid on

%% 模拟退火(核心代码)
for i = 1 : minIter
    for j = 1 : tminIter
        yi = randn(var, 1); % 生成1列var行的N(0,1)随机数
        zi = yi / sqrt(sum(yi .^ 2)); % 新解的产生规则计算z
        xNew = x0 + T * zi; % 新解的产生规则计算xNew的值
        for k = 1 : var
            if xNew(k) < lowBound(k)
                r = rand;
                xNew(k) = r * lowBound(k)+ (1 - r) * x0(k);
            elseif xNew(k) > upBound(k)
                r = rand;
                xNew(k) = r * upBound(k) + (1 - r) * x0(k);
            end
        end
        x1 = xNew;
        y1 = fun(x1);
        if y1 > y0
            x0 = x1;
            y0 = y1;
        else
            p = exp(-(y0 - y1)/ T);
            if rand < p
                x0 = x1;
                y0 = y1;
            end
        end
        if y0 > tempY
           tempY = y0;
           tempX = x0;
        end
    end
    plotMinY(i) = tempY;
    plotMinX(i) = i;
    T = alpha * T;
    title(['当前最大值为',num2str(tempY)])
    plotp.XData = plotMinX;
    plotp.YData = plotMinY;
    pause(0.01)
    h.XData = x0(1);
    h.YData = x0(2);
    h.ZData = fun(x0);
end
title(['最大值为',num2str(tempY)])
fprintf('x1为:%.13f\nx2为:%.13f\ny 为:%.13f\n', tempX(1), tempX(2), tempY)

toc


练习二

书店买书问题中:加入下面这个规则:“如果在同一个书店的消费额不低于88元,那么这个书店包邮,即在这个书店买书不用出运费了”,请你修改代码对该问题重新求解。

只需要重写calcMoney.m即可

注意
要记得修改上方main.m中调用的calcMoney函数参数个数哦

calcMoney.m

function moneyTmp = calcMoney(wayTmp, freight, M, bookStoreNum, buyBookNum)
    % 计算书费
    costBook = zeros(bookStoreNum, 1);
    for i = 1 : buyBookNum
        costBook(wayTmp(i)) = costBook(wayTmp(i)) + M(wayTmp(i), i);
    end
    moneyTmp = sum(costBook);
    % 计算总价
    ind = unique(wayTmp);
    for i = ind
        if costBook(i) < 88
            moneyTmp = moneyTmp + freight(i);
        end
    end
end


练习三

练习四

练习五

40 名同学需要分到 10 个寝室,我们通过调查问卷得到了这 40 名同学的兴趣爱好,完整的数据在下方excel里面

数据

问怎么划分寝室能够保证:不同寝室间,同寝室同学的共同爱好数量分布的尽可能均匀。

思路

理解题意

分布的尽可能均匀:方差(标准差)尽可能小

把每个寝室看成一个整体,先计算这个寝室中 4 名同学的共同爱好数量。例如,记符号(x,y)表示 x 和 y 两名同学的共同爱好数量,如果一个寝室有 ABCD 四名同学,那么这个寝室的共同爱好数量应该等于:(A,B)+(A,C)+(A,D)+(B,C)+(B,D)+(C,D),那么,我们分别计算出 10 个寝室中每个寝室的四名同学的共同爱好数量,记录到 tmp 向量里面,然后再计算这个向量的方差或者标准差即可,我们的目标函数就是让 tmp 的方差或者标准差尽可能小。

数据预处理

  1. 删除兴趣爱好中的其他这一项,因为太宽泛了
  2. 每个同学的兴趣看成一个集合(集合中的元素具有互异性),利用集合之间的“交集∩”运算可以计算出任意两名同学的共同兴趣爱好的集合,计算这个集合里面元素个数即可得到任意两名同学的共同兴趣爱好数。(这里面涉及到基础的字符串处理编程)
import numpy as np 
import pandas as pd 
df = pd.read_csv('./data.csv', encoding='gbk') 
data = df.iloc[:,1].tolist() 
data = [set(i.replace('其他','').strip('、').split('、')) for i in data] 
result = [len(i & j)  for i in data for j in data] 
result = np.array(result).reshape(len(data),-1) 
np.savetxt('./mydata.csv', result, delimiter=',') 

处理好的mycode.csv

最优的划分方案

我们这里可以考虑使用模拟退火算法40 名同学分别用 1-40 替代

  1. 定义初始解:1-40 的一个随机排列
  2. 目标函数:10 个寝室的共同爱好数量的标准差尽可能小
  3. 产生新解的方式:类似于旅行商问题(TSP)

改进

  假设方案 1 中,每个寝室的共同爱好数量都等于 10;方案 2 中每个寝室的共同爱好数量都等于 12,那么方案 1 和方案 2 计算得到的目标函数都等于 0(一组数中各元素都相同意味着没有波动,标准差为 0),但是很明显我们更加偏好方案 2。因此,在保证各寝室间同学们的共同爱好数量分布的尽可能均匀的同时,寝室中的共同爱好的数量也越多越好。(多目标问题)
定义新的目标函数:(std(tmp) + 1) /max(tmp) 这个目标函数越小越好

newPath

getHobby.m

function hobby = getHobby(pathTmp, mydata)
    hobby = 0;
    for i = 2 : 4
            for j = 1 : i - 1
                if i ~= j
                    hobby = hobby + mydata(pathTmp(i), pathTmp(j));
                end
            end
    end
end

calcDist.m

function [ans, tmp] = calcDist(pathTmp, mydata)
    n = length(pathTmp); % 学生个数
    nums = n / 4; % 寝室个数
    pathTmp = reshape(pathTmp, 4, nums); % 每四个人分配到一个寝室
    tmp = zeros(1, nums); % tmp用来保存每个寝室的共同爱好数
    for i = 1 : nums
        tmp(i) = getHobby(pathTmp(:, i), mydata); % 调用子函数单独计算每个寝室四名同学的共同爱好数
    end
%     ans = std(tmp);   % 计算标准差
    ans = (std(tmp) + 1)  / max(tmp);
%     ans = (std(tmp) + 1)  / mean(tmp);
end

main.m

%% 
clear, close, clc
rng('shuffle')  % 控制随机数的生成,否则每次打开matlab得到的结果都一样
mydata = csvread('mydata.csv');
tic

%% 修改此处初始值即可
n = length(mydata); % 城市的数目
T0 = 1000; % 初始温度
T = T0; % 第一次迭代的温度T0
maxIter = 1000; % 最大迭代次数
tmaxIter = 500; % 每个温度下的最大迭代次数
alpha = 0.95; % 温度衰减系数

%% 随机在可行域内生成一个初始解
pathTmp = randperm(n);
resTmp = calcDist(pathTmp, mydata);  % 计算这一组解对应的标准差

%%   中间变量&迭代后找到的最小y的图形
plotMaxY = resTmp;
plotMaxX = 1;
minRes = resTmp;
% 画图可以去掉注释
% plotp = plot(plotMaxX, plotMaxY, 'b', 'LineWidth', 2);
% xlabel('迭代次数'), ylabel('y的值')
% axis([0 maxIter 6000 inf])
% grid on

%% 模拟退火(核心代码)
count = 0;
for i = 1 : maxIter
    for j = 1 : tmaxIter
        path = newPath(pathTmp);
        res = calcDist(path, mydata); 
        if res < resTmp
            resTmp = res; 
            pathTmp = path;
        else
            p = exp(-(res - resTmp) / T); % 根据Metropolis准则计算一个概率
            if rand < p
                resTmp = res;
                pathTmp = path; 
            end
        end
        if resTmp < minRes
            minRes = resTmp;
            bestPath = pathTmp;
        end
    end
    % 如果count次最优值还没有发生变化提前退出循环
    if resTmp == minRes
        count = count + 1;
        if count == 100
           break; 
        end
    end
    plotMinY(i) = minRes;
    plotMinX(i) = i;
    T = alpha * T;
    % 画图可以去掉注释
%     plotp.XData = plotMinX;
%     plotp.YData = plotMinY;
%     pause(1e-8)
end

disp('最佳的方案是:'); disp(num2str(reshape(bestPath, 4, n / 4)))
disp('此时最优值是:'); disp(minRes)

toc

文章作者: Mr.Zhang
文章标题: 模拟退火算法
文章链接: https://www.codewing.cn/index.php/329/simulated-annealing-algorithm/
版权声明: 本博客所有文章除特别声明外,均采用CC BY-NC-SA 4.0 许可协议。转载请注明来自Ying's Blog!
暂无评论

发送评论 编辑评论


				
|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
(╯‵□′)╯︵┴─┴
 ̄﹃ ̄
(/ω\)
∠( ᐛ 」∠)_
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ`)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ( ̄∇ ̄o)
ヾ(´・ ・`。)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò。)
Σ(っ °Д °;)っ
( ,,´・ω・)ノ"(´っω・`。)
╮(╯▽╰)╭
o(*////▽////*)q
>﹏<
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
颜文字
Emoji
小恐龙
花!
上一篇
下一篇