模拟退火算法matlab代码

爱静的龙猫

于 2024-01-10 17:14:59 发布

阅读量1.4k

点赞数 26

文章标签：模拟退火算法 matlab 数学建模

本文链接：https://blog.csdn.net/2302_79394843/article/details/135500535

版权

智能优化算法

智能优化算法又称现代启发式算法，是一种具有全局优化性能、通用性强且适用于并行处理的算法。这种算法一般具有严密的理论依据，而不是单纯凭专家经验，理论上可以在一定的时间内找到最优解或近似最优解。

模拟退火算法（Simulated Annealing，SA）是一种全局优化算法，通常用于解决组合优化问题。它的灵感来源于固体退火过程，其中材料在高温下热化，然后慢慢冷却，最终达到稳定的低能态。在数学建模中，模拟退火被用来在搜索空间中找到全局最优解，尤其在复杂的多维空间中表现出色。

以下是模拟退火算法的基本思想和步骤：

1. **初始解的生成：** 从问题的搜索空间中随机生成一个初始解，作为当前的候选解。

2. **温度初始化：** 设置一个初始温度，表示系统的热度。温度较高时，算法有更大概率接受差于当前解的新解；温度逐渐降低，接受劣解的概率减小。

3. **状态转移：** 在当前解的邻域中随机选择一个新解。如果新解比当前解更好，则接受新解；如果新解比当前解差，以一定概率接受新解，这个概率由温度和能量差决定。接受差解的概率在开始时较高，随着温度的降低而减小。

4. **温度更新：** 通过降低温度的方式逐渐减小接受差解的概率，模拟系统冷却。温度的降低通常遵循一个退火计划。

5. **迭代：** 重复进行状态转移和温度更新，直到满足停止准则，如达到最大迭代次数或温度降低到某个阈值。

6. **输出结果：** 返回搜索过程中找到的最优解，或者在停止准则满足时输出当前最优解。

模拟退火算法在求解组合优化问题时的主要优势在于，它能够在搜索空间中随机游走，从而避免陷入局部最优解。算法中的“退火”过程允许接受劣解，有助于跳出局部最优解，最终达到全局最优解。虽然算法的性能依赖于参数的选择和退火计划的设计，但合理选择这些参数通常能够使算法在多种问题上表现良好。

