随机找一个初始点,在初始点附近找一个新点,对应函数值为,如果新函数值大于旧函数值,则接受点,如果新函数值小于旧函数值,则以一定概率接受,新旧函数值相差越小,则概率越大。搜索前期搜索范围应该尽可能大,搜索后期倾向于局部搜索。概率可以定义为:
Ct可以看做随时间变化,前期Ct较小,搜索的范围也就大,后期Ct大,精准搜索。
流程:
根据退火的理论,将Ct取为温度的倒数:
温度一定时,越小,概率越大,即目标函数相差越小接受的可能性越大。一定时,温度越高,概率越大,即搜索前期温度较高,更有可能接受新解。
下标t看成迭代的次数,为了保证搜索过程的彻底,在同一温度下(t),需要进行多次搜索。所以两层循环。当达到指定迭代次数、达到指定温度、找到的最优解连续M次不变化就可以退出循环过程。
如何生成新的解,在不同问题中,方法不一样。
产生新解方法:
求解函数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
Obj_fun1.m
- function y = Obj_fun1(x)
- y = 11*sin(x) + 7*cos(5*x);
- end
产生新解方法: