退火算法以及多元函数最优值的求解

启发式算法

启发式算法的一个重要的特点就是在搜索最优解的过程中利用到了原来搜索过程中得到的信息，利用之前的信息改进我们的搜索过程。

爬山法

属于启发式算法的一种简单算法，网上有大佬给出了爬山法，退火算法，遗传算法，禁忌搜索的通俗解释，这里借用一下。

为了找出地球上最高的山，一群有志 的兔子们开始想办法。

(1)兔子朝着比现在高的地方跳去。他们找到了不远处的最高山峰。但是这座山不一定是珠穆朗玛峰。这就是爬山法，它不能保证局部最优值就是全局最优值。

(2)兔子喝醉了。他随机地跳了很长时间。这期间，它可能走向高处,也可能踏入平地。但是,他渐渐清醒了并朝他踏过的最方向跳去。这就是模拟退火。

(3)兔子们知道一个兔的力量是渺小的。他们互相转告着，哪里的山已经找过，并组找过的每一座山他们都留 下一只兔子做记号。他们制定了下一步去哪里寻找的策略。这就是禁忌搜索。

(4)兔子们吃了失忆药片，并被发射到太空,然后随机落到了地球上的某些地方。他们不知道自己的使命是什么。但是,如果你过几年就杀死一部分海拔低的兔子， 多产的兔子们自己就会找到珠穆朗玛峰。这就是遗传算法。

爬山法基本思路

以连续一元函数找最大值为例

1. 在解空间上随机初始化一个初始解
2. 根据初始解的位置，向左或者向右走一点，当然这个步距还是有要求的。步距过长，可能会跨过最优值；步距过小，会导致爬山时间过长，效率下降。
3. 如果所在的位置的值比之前大，那么就不断爬山。直到开始走下坡，那么开始走下坡的时候，我们就把上一步的值作为一个最优值。

需要注意的是，爬山法特别容易找到局部最优解，这取决于初始值的位置。在实际应用中存在许多局部最优解的函数，所以爬山法的缺陷很大。

模拟退火

思想

模拟退火算法的思想来源于对固体退火降温过程的模拟。即将固体加温至充分高，再让其徐徐冷却。在加热固体时，固体中原子的热运动不断增强，内能增大，随着温度的不断升高，固体的长程有序被彻底破坏，固体内部粒子随温度的升高而变为无序状。冷却时，粒子逐渐趋于有序，在每个温度下都达到平衡状态，最后在常温下达到基态，同时内能也减为最小。

转化

只看物理意义上的退火，可能有点无法理解，我们需要将物理退火的模型转化到模拟退火的模型

物理退火	模拟退火
粒子状态	解
能量最低态	最优解
熔解过程	设定初温
冷却	控制参数的下降
能量	目标函数

与爬山法相比

我们知道爬山法对于小于当前的位置的地方进行直接抛弃，这也是其容易找到局部最优值的原因。退火算法与其不同的地方就是按照一定的概率来接受那个比当前位置小的地方，向那里爬山。这就是退火算法中的接受准则。

Metropolis准则

下面为下一步的值减去当前值的数值
$$△E=NT-OT$$

如果E>0的话，我们接受这个位置，否则我们以概率exp(-E/T) 的概率来接受这个位置（T指的是温度）。

我们可以注意到，温度是在变化的，也就是说，接受的概率是在变化的，随着温度的降低，接受的概率也就会变小。即随着算法的运行，温度不断下降，对于小于当前位置的值的位置接受概率逐渐变小。

流程

实际应用

我们现在应用退火算法来求解有三个变量的函数最优值
$$y=5\ast sin(x1)-6\ast cos(5\ast x2)+3\ast x3;$$
变量约束范围如下：
$$\begin{gathered} -10\le x1\le 3 \\ -14\le x2\le 5 \\ -7\le x3\le 8 \end{gathered}$$

完整代码如下：

fum.m

function y = fun(x)
    y = 5*sin(x(1)) - 6*cos(5*x(2)) + 3*x(3);
end

SA.m

clear; clc
%% 参数初始化
narvs = 3; % 变量个数
T0 = 100;   % 初始温度
T = T0; 
maxgen = 1000;  % 最大迭代次数
Lk = 100;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数
x_lb = [-10,-14, -7]; % x的下界
x_ub = [3,5, 8]; % x的上界

%%  随机生成一个初始解
x0 = zeros(1,narvs);
for i = 1: narvs
    x0(i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(1);    
end
y0 = fun(x0); 


%% 定义一些保存中间过程的量，方便输出结果和画图
max_y = y0;     % 初始化找到的最佳的解对应的函数值为y0
MAXY = zeros(maxgen,1); % 记录每一次外层循环结束后找到的max_y (方便画图）

%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  % 内循环，在每个温度下开始迭代
        % 产生随机新解
        y = randn(1,narvs);  
        z = y / sqrt(sum(y.^2)); 
        x_new = x0 + z*T; 
        % 对新解的位置进行调整
        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;    
        y1 = fun(x1); 
        if y1 > y0    % 如果新解函数值大于当前解的函数值
            x0 = x1; 
            y0 = y1;
        else
            % 根据Metropolis准则判断是否接受这个解
            p = exp(-(y0 - y1)/T); 
            if rand(1) < p  
                x0 = x1; 
                y0 = y1;
            end
        end
        % 判断是否有更，好值更新解
        if y0 > max_y  
            max_y = y0; 
            best_x = x0;  
        end
    end
    MAXY(iter) = max_y; 
    T = alfa*T;   
    pause(0.01)  
end

disp('最佳的位置是：'); disp(best_x)
disp('此时最优值是：'); disp(max_y)

运行如下：

最佳的位置是：
    1.5702   -5.6550    8.0000

此时最优值是：
   35.0000