基于模拟退火算法的函数寻优算法

一，一元函数寻优

以最小值问题为例，例如计算以下函数的最小值：$$f(x)=\frac{sin(10\pi x)}{x}$$$$s.t.x\in [1,3]$$MATLAB程序实现：

%%% 求一元函数 y=sin(10*pi*x) / x ,x∈[1, 3] 的最小值
%% main.m 
clear all
clc
figure(1);
hold on
ezplot('sin(10*pi*x)/x', [1, 3]);
%% I. 初始化参数
T0 = 1e10;             % 初始温度
Tend = 1e-10;       % 终止温度
L = 2;                    % 各温度下的迭代次数(链长)
q = 0.9;                 % 降温速率
syms x;
Time = ceil(double(solve(T0*(0.9)^x == Tend)));    % 计算迭代的次数
count = 0;        %迭代计数
N = 1;           % 初始解个数
Obj = zeros(Time,1);           % 目标值矩阵初始化
track = zeros(Time, N);       % 每代的最优解矩阵初始化

%% II. 随机得到初始解
x1 = 2 * rand + 1;          

%% III.迭代优化
while T0 > Tend
    count = count + 1;       % 更新迭代次数
    temp = zeros(L, N+1);
    for k = 1:L
        % 1.扰动产生新解
        while 1
            x2 = x1 + 2 * rand - 1;
            % 重复得到新解，直到产生的新解满足约束条件
            if (x1>=1 && x1 <=3) && (x2>=1 && x2<=3)
                break;
            end
        end
        % 2.Metropolis准则
        fx1 = f(x1);     % 计算当前解函数值
        fx2 = f(x2);     % 计算新解函数值
        dC = fx2 - fx1;     % 计算能力之差
        if dC < 0           % 如果能力降低，接受新解
            fx1 = fx2;
            x1 = x2;
        elseif exp(-dC/T0) >= rand  % 以exp(-dC/T)概率接受新解
            fx1 = fx2;
            x1 = x2;
        end
        temp(k, :) = [x1 fx1];      % 记录下一次优化的解和函数值
    end
    %% 3. 记录每次迭代过程的最优路线
    [fmin, index] = min(temp(:, end));    % 找出当前温度下最优解
    if count == 1 || fmin <= Obj(count-1)
        Obj(count) = fmin;                      % 如果当前温度下最优解小于上一过程则记录当前过程
    else
        Obj(count) = Obj(count-1);         % 如果当前温度下最优解大于上一过程则记录上一过程
    end
    track(count, :) = temp(index, 1:end-1);  % 记录当前温度的最优解
    % 
    figure(1);
    plot(track(count, :), fmin, '.', 'color', [count/Time, 0, 0]);
    hold on;
    %降温
    T0 = q * T0;
end
%% IV.结果显示
disp(['自变量：', num2str(track(end))]);
disp(['最小值：', num2str(Obj(end))]);
%% V.绘图
plot(track(end, :), Obj(end), 'ro', 'MarkerSize', 100);
xlabel 'x';
ylabel 'y';
title 'SA算法迭代过程中最优坐标移动';
figure(2);
plot(1:count,Obj)
xlabel('迭代次数')
ylabel('函数值')
title('优化过程')
%% 目标函数
function y = f(x)
y=sin(10*pi*x) / x;
end

Command Window显示的结果为：

自变量：1.1477
最小值：-0.86908

绘图的结果为：

二，二元函数寻优

以最小值问题为例，例如计算以下函数的最小值：$$z=x_2+y_2-10cos(2\pi x)-10cos(xy)+20$$$$s.t.x\in [-5,5],y\in [-5,5]$$MATLAB程序实现：

%%% 求二元函数y=x^2 + y^2 - 10*cos(2*pi*x) - 10*cos(x*y) + 20;
%%% (x, y)∈[-5,5]×[-5,5]
%% main.m 
clear all
clc
figure(1);
[x,y] = meshgrid(-5:0.1:5,-5:0.1:5);
z = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(x*y) + 20;
figure(1);
mesh(x,y,z)
hold on
xlabel('x')
ylabel('y')
zlabel('z')
title('z =  x^2 + y^2 - 10cos(2\pix) - 10cos(xy) + 20')
%% I. 初始化参数
T0 = 1e10;             % 初始温度
Tend = 1e-30;       % 终止温度
L = 2;                    % 各温度下的迭代次数(链长)
q = 0.9;                 % 降温速率
syms x;
Time = ceil(double(solve(T0*(0.9)^x == Tend)));    % 计算迭代的次数
count = 0;        %迭代计数
N = 2;           % 初始解个数
Obj = zeros(Time,1);           % 目标值矩阵初始化
track = zeros(Time, N);       % 每代的最优解矩阵初始化

%% II. 随机得到初始解
x1 = 10 * rand(1, 2) - 5;    

%% III.迭代优化
while T0 > Tend
    count = count + 1;       % 更新迭代次数
    temp = zeros(L, N+1);
    for k = 1:L
        % 1.扰动产生新解
        x2 = x1./5 + 8*rand(1, 2) - 4;
        % 2.Metropolis准则
        fx1 = f(x1(1), x1(2));     % 计算当前解函数值
        fx2 = f(x2(1), x2(2));     % 计算新解函数值
        dC = fx2 - fx1;     % 计算能力之差
        if dC < 0           % 如果能力降低，接受新解
            fx1 = fx2;
            x1 = x2;
        elseif exp(-dC/T0) >= rand  % 以exp(-dC/T)概率接受新解
            fx1 = fx2;
            x1 = x2;
        end
        temp(k, :) = [x1 fx1];      % 记录下一次优化的解和函数值
    end
    %% 3. 记录每次迭代过程的最优路线
    [fmin, index] = min(temp(:, end));    % 找出当前温度下最优解
    if count == 1 || fmin <= Obj(count-1)
        Obj(count) = fmin;                      % 如果当前温度下最优解小于上一过程则记录当前过程
    else
        Obj(count) = Obj(count-1);         % 如果当前温度下最优解大于上一过程则记录上一过程
    end
    track(count, :) = temp(index, 1:end-1);  % 记录当前温度的最优解
    % 
    figure(2);
    plot(track(count, 1), track(count, 2), '.', 'color', [count/Time, 0, 0]);
    hold on;
    %降温
    T0 = q * T0;
end
%% IV.结果显示
disp(['自变量：', num2str(track(end, :))]);
disp(['最小值：', num2str(Obj(end))]);
%% V.绘图
plot(track(end, 1), track(end, 2), 'ro', 'MarkerSize', 50);
xlabel 'x';
ylabel 'y';
title 'SA算法迭代过程中最优坐标移动';
figure(3);
plot(1:count,Obj)
xlabel('迭代次数')
ylabel('函数值')
title('优化过程')
%% 目标函数
function z = f(x, y)
z = x^2 + y^2 - 10*cos(2*pi*x) - 10*cos(x*y) + 20;
end

Command Window中显示的结果为：

自变量：-0.033329   -0.048644
最小值：0.000386

绘图显示为：