傻瓜学建模-模拟退火算法

此博客根据清风老师的视频编写(https://www.bilibili.com/video/BV1hK41157JL?t=2141)
博客中会用到清风老师教案里的相关图片

1.问题的描述:

我们先看看以下几个问题

在这里插入图片描述

学过dp的同学会发现,上面四个问题中的后三个问题好像基本都能通过dp求解,没错确实可以,但这里我们看看另一种做法.

我们发现这些问题都可以转化为:求某个目标函数最值的问题,一般问题中我们求的都是最小值,如果碰到球最大值的情况也能加一个负号转化为求最小值.

2.初步的想法:

我们既然是求一个函数的最小值,那么我们可以遍历这个函数的定义域,对于一些考试中的数学题,我们往往只需要把枚举的元素差定义到0.01的数量级就足以解决问题,但放在实际问题中0.01这个数量级就先的不够精确,但当我们提升精确度时,算法的复杂度也就上去了,由此可见暴力的算法是不可行的.

我们先介绍一个概念:
启发式算法:在搜索最优解的过程中利用到原来搜索过程中得到的信息,且这个信息会优化我们的搜索过程.(讲的通俗一点就是用内存换时间)

下面我们介绍的爬山法就是一种启发式算法:
先看看如下的函数
在这里插入图片描述

1.在解空间里随机找一初始解

2.根据初始解的位置我们可以选择向左走还是向右走(走的距离越小越好,但会增加计算量),如果向左走函数值减小(假设我们要找最小值),那么我们就向左走.

扫描二维码关注公众号，回复： 12559018 查看本文章

3.重复2过程直到我们找到一个极小值

方法本身不难,但方法有很大的局限性:容易找到局部最优解.

要解决这个问题我们就需要引入模拟退火算法.

3.模拟退火算法:

3.1:模拟退火算法简介:

模拟退火算法来源于固体退火原理，是一种基于概率的算法，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。

3.2与爬山法的异同:
我们上面提出爬山法就是方便引出模拟退火,爬山法的主要缺点就在容易找到局部最优解,造成它的根本原因在于每次我们走的方向被函数值的大小限定的太死,导致函数值减小的方向我们根本不考虑.模拟退货算法就是在这一点上做出了改进.

3.3具体实现

在这里插入图片描述

如果我们拒绝了,不就是登山法了吗,那么肯定是不能拒绝的.
在这里插入图片描述

为了确定p,我们需要找到一个范围是0~1的函数
y=e^(-x)就是一个值域在0~1的函数,且随着x的增大而减小那么我们可以把函数变成
在这里插入图片描述
但这仍然不是最终形式,因为我们在搜索前期,为了使得搜索的范围更广,我们需要p在前期比较大,反之在后期比较小,那么p受到时间的影响.我们需要将函数在变形

在这里插入图片描述
Ct是一个受时间影响的变量

这个公式是根据一些物理知识推导而来

3.4代码:

代码中用到的方法(非初学者可跳过)

rand(x,y);//产生由随机生成的0~1的数组成的x*y的矩阵

scatter(x,y,'*r')//在x,y处绘制一个*r形式的点

title()//给函数设置标题

num2str()//数字转化成字符串

zeros(x,y);//生成一个x,y型的0矩阵

总代码:


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

%% 参数初始化
narvs = 1; % 变量个数
T0 = 100;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 300;  % 最大迭代次数
Lk = 200;  % 每个温度下的迭代次数
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的值
        %%此生成新解的方法参考matlab模拟退火
        % 如果这个新解的位置超出了定义域，就对其进行调整
        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

傻瓜学建模-模拟退火算法

猜你喜欢