%% SA 模拟退火: 求解函数y = 11*sin(x) + 7*cos(5*x)在[-3,3]内的最大值(动画演示）
tic
clear; clc

%% 绘制函数的图形
x = -3:0.1:3;
y = 11*sin(x) + 7*cos(5*x);
figure
plot(x,y,'b-')
hold on  % 不关闭图形，继续在上面画图

%% 参数初始化
narvs = 1; % 变量个数
T0 = 100;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 200;  % 最大迭代次数
Lk = 100;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数
x_lb = -3; % x的下界
x_ub = 3; % x的上界

%%  随机生成一个初始解
x0 = zeros(1,narvs);
for i = 1: narvs
    x0(i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(1);    
end
y0 = Obj_fun1(x0); % 计算当前解的函数值
h = scatter(x0,y0,'*r');  % scatter是绘制二维散点图的函数（这里返回h是为了得到图形的句柄，未来我们对其位置进行更新）

%% 定义一些保存中间过程的量，方便输出结果和画图
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) < 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;    % 将调整后的x_new赋值给新解x1
        y1 = Obj_fun1(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;   % 温度下降
    pause(0.01)  % 暂停一段时间(单位：秒)后再接着画图
    h.XData = x0;  % 更新散点图句柄的x轴的数据（此时解的位置在图上发生了变化）
    h.YData = Obj_fun1(x0); % 更新散点图句柄的y轴的数据（此时解的位置在图上发生了变化）
end

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

pause(0.5)
h.XData = [];  h.YData = [];  % 将原来的散点删除
scatter(best_x,max_y,'*r');  % 在最大值处重新标上散点
title(['模拟退火找到的最大值为', num2str(max_y)])  % 加上图的标题

%% 画出每次迭代后找到的最大y的图形
figure
plot(1:maxgen,MAXY,'b-');
xlabel('迭代次数');
ylabel('y的值');
toc

hold on 是 MATLAB 中的一个命令，它用于将当前图形保持在打开状态，以便继续在上面绘制更多的图形，而不会覆盖已有的图形。

在你的例子中：
matlabCopy code

x = -3:0.1:3; y = 11 * sin(x) + 7 * cos(5 * x); figure plot(x, y, 'b-') hold on

这段代码首先创建了一个包含函数 11sin⁡(�)+7cos⁡(5�)11sin(x)+7cos(5x) 图形的新窗口（figure 命令）。然后，使用 plot 命令在这个窗口中绘制了蓝色实线图形，表示函数在给定范围内的图形。

hold on 使得当前图形保持打开状态，这样你可以在这个图形上继续添加更多的图形，而不会清除已有的内容。如果没有使用 hold on，在下一次使用 plot 命令时，新的图形将会替代当前的图形。

例如，你可以在这个图形上添加更多的曲线、数据点、标签等。当你完成所有需要绘制的内容后，你可以使用 hold off 命令来关闭保持状态。

这段 MATLAB 代码用于初始化模拟退火算法的一些参数。让我们逐行解释：

1. **`narvs = 1;`：** 变量个数。在这个上下文中，`narvs` 表示问题中的变量个数。在这里，变量的个数被设置为 1。

2. **`T0 = 100;`：** 初始温度。模拟退火算法中的温度参数，表示系统的热度。在算法开始时，系统会以这个初始温度开始搜索。

3. **`T = T0;`：** 迭代中温度会发生改变，第一次迭代时温度就是 `T0`。这个变量用于记录当前迭代的温度。

4. **`maxgen = 200;`：** 最大迭代次数。表示算法最多进行多少次迭代，防止算法无限制地执行。

5. **`Lk = 100;`：** 每个温度下的迭代次数。表示在每个温度下，算法会进行多少次状态转移。

6. **`alfa = 0.95;`：** 温度衰减系数。在每次迭代结束后，当前温度会以 `alfa` 的比例进行衰减，模拟系统冷却的过程。

7. **`x_lb = -3;`：** 变量 x 的下界。在问题中，x 的取值范围被限定在 -3 到 3 之间。

8. **`x_ub = 3;`：** 变量 x 的上界。表示 x 的取值上限。

这些参数主要用于配置模拟退火算法的初始条件和运行参数，以便在解决具体问题时进行调整。在你的代码中，这些参数可能是用于解决一个特定问题的模拟退火算法。

这部分代码是为了确保生成的新解 `x_new` 在定义域内。让我们逐行解释它的含义：

1. **`for j = 1: narvs`：** 这是一个循环，用于遍历每个解向量的分量（变量）。

2. **`if x_new(j) < x_lb(j)`：** 这个条件判断是否新生成的解向量中的第 `j` 个分量小于定义域的下界 `x_lb(j)`。

   - 如果是，就执行下面的语句块。这意味着新解的某个分量小于下界，需要进行调整。

3. **`r = rand(1);`：** 生成一个在 `[0, 1)` 范围内的随机数 `r`。

4. **`x_new(j) = r*x_lb(j) + (1-r)*x0(j);`：** 通过线性插值的方式调整新解的第 `j` 个分量，确保它大于等于下界 `x_lb(j)`。具体地，`r*x_lb(j)` 是下界的一部分，`(1-r)*x0(j)` 是当前解的一部分。通过这种方式，新解的第 `j` 个分量会在下界和当前解之间取值。

   - 这个调整确保了新解在该分量上不会小于定义域的下界。

5. **`elseif x_new(j) > x_ub(j)`：** 这个条件判断是否新生成的解向量中的第 `j` 个分量大于定义域的上界 `x_ub(j)`。

   - 如果是，就执行下面的语句块。这意味着新解的某个分量大于上界，需要进行调整。

6. **`r = rand(1);`：** 同样，生成一个在 `[0, 1)` 范围内的随机数 `r`。

7. **`x_new(j) = r*x_ub(j) + (1-r)*x0(j);`：** 通过线性插值的方式调整新解的第 `j` 个分量，确保它小于等于上界 `x_ub(j)`。具体地，`r*x_ub(j)` 是上界的一部分，`(1-r)*x0(j)` 是当前解的一部分。通过这种方式，新解的第 `j` 个分量会在上界和当前解之间取值。

   - 这个调整确保了新解在该分量上不会大于定义域的上界。

综合起来，这段代码的目的是确保生成的新解 `x_new` 在每个分量上都处于定义域 `[x_lb(j), x_ub(j)]` 内。这有助于保持搜索空间的有效性。

这样的调整是为了确保生成的新解 `x_new` 落在定义域内。在模拟退火算法中，通过随机扰动和接受劣质解的策略来遍历解空间，但在这个过程中，生成的新解可能会超出问题定义域。

通过线性插值的方式调整，是为了在处理超出边界的情况时保持随机性。具体地，随机生成一个 `[0, 1)` 范围内的值 `r`，然后用这个值对当前解 `x0` 和边界值 `x_lb` 或 `x_ub` 进行线性插值。这个过程确保了调整后的新解在当前解和边界值之间取值，并且这种插值的方式不仅保持了一定的随机性，而且确保了解的合理性。

在算法中，这种调整确保了搜索过程中生成的新解保持在问题的定义域内，使得算法能够在解空间内自由搜索，并通过随机性来尽可能避免陷入局部最小值。

这部分代码是模拟退火算法的核心迭代逻辑，用于接受或拒绝新的解，并在算法中进行状态更新。下面逐行解释每句代码：

1. **`x1 = x_new;`：** 将经过调整的新解 `x_new` 赋值给 `x1`。这一步是为了在下面的判断和更新中使用新的解。

2. **`y1 = Obj_fun1(x1);`：** 计算新解 `x1` 的函数值，得到 `y1`。这一步是为了比较新解和当前解在目标函数上的表现。

3. **`if y1 > y0`：** 判断新解的函数值是否大于当前解的函数值。

- 如果是，表示找到了更好的解，执行下面的语句块。

4. **`x0 = x1;`：** 更新当前解为新解，将 `x1` 赋值给 `x0`。

5. **`y0 = y1;`：** 更新当前解的函数值为新解的函数值。

- 这表示算法接受了更好的解。

6. **`else`：** 如果新解的函数值不大于当前解的函数值，执行下面的语句块。

7. **`p = exp(-(y0 - y1)/T);`：** 根据Metropolis准则计算一个接受劣质解的概率 `p`。Metropolis准则是模拟退火算法中用于接受劣质解的概率计算方法。

8. **`if rand(1) < p`：** 生成一个随机数和接受劣质解的概率 `p` 进行比较。

- 如果生成的随机数小于概率 `p`，表示接受劣质解，执行下面的语句块。

9. **`x0 = x1;`：** 更新当前解为新解，将 `x1` 赋值给 `x0`。

10. **`y0 = y1;`：** 更新当前解的函数值为新解的函数值。

- 这表示算法以一定的概率接受劣质解。

总体而言，这段代码实现了模拟退火算法中的接受或拒绝新解的逻辑。如果新解更好，则直接接受；如果新解较差，根据Metropolis准则以一定概率接受。这个策略有助于算法在搜索过程中在全局范围内探索，并以较小的概率接受劣质解，以避免陷入局部最小值。

这部分代码用于判断是否要更新已找到的最佳解，即是否有更好的解出现。下面逐行解释每句代码：

1. **`if y0 > max_y`：** 判断当前解的函数值 `y0` 是否大于已知的最大函数值 `max_y`。

- 如果是，表示当前解更好，执行下面的语句块。

2. **`max_y = y0`：** 更新最大的函数值为当前解的函数值。

3. **`best_x = x0`：** 更新找到的最佳解 `best_x` 为当前解 `x0`。

- 这表示算法发现了一个更好的解，因此更新已知的最佳解和最大函数值。

这段代码的作用是在每次内循环结束后，检查当前解是否比已知的最佳解更好。如果是，则更新最佳解的信息，确保最佳解是目前已搜索到的最优解。这有助于在算法运行过程中跟踪并记录找到的最佳解。

这部分代码涉及到模拟退火算法的一些辅助操作和图形更新，下面逐行解释每句代码：

1. **`MAXY(iter) = max_y;`：** 将当前外循环迭代结束后找到的最大函数值 `max_y` 记录到数组 `MAXY` 的第 `iter` 个位置。这个数组用于跟踪每轮迭代结束后找到的最大函数值。

2. **`T = alfa*T;`：** 温度 `T` 下降，根据模拟退火算法的温度下降策略。温度下降是模拟退火算法的核心机制之一，用于控制搜索空间的随机性，使得在搜索开始时更容易接受劣质解，随着搜索进行逐渐降低接受劣质解的概率。

3. **`pause(0.01)`：** 暂停一段时间，单位是秒。这个操作主要是为了在图形更新中观察到算法的迭代过程，以便更好地理解算法的行为。这是一个可选的操作，对算法本身的优化并没有影响。

4. **`h.XData = x0;`：** 更新散点图句柄 `h` 的 x 轴数据为当前解 `x0`。这表示在图形上更新当前解的位置。

5. **`h.YData = Obj_fun1(x0);`：** 更新散点图句柄 `h` 的 y 轴数据为当前解的函数值。这表示在图形上更新当前解的函数值。

这段代码主要用于记录每轮外循环结束后找到的最大函数值，同时通过图形更新操作观察算法的迭代过程。

这部分代码主要是用于在算法运行结束后展示最终的结果，包括输出最佳位置和最优值，并通过图形的方式将最佳位置标记出来。下面逐行解释每句代码：

1. **`disp('最佳的位置是：'); disp(best_x)`：** 在命令行窗口输出字符串 "最佳的位置是："，然后输出变量 `best_x` 的值，即找到的最佳解的位置。

2. **`disp('此时最优值是：'); disp(max_y)`：** 在命令行窗口输出字符串 "此时最优值是："，然后输出变量 `max_y` 的值，即找到的最优解的函数值。

3. **`pause(0.5)`：** 暂停 0.5 秒，给用户一定时间查看命令行中的输出。

4. **`h.XData = []; h.YData = [];`：** 将原来图形上的散点数据清空。这是为了清除之前迭代过程中的散点，为最终结果的展示做准备。

5. **`scatter(best_x,max_y,'*r');`：** 在最佳位置处重新标上一个红色星号的散点。这样图形上就只有最终找到的最佳位置标记。

6. **`title(['模拟退火找到的最大值为', num2str(max_y)])`：** 添加图的标题，显示模拟退火算法找到的最大值。标题中包含最大值的数值。

这段代码的作用是在算法运行结束后，清空原来的图形上的散点数据，将最佳位置用红色星号标记出来，并在命令行窗口输出最佳位置和最优值的信息。

这段代码用于绘制模拟退火算法每次外循环结束后找到的最大函数值随迭代次数的变化图。下面逐行解释每句代码：

1. **`figure`：** 创建一个新的图形窗口，用于绘制图形。

2. **`plot(1:maxgen, MAXY, 'b-')`：** 绘制折线图，横轴是迭代次数，纵轴是每次迭代后找到的最大函数值。`1:maxgen` 表示横轴的范围是从 1 到 `maxgen`，`MAXY` 是纵轴上的数值，`'b-'` 表示使用蓝色的实线连接数据点。

3. **`xlabel('迭代次数')`：** 设置横轴的标签为 "迭代次数"。

4. **`ylabel('y的值')`：** 设置纵轴的标签为 "y的值"。

5. **`toc`：** 输出程序运行的总时间，即自开始计时以来经过的时间。

这段代码的作用是通过折线图展示模拟退火算法每次外循环结束后找到的最大函数值随迭代次数的变化趋势。