17-1利用模拟退火算法进行一元函数和二元函数的极值问题优化（附matlab程序）

1.简述

      

        模拟退火算法来源于固体退火原理，是一种基于概率的算法。模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。

物理退火过程

将固体加温至充分高，再让其徐徐冷却（退火），加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。

算法思想

物理中，淬火（快速降温），会导致不是最低能态的非晶形；退火，可达到最低能量状态。模仿自然界退火现象而得，利用了物理中固体物质的退火过程与一般优化问题的相似性。

模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是一种通用的优化算法，理论上算法具有概率的全局优化性能。

模拟退火算法包含两个部分即Metropolis算法和退火过程。Metropolis算法就是如何在局部最优解的情况下让其跳出来，是退火的基础。1953年Metropolis提出重要性采样方法，即以概率来接受新状态，而不是使用完全确定的规则，称为Metropolis准则，计算量较低。

模拟退火算法新解的产生和接受的步骤：

第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响（冷却进度表(Cooling Schedule)包括控制参数的初值 t 及其衰减因子Δt 、每个 t 值时的迭代次数L和停止条件S）。

第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。

第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若ΔT<0则接受S′作为新的当前解S，否则以概率exp(-ΔT/T)接受S′作为新的当前解S。

第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。

流程图

应用场景

在已知的某个定义域内求函数最优值的问题通常有以下三种情况，以求最小值问题为例（求最大值可以转换为求最小值）：

· 如果是离散取值，则可以用穷举法获得问题最优解

· 如果取值连续，且函数是凸函数，则可以用梯度下降等方法获得最优解

· 如果连续非凸，虽说根据已有的近似求解法能够找到问题解，可解是否是最优的还有待考量，很多时候若初始值选择的不好，非常容易陷入局部最优值

2.代码

%% I. 清空环境变量   optimtool    打开自带的模拟退火工具箱   1.5     100次迭代最大次数    [1 1]    [-5 -5]   [5 5]
clear all
clc

%% II. 一元函数优化
x = 1:0.01:2;
y = sin(10*pi*x) ./ x;
figure
plot(x,y,'linewidth',1.5)
ylim([-1.5, 1.5])
xlabel('x')
ylabel('y')
title('y = sin(10*pi*x) / x')
hold on

%%
% 1. 标记出最大值点
[maxVal,maxIndex] = max(y);
plot(x(maxIndex), maxVal, 'r*','linewidth',2)
text(x(maxIndex), maxVal, {['    X: ' num2str(x(maxIndex))];['    Y: ' num2str(maxVal)]})
hold on

%%
% 2. 标记出最小值点
[minVal,minIndex] = min(y);
plot(x(minIndex), minVal, 'ks','linewidth',2)
text(x(minIndex), minVal, {['    X: ' num2str(x(minIndex))];['    Y: ' num2str(minVal)]})

%% III. 二元函数优化
[x,y] = meshgrid(-5:0.1:5,-5:0.1:5);
z = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20;
figure
mesh(x,y,z)
hold on
xlabel('x')
ylabel('y')
zlabel('z')
title('z =  x^2 + y^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20')

%%
% 1. 标记出最大值点
maxVal = max(z(:));
[maxIndexX,maxIndexY] = find(z == maxVal);
for i = 1:length(maxIndexX)
    plot3(x(maxIndexX(i),maxIndexY(i)),y(maxIndexX(i),maxIndexY(i)), maxVal, 'r*','linewidth',2)
     text(x(maxIndexX(i),maxIndexY(i)),y(maxIndexX(i),maxIndexY(i)), maxVal, {['    X: ' num2str(x(maxIndexX(i),maxIndexY(i)))];['    Y: ' num2str(y(maxIndexX(i),maxIndexY(i)))];['    Z: ' num2str(maxVal)]})
    hold on
end
 

3.运行结果