• 数学建模——模拟退火


    随机找一个初始点x_{i},在初始点附近找一个新点x_{j},对应函数值为f(x_{j}),如果新函数值大于旧函数值,则接受点x_{j},如果新函数值小于旧函数值,则以一定概率接受x_{j},新旧函数值相差越小,则概率越大。搜索前期搜索范围应该尽可能大,搜索后期倾向于局部搜索。概率可以定义为:

     Ct可以看做随时间变化,前期Ct较小,搜索的范围也就大,后期Ct大,精准搜索。

    流程:

     根据退火的理论,将Ct取为温度的倒数:

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

    下标t看成迭代的次数,为了保证搜索过程的彻底,在同一温度下(t),需要进行多次搜索。所以两层循环。当达到指定迭代次数、达到指定温度、找到的最优解连续M次不变化就可以退出循环过程。

    如何生成新的解,在不同问题中,方法不一样。

    求函数最大最小值

    产生新解方法:

     求解函数y = 11*sin(x) + 7*cos(5*x)在[-3,3]内的最大值

    1. tic
    2. clear; clc
    3. %% 绘制函数的图形
    4. x = -3:0.1:3;
    5. y = 11*sin(x) + 7*cos(5*x);
    6. figure
    7. plot(x,y,'b-')
    8. hold on % 不关闭图形,继续在上面画图
    9. %% 参数初始化
    10. narvs = 1; % 变量个数
    11. T0 = 100; % 初始温度
    12. T = T0; % 迭代中温度会发生改变,第一次迭代时温度就是T0
    13. maxgen = 200; % 最大迭代次数
    14. Lk = 100; % 每个温度下的迭代次数
    15. alfa = 0.95; % 温度衰减系数
    16. x_lb = -3; % x的下界
    17. x_ub = 3; % x的上界
    18. %% 随机生成一个初始解
    19. x0 = zeros(1,narvs);
    20. for i = 1: narvs
    21. x0(i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(1);
    22. end
    23. y0 = Obj_fun1(x0); % 计算当前解的函数值
    24. h = scatter(x0,y0,'*r'); % scatter是绘制二维散点图的函数(这里返回h是为了得到图形的句柄,未来我们对其位置进行更新)
    25. %% 定义一些保存中间过程的量,方便输出结果和画图
    26. max_y = y0; % 初始化找到的最佳的解对应的函数值为y0
    27. MAXY = zeros(maxgen,1); % 记录每一次外层循环结束后找到的max_y (方便画图)
    28. %% 模拟退火过程
    29. for iter = 1 : maxgen % 外循环, 我这里采用的是指定最大迭代次数
    30. for i = 1 : Lk % 内循环,在每个温度下开始迭代
    31. y = randn(1,narvs); % 生成1行narvs列的N(0,1)随机数
    32. z = y / sqrt(sum(y.^2)); % 根据新解的产生规则计算z
    33. x_new = x0 + z*T; % 根据新解的产生规则计算x_new的值
    34. % 如果这个新解的位置超出了定义域,就对其进行调整
    35. for j = 1: narvs
    36. if x_new(j) < x_lb(j)
    37. r = rand(1);
    38. x_new(j) = r*x_lb(j)+(1-r)*x0(j);
    39. elseif x_new(j) > x_ub(j)
    40. r = rand(1);
    41. x_new(j) = r*x_ub(j)+(1-r)*x0(j);
    42. end
    43. end
    44. x1 = x_new; % 将调整后的x_new赋值给新解x1
    45. y1 = Obj_fun1(x1); % 计算新解的函数值
    46. if y1 > y0 % 如果新解函数值大于当前解的函数值
    47. x0 = x1; % 更新当前解为新解
    48. y0 = y1;
    49. else
    50. p = exp(-(y0 - y1)/T); % 根据Metropolis准则计算一个概率
    51. if rand(1) < p % 生成一个随机数和这个概率比较,如果该随机数小于这个概率
    52. x0 = x1; % 更新当前解为新解
    53. y0 = y1;
    54. end
    55. end
    56. % 判断是否要更新找到的最佳的解
    57. if y0 > max_y % 如果当前解更好,则对其进行更新
    58. max_y = y0; % 更新最大的y
    59. best_x = x0; % 更新找到的最好的x
    60. end
    61. end
    62. MAXY(iter) = max_y; % 保存本轮外循环结束后找到的最大的y
    63. T = alfa*T; % 温度下降
    64. pause(0.01) % 暂停一段时间(单位:秒)后再接着画图
    65. h.XData = x0; % 更新散点图句柄的x轴的数据(此时解的位置在图上发生了变化)
    66. h.YData = Obj_fun1(x0); % 更新散点图句柄的y轴的数据(此时解的位置在图上发生了变化)
    67. end
    68. disp('最佳的位置是:'); disp(best_x)
    69. disp('此时最优值是:'); disp(max_y)
    70. pause(0.5)
    71. h.XData = []; h.YData = []; % 将原来的散点删除
    72. scatter(best_x,max_y,'*r'); % 在最大值处重新标上散点
    73. title(['模拟退火找到的最大值为', num2str(max_y)]) % 加上图的标题
    74. %% 画出每次迭代后找到的最大y的图形
    75. figure
    76. plot(1:maxgen,MAXY,'b-');
    77. xlabel('迭代次数');
    78. ylabel('y的值');
    79. toc

     Obj_fun1.m

    1. function y = Obj_fun1(x)
    2. y = 11*sin(x) + 7*cos(5*x);
    3. end

    旅行商问题

    产生新解方法:

  • 相关阅读:
    QT静态成员函数访问和操作UI对象
    CodeForces每日好题10.14
    基于SSM的学生(成绩)管理系统
    数组排序简介-计数排序(Counting Sort)
    EncodedResource类解读
    【QT】QRadioButton的使用(17)
    拿捏大厂面试,2022最新版的Java面试突击班手册
    UNCTF-日常训练-reverse-ezRust
    DockerFile的使用
    React 之 react-router-dom
  • 原文地址:https://blog.csdn.net/m0_51311105/article/details/126090506