matlab实现模拟退火算法

                              模拟退火算法

在这里插入图片描述
物理退火原理概述:

对于固体,微观粒子的无序度较大,为了使粒子排列更为有序,就需要对其加热,加高到一定温度的过程中,固体的内能增大,粒子热运动更加剧烈。然后逐渐降温,在每一温度梯度中,粒子逐渐变得有序,随着温度的下降,粒子运动不再那么剧烈,最后变为趋于有序的状态。

模拟物理的退火过程,创造了一种寻找全局最优解的近似解的一种新算法-------模拟退火算法。

我们先来看一个具体的问题:TSP(旅行商问题)

有n个城市,一个商人在A城市,为了能到每一个城市宣传自己的商品,他要寻找一个最短的路径,并且最后回到A城市。

问题分析:在A城市,下一步有(n-1)种选择,假设到达B城市,那么下一步还有(n-2)种选择…最后回到A城市,总共有(n-1)!种路径可供选择,那么哪一种最短呢?最直接的方法就是把每一种路径的距离都算一下,再做个比较取最小值不就好了,可是当n取比较大的数值(比如北京,上海,杭州,承德.等34个城市,33!=8乘10^36,也许你会觉得这个数并不‘大’,那么我们来具体算一下,假设你每秒钟能够计算出一种方案的结果,并且一天24小时一直以这种速度计算,那么在有生之年,假设还能活100年,能够计算完的方案数是:100乘365乘24乘60乘60=3乘10的9次方,如果计算机每秒能完成24个城市的枚举,在这一百年里,也只不过算出了7乘10的10次方个)时,显然此方案不可行。那么又该如何做呢?目前还没有一种算法能够精确求出最优解,只能通过模拟退火算法求出最优解的近似解。

模拟退火算法

Metropolis准则

为了能够找到全局的最优解,需要在解空间中遍历搜索,而通常情况下当在已知解的左右找到的新解不如已知解好的话,会认为已知的解就是最优的解,但比较的空间是在一个很小的邻域里面,这个不一定是全局的最优解,所以为了跳出原来的那个小邻域,就要以一定的概率接受这个不如已知解好的新解,这就是Metropolis准则。这个概率是(见图二下边),K是玻尔兹曼常数,K=1.38乘10的-32次方。可以看出,当新解大于已知解时,随着温度的逐渐降低,这个概率逐渐下降趋于零,也就是温度越低,越不太可能接受非最优解。

算法描述:
初始温度设定,初始解的设定,恒温迭代次数的确定,温度最低值的确定,
恒温搜索新解的结果要记录下来,当完成内层迭代后,需要计算一下保存的这几个结果的方差,判断是否趋于稳定(局部最优),是的话跳出循环,执行降温,不是的话继续本层温度的迭代。实际上内层与外层都应用到了Metropolis准测接受非最优解,只不过接受的概率不一样,温度高的接受的概率大,温度低的接受的概率小。

% 模 拟 退 火 算 法 ( Simulated Annealing Algorithm ) MATLAB 程 序
%模拟火算法(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 ; % 初 温 t0
iLk = 20 ; % 内 循 环 最 大 迭 代 次 数 iLk
oLk = 50 ; % 外 循 环 最 大 迭 代 次 数 oLk
lam = 0.95 ; % λ lambda
istd = 0.001 ; % 若 内 循 环 函 数 值 方 差 小 于 istd 则 停 止
ostd = 0.001 ; % 若 外 循 环 函 数 值 方 差 小 于 ostd 则 停 止
ilen = 5 ; % 内 循 环 保 存 的 目 标 函 数 值 个 数
olen = 5 ; % 外 循 环 保 存 的 目 标 函 数 值 个 数
% 程 序 主 体
m = length( Coord ) ; % 城 市 的 个 数 m
fare = distance( Coord ) ; 
path = 1 : m ; % 初 始 路 径 path
pathfar = pathfare( fare , path ) ; % 路 径 费 用 path fare
ores = zeros( 1 , olen ) ; % 外 循 环 保 存 的 目 标 函 数 值
e0 = pathfar ; % 能 量 初 值 e0
t = t0 ; % 温 度 t
for 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 ;
%模拟火算法(MATLAB 实现)
end
pathfar = e0 ;
% 输 入 结 果
fprintf( '近似最优路径为:\n ' )
%disp( char( [ path , path(1) ] + 64 ) ) ;
disp(path)
fprintf( '近似最优路径费用\tpathfare=' ) ;
disp( pathfar ) ;
myplot( path , Coord , pathfar)

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 ) ; % 距 离 矩 阵 对 称主对角线为零
end
end
end

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

end


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( 3 ) ; 
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 ) ; 
end 
plot( 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 ; 
end


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
end


在这里插入图片描述

  • 8
    点赞
  • 29
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
以下是一个简单的 Matlab 实现模拟退火算法代码: ``` function [x,fval] = simulanneal(fun,x0,options) % fun: 目标函数 % x0: 初始解 % options: 选项 % 默认选项 defaultOptions = struct('MaxIter',1000,'MaxFunEvals',10000,'TolFun',1e-6,'TolX',1e-6,'InitialTemperature',100,'TemperatureFcn',@temperatureboltz,'AnnealingFcn',@annealingfast); % 合并选项 if nargin < 3 options = []; end options = mergeoptions(defaultOptions,options); % 初始化 x = x0; fval = feval(fun,x); T = options.InitialTemperature; iter = 0; funccount = 1; % 主循环 while funccount < options.MaxFunEvals && iter < options.MaxIter && fval > options.TolFun && T > options.TolX % 生成新解 xnew = feval(options.AnnealingFcn,x,T); fnew = feval(fun,xnew); funccount = funccount + 1; % 接受新解 if fnew < fval || rand < exp(-(fnew-fval)/T) x = xnew; fval = fnew; end % 更新温度 T = feval(options.TemperatureFcn,T,iter); % 更新迭代次数 iter = iter + 1; end end function T = temperatureboltz(T,iter) % Boltzmann 降温方案 T = T / log(iter+1); end function xnew = annealingfast(x,T) % 快速模拟退火方案 xnew = x + T*randn(size(x)); end function options = mergeoptions(defaultOptions,options) % 合并选项 if isempty(options) options = defaultOptions; else names = fieldnames(defaultOptions); for i = 1:length(names) if ~isfield(options,names{i}) options.(names{i}) = defaultOptions.(names{i}); end end end end ``` 这个代码实现了一个简单的模拟退火算法,可以用于解决优化问题。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值