模拟退火算法

模拟退火算法（智能算法）

一、算法的思想

​ 模拟退火的简单流程：

接受概率：$p_t=e^{-|f(B)-f(A)|C_t}$ （metropolis准则）

1、优化问题有约束时——》建立惩罚机制（将目标函数增大或减小）

2、$C_t$的设置 ：

模拟退火的概念

$C_t=\frac{1}{T_t};T_{t+1}=\alpha T_t$ 将 $C_t$ 用温度来表示，其中$\alpha$是常数，取0.95

两次循环，不同温度进行搜索（在同一温度下进行搜索，改变温度再进行搜索）

停止搜索的条件：（1）达到设定的迭代次数、（2）温度小于0.0000001、（3）找到最优解重复300次

3、怎么在旧解的附近找到新的解。

（1）求解函数的最优解时，新解的产生

​ ——code2.m

（2）旅行商问题

三种的方法：（1）交换法【随机选择两个点，进行交换】（2）移位法【随机选择3个点，将2个点之间的点移位到第三个点】（3）倒置法【随机选择两个点，将两点之间的点全部颠倒顺序】

%%  例子：TSP问题旅行商问题
Coord = [ 66.83 61.95 40    24.39 17.07 22.93 51.71 87.32 68.78 84.88 50 40 25 ;
    25.36 26.34 44.39 14.63 22.93 76.1  94.14 65.36 52.19 36.09 30 20 26] ;
%程序参数设定
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 ) ; % 路 径 费 用 fare
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 , v ] = 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 ;
end
pathfar = e0 ;
% 输 入 结 果
fprintf( '近似最优路径为：\n ' )
%disp( char( [ path , path(1) ] + 64 ) ) ;
disp(path)
fprintf( '近似最优路径路程\tpathfare=' ) ;
disp( pathfar ) ;
myplot( path , Coord , pathfar ) ;

新解的产生要旧解有关（新解的产生没有同意的定义，第一种是matlab内部函数的定义）。 ——code3.m

（3）书店买书问题

​ 找新解的方法：随机找到一本书改变买书的书店

（3）背包问题 ——code4.m

​ metropolis准则

​ 产生新解的方法

产生新解的三种情况：

价值变化$\Delta f=\{\begin{array}{l}c(i)，将物品i直接装入\\ c(i)-c(j)，将物品i装入且物品j取出\\ c(j)-c(i)，将物品i取出且装入物品j\end{array}$

模拟退火算法更倾向于解决组合优化问题

问题一

代码：

%% 参数初始化
narvs = 2; % 变量个数
T0 = 100;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 300;  % 最大迭代次数
Lk = 200;  % 每个温度下的迭代次数
alfa = 0.88;  % 温度衰减系数   最终
lb =[-3 4.1]; % x的下界
ub =[12.1 5.8]; % x的上界
%%  随机生成一个初始解
x0 = zeros(1,narvs);
for i = 1: narvs
    x0(i) = lb(i) + (ub(i)-lb(i))*rand(1);    
end
y0 = Obj_fun(x0); % 计算当前解的函数值

%% 定义一些保存中间过程的量，方便输出结果和画图
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) < lb(j)
                r = rand(1);
                x_new(j) = r*lb(j)+(1-r)*x0(j);
            elseif x_new(j) > ub(j)
                r = rand(1);
                x_new(j) = r*ub(j)+(1-r)*x0(j);
            end
        end
        x1 = x_new;    % 将调整后的x_new赋值给新解x1
        y1 = Obj_fun(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;   % 温度下降
end
disp('最佳的位置是：'); disp(best_x)
disp('此时最优值是：'); disp(max_y)
%% 画出每次迭代后找到的最大y的图形
figure
plot(1:maxgen,MAXY,'b-');
xlabel('迭代次数');
ylabel('y的值');