国赛培训——最优化智能算法——模拟退火

引入

1.1：模拟退火与蒙特卡洛模拟的区别

1.1 蒙特卡洛模拟：
① 随机生成多组解–>② 验证是否满足约束条件–>③ 将满足约束条件的加入“可行集”中，并计算出可行集中每个解对应的目标函数值–>④ 找出最优解

1.2 区别
①蒙特卡洛模拟法在搜索过程中不改进搜索策略，称为盲目搜索。
②启发式搜索算法会利用中间信息改进搜索策略。

1.2：爬山法与模拟退火法的区别

2.1 爬山法：
① 在解空间随机生成一个初始解
② 向左右邻域走一步，比较大小选择移动方向
③ 不断重复步骤，直到找到一个极大值点
缺点：容易找到局部最优解（找到不符合优化目标的就退出循环）

1.3：基于爬山法的解决方案

为了解决爬山法容易找到局部最优解的情况，设定概率p(0<p<1)，表示接受非当前优化目标解的概率。
（1）f(xi)与f(xj)越接近，概率p越大

（2）接受新解的概率p越大意味着解空间搜索范围越大
（3）搜索前期搜索范围较大，搜索后期搜索范围较小。即前期概率p大，后期概率p小

（4）C_t表示时间相关的系数，时间越长，p越小，C_t越大。C_t关于t递增。

2：模拟退火算法

2.1 算法流程

假设求最大值的问题
（1）随机生成解A，计算A对应的目标函数值f(A)
（2）在A附近随机生成B，验证B是否满足约束条件后，计算B对应的目标函数值f(B)
（3）若f(B) > f(A)，f(A) = f(B)，作为当前最优解
（4）若f(B) < f(A)，计算接受概率P。生成随机数r∈[0,1]，若r < P则接受解B并将B值赋值给A

2.2 C_t的设置

其中T_0=100与α=0.95为初始数据设置，C_t = 1 / T_t

2.3 时间t的实现与停止搜索的条件

（1）通过循环实现t的迭代
通过迭代次数表示时间t的变化，并在同一个t下进行多次搜索。可以通过两层循环实现
for i=1:Max %时间t
for j=1:N %表示相同t内的搜索次数

（2）什么时候停止搜索
① 达到指定迭代次数（循环执行完毕）
② 达到指定温度，比如T < 0.000001
③ 找到最优解连续M次不改变

2.4 如何产生新解（最关键，可通过查询论文或自己制定）

（1）matlab内置函数定义方法

y = randn(1,narvs);			%生成N(0,1)的随机数
z = y / sqrt(sum(y .^2));	%生成z
x_new = x0 + T*z;			%根据规则生成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

3. 例题学习

步骤一：绘制图形（限制在三维以下）
步骤二：参数初始化
步骤三：随机生成初始解并scatter标点
步骤四：定义中间变量
步骤五：模拟退火算法
步骤六：画出最大y值图形

3.1 含上下界约束的规非线性划问题

例题1：

SA 模拟退火: 求解函数y = 11sin(x) + 7cos(5*x)在[-3,3]内的最大值

产生新解的方法：matlab内置函数生成方法 x_new = x0+z*T
① code.m

clc,clear;
%% 第一步：绘制函数图形
x = -3:0.1:3;
y = 11*sin(x) + 7*cos(5*x);
plot(x,y,'b-')
hold on					%不关闭图形

%% 第二步：参数初始化
%需修改项
narvs = 1;				%变量个数（需修改）
maxgen = 200;			%最大迭代次数
Lk = 100;				%每个温度下的搜索次数
x_lb = -3;				%变量下界
x_ub = 3;				%变量上界

%不需修改项
T0 = 100;				%初始温度
T = T0；				%第一次迭代时的温度
alfa = 0.95;

%% 第三步：随机生成初始解 unifrnd(lb,ub)  标点scatter
x0 = zeros(1,narvs);
for i=1:narvs
	x0(i) = unifrnd(x_lb(i),x_ub(i));
end
y0 = Obj_fun1(x0);
h = scatter(x0,y0,'*r');  % scatter是绘制二维散点图的函数

%% 第四步：定义中间变量
max_y = y0;					%初始化找到的最佳解对应的函数值y0
MAXY = zeros(maxgen,1);		%每一次外层循环求得的max_y

%% 第五步：模拟退火过程
% 双重循环迭代查找
for iter= 1:maxgen
	for j=1:Lk
		%(1)产生新解，此处的方法是matlab内置函数方法
		y = randn(1,narvs);		%生成N(0,1)的随机数
		z = y / sqrt(sum(y .^2));
		x_new = x0 + T*z;		%生成新解公式
		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
		%(2)生成新解函数值
		x1 = x_new;
		y1 = Obj_fun1(x1);
		%(3)搜索，当满足优化目标更新；不满足则求解接受概率
		if y1 > y0
			x0 = x1;
			y0 = y1;
		else
			p = exp(-(y0-y1)/T); % 根据Metropolis准则计算一个概率
			if rand(1) < p
				x0 = x1;
				y0 = y1;
			end
		end
		%(4)判断本轮内部搜索是否为最佳解
			if y0 > max_y
				max_y = y0;
				best_x = x0;
			end
	end
	%(5)一个温度情况结束，更新温度与最优值并画图
	MAXY(iter) = max_y;
	T = T*alfa;
	pause(0.01)
	h.XData = x0;
	h.YData = Obj_fun1(x0);
end
%(6)模拟退火结束后的可视化表述
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,'r-');
xlabel("迭代次数");
ylabel("y值");

② Obj_fun1.m

function y = Obj_fun1(x)
	y = 11*sin(x) + 7*cos(5*x);
end

结果：

最佳的位置是：<