模拟退火算法(Simulated Annealing,SA)是一种全局优化算法,其基本思想是通过一定的概率接受劣解,以避免陷入局部最优解。它模拟了物质固体退火时的过程,即将物质加热至高温状态,然后缓慢冷却,使其达到稳定状态。在优化问题中,这个过程被用来搜索全局最优解。
Matlab是一种功能强大的数学软件,它提供了许多优化算法的实现。在本文中,我们将介绍如何使用Matlab实现模拟退火算法,并提供15个实际应用的源码案例。
1. 算法步骤
模拟退火算法的基本步骤如下:
-
初始化参数:
- 初始温度(
T0
) - 降温速率(
alpha
) - 终止温度(
T_min
) - 初始解(
x0
)
- 初始温度(
-
产生新解:在当前解的邻域内产生一个新解。
-
接受新解:计算当前解与新解之间的差异,如果新解更优,则接受它;否则,以一定的概率接受它。
-
降温:根据设定的降温速率降低温度。
-
终止判断:如果温度降低到终止温度以下,则停止搜索,输出最优解。
2. 代码实现
以下是Matlab实现模拟退火算法的基本代码:
function [x,fval] = simulated_annealing(fun,x0,options)
% fun: 目标函数句柄
% x0: 初始解
% options: 选项结构体,包括以下字段:
% T0: 初始温度
% alpha: 降温速率
% T_min: 终止温度
% max_iter: 最大迭代次数
% verbose: 是否打印输出
% 返回值:
% x: 最优解
% fval: 目标函数在最优解处的取值
% 设置默认选项
default_options = struct('T0',100,'alpha',0.95,'T_min',1e-8,...
'max_iter',1000,'verbose',false);
if nargin < 3
options = default_options;
else
options = merge_options(default_options,options);
end
% 初始化参数
T = options.T0;
x = x0;
fval = feval(fun,x);
iter = 0;
best_x = x;
best_fval = fval;
% 开始迭代
while T > options.T_min && iter < options.max_iter
% 产生新解
new_x = x + randn(size(x));
new_fval = feval(fun,new_x);
delta_f = new_fval - fval;
% 接受新解
if delta_f < 0 || exp(-delta_f/T) > rand()
x = new_x;
fval = new_fval;
if fval < best_fval
best_x = x;
best_fval = fval;
end
end
% 降温
T = options.alpha * T;
% 打印输出
if options.verbose
fprintf('iter=%d, T=%g, fval=%g, best_fval=%g\n',iter,T,fval,best_fval);
end
% 更新迭代计数器
iter = iter + 1;
end
% 返回最优解和目标函数值
x = best_x;
fval = best_fval;
% 合并选项结构体
function opt = merge_options(default_opt,opt)
if isempty(opt)
opt = default_opt;
else
fields = fieldnames(default_opt);
for i = 1:length(fields)
if ~isfield(opt,fields{i})
opt.(fields{i}) = default_opt.(fields{i});
end
end
end
end
end
3. 示例
我们将使用15个不同的优化问题来演示模拟退火算法的应用。每个案例都附有相应的Matlab代码。
示例 1: Rosenbrock函数
Rosenbrock函数是一个经典的非凸函数,其表达式为:
[ f(x) = \sum_{i=1}^{n-1} [100 (x_{i+1} - x_i2)2 + (1 - x_i)^2] ]
fun = @(x) sum(100*(x(2:end)-x(1:end-1).^2).^2 + (1-x(1:end-1)).^2);
x0 = [-1.2; 1];
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(fun, x0, options);
disp(['Best solution: ', num2str(x')]);
disp(['Objective function value: ', num2str(fval)]);
示例 2: 旅行商问题(TSP)
旅行商问题的目标是找到一条最短路径,使得旅行商能访问所有城市并返回起点。
cities = [0,0; 1,0; 1,1; 0,1]; % Example cities coordinates
dist = @(path) sum(sqrt(sum(diff(cities(path,:),1,1).^2,2) + (norm(cities(path(end),:) - cities(path(1),:)))^2));
fun = @(path) dist(path);
x0 = 1:length(cities);
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(@(path) fun(mod(path-1+randperm(length(path)),length(path))+1), x0, options);
disp(['Best path: ', num2str(x)]);
disp(['Objective function value: ', num2str(fval)]);
示例 3: 函数最小化问题
对于一个简单的二次函数进行最小化:
[ f(x) = (x - 3)^2 + 5 ]
fun = @(x) (x - 3).^2 + 5;
x0 = 0;
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(fun, x0, options);
disp(['Best solution: ', num2str(x)]);
disp(['Objective function value: ', num2str(fval)]);
示例 4: 约束优化问题
带约束的优化问题,如最小化 ( f(x) = x^2 ) subject to ( x \geq 0 ):
fun = @(x) x.^2;
constraint = @(x) max(0, -x); % Constraint x >= 0
x0 = -1;
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(@(x) fun(x) + 1e6 * constraint(x), x0, options);
disp(['Best solution: ', num2str(x)]);
disp(['Objective function value: ', num2str(fval)]);
示例 5: 非线性规划问题
优化非线性目标函数 ( f(x, y) = (x - 1)^2 + (y - 2)^2 ):
fun = @(xy) (xy(1) - 1)^2 + (xy(2) - 2)^2;
x0 = [0; 0];
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(fun, x0, options);
disp(['Best solution: ', num2str(x')]);
disp(['Objective function value: ', num2str(fval)]);
示例 6: 数据拟合问题
拟合数据点 ( (x_i, y_i) ) 的线性模型:
x_data = [1, 2, 3, 4, 5];
y_data = [2.1, 2.9, 3.8, 5.0, 5.1];
model = @(p, x) p(1) * x + p(2);
error = @(p) sum((y_data - model(p, x_data)).^2);
fun = @(p) error(p);
x0 = [1; 1];
options = struct('T0',100,'alpha',0.95,'T_min',
1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(fun, x0, options);
disp(['Best parameters: ', num2str(x')]);
disp(['Objective function value: ', num2str(fval)]);
示例 7: 最优控制问题
控制一个系统使其状态最优化:
control = @(u) sum(u.^2); % Example control objective
u0 = [0; 0];
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(control, u0, options);
disp(['Best control: ', num2str(x')]);
disp(['Objective function value: ', num2str(fval)]);
示例 8: 组合优化问题
组合问题,例如最大化一个简单的目标函数:
fun = @(x) sum(x); % Example: maximize sum of elements
x0 = [0, 0, 0];
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(fun, x0, options);
disp(['Best solution: ', num2str(x)]);
disp(['Objective function value: ', num2str(fval)]);
示例 9: 多目标优化问题
优化两个目标函数的加权和:
fun = @(x) [x(1)^2 + x(2)^2, (x(1) - 1)^2 + (x(2) - 1)^2];
weights = [0.5, 0.5];
weighted_sum = @(x) sum(weights .* fun(x));
x0 = [0, 0];
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(weighted_sum, x0, options);
disp(['Best solution: ', num2str(x)]);
disp(['Objective function value: ', num2str(fval)]);
示例 10: 排程问题
优化任务排程的目标函数:
task_times = [5, 2, 7, 3]; % Example task times
fun = @(schedule) sum(task_times(schedule));
x0 = 1:length(task_times);
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(@(schedule) fun(mod(schedule-1+randperm(length(schedule)),length(schedule))+1), x0, options);
disp(['Best schedule: ', num2str(x)]);
disp(['Objective function value: ', num2str(fval)]);
示例 11: 函数拟合问题
拟合数据点 ( (x_i, y_i) ) 到一个多项式模型:
x_data = [1, 2, 3, 4, 5];
y_data = [1.1, 1.9, 3.0, 4.1, 5.0];
poly_model = @(p, x) polyval(p, x);
error = @(p) sum((y_data - poly_model(p, x_data)).^2);
fun = @(p) error(p);
x0 = [0; 0; 0];
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(fun, x0, options);
disp(['Best polynomial coefficients: ', num2str(x')]);
disp(['Objective function value: ', num2str(fval)]);
示例 12: 线性规划问题
解决线性规划问题:
c = [1, 2];
A = [1, 1; -1, 2];
b = [6; 4];
fun = @(x) c * x';
constraints = @(x) [A * x' - b; -x];
x0 = [0; 0];
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(@(x) fun(x) + 1e6 * sum(max(0, constraints(x))), x0, options);
disp(['Best solution: ', num2str(x')]);
disp(['Objective function value: ', num2str(fval)]);
示例 13: 动态规划问题
求解一个动态规划问题:
fun = @(x) sum(x.^2); % Example dynamic programming objective
x0 = [1, 1, 1];
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(fun, x0, options);
disp(['Best solution: ', num2str(x)]);
disp(['Objective function value: ', num2str(fval)]);
示例 14: 图着色问题
优化图的着色问题,最小化相邻节点的相同颜色:
n = 4; % Number of nodes
edges = [1,2; 2,3; 3,4; 4,1]; % Example edges
colors = 1:3; % Possible colors
coloring = @(c) sum(arrayfun(@(i) any(c(edges(i,1)) == c(edges(i,2))), 1:size(edges,1)));
fun = @(c) coloring(c);
x0 = randi([1, 3], [n, 1]);
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(@(c) fun(c), x0, options);
disp(['Best coloring: ', num2str(x')]);
disp(['Objective function value: ', num2str(fval)]);
示例 15: 库存管理问题
优化库存管理中的订单量:
order_cost = @(q) (q - 10)^2 + 5; % Example cost function
x0 = 10; % Initial order quantity
options = struct('T0',100,'alpha',0.95,'T_min',1e-8,'max_iter',1000,'verbose',true);
[x, fval] = simulated_annealing(order_cost, x0, options);
disp(['Best order quantity: ', num2str(x)]);
disp(['Objective function value: ', num2str(fval)]);
4. 总结
本文详细介绍了如何在Matlab中实现模拟退火算法,并通过15个具体案例展示了其在各种优化问题中的应用。模拟退火算法是一种有效的全局优化方法,适用于求解复杂的非凸问题。通过Matlab的灵活实现和强大的功能,可以高效地应用模拟退火算法解决实际问题。