2024国赛数学建模-模仿火算法（MATLAB 实现)

• 模仿退火算法
  

1.1 算法原理 模仿退火算法的基本思想是从一给定解开始 ,从邻域 中随机产生另一个解 ,接受 Metropolis准则允许目的函数在 有限范围内变坏 ,它由一控制参数 t决定 ,其作用类似于物 理过程中的温度 T,对于控制参数的每一取值 ,算法持续进 行“产生 —判定 —接受或舍去 ”的迭代过程 ,对应着固体在 某一恒定温度下的趋于热平衡的过程 ,当控制参数逐渐减 小并趋于 0时 ,系统越来越趋于平衡态 ,末了系统状态对应于优化问题的全局最优解 ,该过程也称为冷却过程 ,由于固 体退火必须痴钝降温 ,才能使固体在每一温度下都达到热 平衡 ,最终趋于平衡状态 ,因此控制参数 t经痴钝衰减 ,才 能确保模仿退火算法最终优化问题的团体最优解。

• 2 算法详细步骤
  

(1)给定模子每一个参数变革范围 ,在这个范围内随 机选择一个初始模子 m0 ,并计算相应的目的函数值 E (m0 )。
(2)对当前模子进行扰动产生一个新模子 m,计算相应 的目的函数值 E (m) ,得到 ΔE = E (m) - E (m0 )。
(3)若 ΔE < 0,则新模子被接受 ;若 ΔE > 0,则新模子 m 按概率 P = exp ( -ΔE/T)进行接受 , T为温度。当模子被接 受时 ,置 m0 =m, E (m0 ) = E (m)。
(4)在温度 T下 ,重复一定次数的扰动和接受过程 ,即 重复步骤 (2)、(3)。
(5)痴钝低落温度 T。 
(6)重复步骤 (2)、(5) ,直至收敛条件满足为止。
算法的实质分两次循环 ,随机扰动产生新模子并计算 目的函数值 (或称能量 )的变革 ,决定是否被接受。由于算 法初始温度筹划在高温条件 ,这使得 E增大的模子大概被 接受 ,因而能舍去局部极小值 ,通过痴钝地低落温度 ,算法 最终能收敛到全局最优点。

实行用例：用模仿退火算法办理如下 10 个城市的 TSP 问题（Traveling Salesman Problem，旅行商问题），由威廉哈密顿爵士和英国数学家克克曼T.P.Kirkman于19世纪初提出。 问题形貌如下： 有多少个城市，任何两个城市之间的间隔都是确定的，现要求一旅行商从某城市出发必须经过每一个城市且只在一个城市逗留一次，末了回到出发的城市，问怎样事先确定一条最短的线路已包管其旅行的费用最少？)，该问题最优解为 f_opt = 2.691。

编程实现 
用 MATLAB 实现模仿退火算法时，共编制了 5 个 m 文件，分别如下

• swap.m
  

1. function [ newpath , position ] = swap( oldpath , number )% 对 oldpath 进 行 互 换 操 作% number 为 产 生 的 新 路 径 的 个 数% position 为 对 应 newpath 互 换 的 位 置m = length( oldpath ) ; % 城 市 的 个 数newpath = zeros( number , m ) ;position = sort( randi( m , number , 2 ) , 2 ); % 随 机 产 生 交 换 的 位 置for i = 1 : number newpath( i , : ) = oldpath ;% 交 换 路 径 中 选 中 的 城 市 newpath( i , position( i , 1 ) ) = oldpath( position( i , 2 ) ) ; newpath( i , position( i , 2 ) ) = oldpath( position( i , 1 ) ) ;end

复制代码

2.pathfare.m

1. function [ objval ] = pathfare( fare , path )% 计 算 路 径 path 的 代 价 objval% path 为 1 到 n 的 排 列 ，代 表 城 市 的 访 问 顺 序 ；% fare 为 代 价 矩 阵 ， 且 为 方 阵 。[ m , n ] = size( path ) ;objval = zeros( 1 , m ) ;for i = 1 : m for j = 2 : n  objval( i ) = objval( i ) + fare( path( i , j - 1 ) , path( i , j ) ) ; end objval( i ) = objval( i ) + fare( path( i , n ) , path( i , 1 ) ) ;end​

复制代码

3、distance.m

1. function [ fare ] = distance( coord )% 根 据 各 城 市 的 距 离 坐 标 求 相 互 之 间 的 距 离% fare 为 各 城 市 的 距 离 ， coord 为 各 城 市 的 坐 标[ ~ , m ] = size( coord ) ; % m 为 城 市 的 个 数fare = zeros( m ) ;for i = 1 : m % 外 层 为 行 for j = i : m % 内 层 为 列 fare( i , j ) = ... ( sum( ( coord( : , i ) - coord( : , j ) ) .^ 2 ) ) ^ 0.5 ; fare( j , i ) = fare( i , j ) ; % 距 离 矩 阵 对 称 endend

复制代码

4、myplot.m

1. function [ ] = myplot( path , coord , pathfar )% 做 出 路 径 的 图 形% path 为 要 做 图 的 路 径 ，coord 为 各 个 城 市 的 坐 标% pathfar 为 路 径 path 对 应 的 费 用len = length( path ) ;clf ;hold on ;title( [ '近似最短路径如下，费用为' , num2str( pathfar ) ] ) ;plot( coord( 1 , : ) , coord( 2 , : ) , 'ok');pause( 0.4 ) ;for ii = 2 : len plot( coord( 1 , path( [ ii - 1 , ii ] ) ) , coord( 2 , path( [ ii - 1 , ii ] ) ) , '-b'); x = sum( coord( 1 , path( [ ii - 1 , ii ] ) ) ) / 2 ; y = sum( coord( 2 , path( [ ii - 1 , ii ] ) ) ) / 2 ; text( x , y , [ '(' , num2str( ii - 1 ) , ')' ] ) ; pause( 0.4 ) ;endplot( coord( 1 , path( [ 1 , len ] ) ) , coord( 2 , path( [ 1 , len ] ) ) , '-b' ) ;x = sum( coord( 1 , path( [ 1 , len ] ) ) ) / 2 ;y = sum( coord( 2 , path( [ 1 , len ] ) ) ) / 2 ;text( x , y , [ '(' , num2str( len ) , ')' ] ) ;pause( 0.4 ) ;hold off ;​

复制代码

5、mySAA.m

1. % 模 拟 退 火 算 法 ( Simulated Annealing Algorithm ) MATLAB 程 序clear ;% 程 序 参 数 设 定Coord = ... % 城 市 的 坐 标 Coordinates [ 0.6683 0.6195 0.4 0.2439 0.1707 0.2293 0.5171 0.8732 0.6878 0.8488 ; ... 0.2536 0.2634 0.4439 0.1463 0.2293 0.761 0.9414 0.6536 0.5219 0.3609 ] ;t0 = 1 ; % 初 温 t0iLk = 20 ; % 内 循 环 最 大 迭 代 次 数 iLkoLk = 50 ; % 外 循 环 最 大 迭 代 次 数 oLklam = 0.95 ; % λ lambdaistd = 0.001 ; % 若 内 循 环 函 数 值 方 差 小 于 istd 则 停 止ostd = 0.001 ; % 若 外 循 环 函 数 值 方 差 小 于 ostd 则 停 止ilen = 5 ; % 内 循 环 保 存 的 目 标 函 数 值 个 数olen = 5 ; % 外 循 环 保 存 的 目 标 函 数 值 个 数% 程 序 主 体m = length( Coord ) ; % 城 市 的 个 数 m fare = distance( Coord ) ; % 路 径 费 用 farepath = 1 : m ; % 初 始 路 径 pathpathfar = pathfare( fare , path ) ; % 路 径 费 用 path fareores = zeros( 1 , olen ) ; % 外 循 环 保 存 的 目 标 函 数 值e0 = pathfar ; % 能 量 初 值 e0t = t0 ; % 温 度 tfor out = 1 : oLk % 外 循 环 模 拟 退 火 过 程 ires = zeros( 1 , ilen ) ; % 内 循 环 保 存 的 目 标 函 数 值 for in = 1 : iLk % 内 循 环 模 拟 热 平 衡 过 程 [ newpath , ~ ] = swap( path , 1 ) ; % 产 生 新 状 态 e1 = pathfare( fare , newpath ) ; % 新 状 态 能 量 % Metropolis 抽 样 稳 定 准 则 r = min( 1 , exp( - ( e1 - e0 ) / t ) ) ; if rand < r path = newpath ; % 更 新 最 佳 状 态 e0 = e1 ; end ires = [ ires( 2 : end ) e0 ] ; % 保 存 新 状 态 能 量 % 内 循 环 终 止 准 则 ：连 续 ilen 个 状 态 能 量 波 动 小 于 istd if std( ires , 1 ) < istd break ; end end ores = [ ores( 2 : end ) e0 ] ; % 保 存 新 状 态 能 量% 外 循 环 终 止 准 则 ：连 续 olen 个 状 态 能 量 波 动 小 于 ostd if std( ores , 1 ) < ostd break ; end t = lam * t ; endpathfar = e0 ;% 输 入 结 果fprintf( '近似最优路径为：\n ' )%disp( char( [ path , path(1) ] + 64 ) ) ;disp(path)fprintf( '近似最优路径费用\tpathfare=' ) ;disp( pathfar ) ;myplot( path , Coord , pathfar ) ;​

复制代码

我试着运行了几次（只是改变了一下初温，也可以更改一下其他参数），发现初始温度 t0=1 时程序的末了结果与最优解差距小的概率比较大。 希望对大家有用！！​


免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。