模拟退火算法实战
案例一：求解y=x^(2) cos(x)+x^(2) cos(x) sin(x)
1、解的编码
由于该解的表示简单，采用实数编码即可。
2、确定初始解
根据定义域范围，随机生成初始解
#[start,end]为定义域
def initialization(start, end):
return random.uniform(start, end)
3、邻域新解的生成
① 判断生成的解是否在定义域范围内
def in_range(x, start, end):
return True if start <= x <= end else False
② 生成邻域新解
def generate_new(x, start, end):
while True:
#采用高斯分布生成新解
upper_bound = end - x
lower_bound = start - x
sigma = max(upper_bound, lower_bound) / 3
new_x = random.gauss(x, sigma)
#判断是否在定义域内，在则返回；否则重复生成
if in_range(new_x, start, end):
return new_x
4、确定能量函数
能量越低，被接受的概率越高。因此，当求解最大值目标函数时，可以以目标函数的倒数作为能量函数；当求解最小值目标函数时，可以直接以目标函数作为能量函数
def e(x):
return pow(x, 2) * math.cos(x)+ pow(x, 2) * math.cos(x)*math.sin(x)
5、新解接受准则
新解接受准则的三条原则
① 在固定温度下，接受使目标函数下降的候选解的概率要大于使目标函数上升的候选解概率；
② 随着温度的下降，接受使目标函数上升的解的概率要逐渐减小；
③ 当温度趋于零时，只能接受目标函数下降的解。
Metropolis准则
新解接受常采用标准的Metropolis准则
若，接受新解作为当前解；否则，按照概率判断是否接受新解。
def metropolis(e, new_e, t):
if new_e <= e:
return True
else:
p = math.exp((e - new_e) / t)
return True if random.random() < p else False
6、搜索最优解
def search(all_x):
best_e = 0xffff
best_x = all_x[0]
for x in all_x:
ex = e(x)
if ex < best_e:
best_e = ex
best_x = x
return best_x
7、主函数
#t0为初始温度，t_final为终止温度，alpha为冷却系数，inner_iter为内层迭代次数，[start,end]为定义域
def sa(t0, t_final, alpha, inner_iter, start, end):
#记录每次内层迭代生成的解
all_x = []
#生成初始解并加入解集
init = initialization(start, end)
all_x.append(init)
#设置初始温度
t = t0
#外层循环
while t > t_final:
x = search(all_x)
ex = e(x)
#内层循环
for i in range(inner_iter):
new_x = generate_new(x, start, end)
new_ex = e(new_x)
if metropolis(ex, new_ex, t):
x = new_x
ex = new_ex
all_x.append(x)
t = alpha * t
return all_x
8、问题求解
all_x = sa(3000, 10, 0.98, 200, -1.57, 20.18)
print(round(search(all_x), 3), round(e(search(all_x)), 3))
iteration = len(all_x)
all_x = np.array(all_x)
all_ex = all_x ** 3 * np.cos(all_x)
plt.xlabel("Iteration")
plt.ylabel("f(x)")
plt.plot(range(iteration), all_ex)
plt.show()
输出结果如下
模拟退火算法求解TSP模型
以下图为例，求从0号城市出发访问每一座城市并回到0号城市的最短回路
1、基本模型
① 解的编码
求解TSP模型采用自然数编码，0节点代表初始城市，1-7节点表示其他7座城市。如 1，7，6，3，2，5，4 表示从0出发依次经过城市1、7、6、3、2、5、4最后回到0
② 确定初始解
随机生成一个由数字1-7各出现一次的数列。
③ 邻域新解的生成
采用两点逆转的方法，即对数列设置两个点，然后数列在这两点中发生逆转，改变数列。逆转操作如下：假设数列为A=1,4,5,2,3,6,7，随机生成两个点2和6，那么A1=4,5,2,3,6，逆转后A1=6,3,2,5,4，A=1,6,3,2,5,4,7。
④ 确定能量函数
TSP模型以最短路径为目标函数，因此直接以路径长度为能量函数即可。
⑤ 新解接受准则
采用标准的Metropolis准则。
⑥ 搜索最优解
根据已生成的所有解的能量，搜索能量最低的解，即为模拟退火算法的最优解
2、代码求解
① 确定初始解
#numbers为需要经过的其他城市的数量
def initialization(numbers):
path_random = np.random.choice(list(range(1, numbers + 1)), replace=False, size=numbers)
path_random = path_random.tolist()
return path_random
② 邻域新解的生成
def generate_new(path):
numbers = len(path)
#随机生成两个不重复的点
positions = np.random.choice(list(range(numbers)), replace=False, size=2)
lower_position = min(positions[0], positions[1])
upper_position = max(positions[0], positions[1])
#将数列中段逆转
mid_reversed = path[lower_position:upper_position + 1]
mid_reversed.reverse()
#拼接生成新的数列
new_path = path[:lower_position]
new_path.extend(mid_reversed)
new_path.extend(path[upper_position + 1:])
return new_path
③ 确定能量函数
#length_mat为各个节点之间的最短距离矩阵
def e(path, length_mat):
numbers = len(path) - 1
length = length_mat[0][path[0]] + length_mat[path[-1]][0]
for i in range(numbers):
length += length_mat[path[i]][path[i + 1]]
return length
④ Metropolis准则
def metropolis(e, new_e, t):
if new_e <= e:
return True
else:
p = math.exp((e - new_e) / t)
return True if random.random() < p else False
⑤ 搜索最优解
def search(all_path, length_mat):
best_e = 0xffff
best_path = all_path[0]
for path in all_path:
ex = e(path, length_mat)
if ex < best_e:
best_e = ex
best_path = path
return best_path
⑥ 模拟退火算法主函数
#t0为初始温度，t_final为终止温度，alpha为冷却系数，inner_iter为内层迭代次数
#city_numbers为所需要经过的其他城市，length_mat为各个节点之间的最短距离矩阵
def sa(t0, t_final, alpha, inner_iter, city_numbers, length_mat):
all_path = []
all_ex = []
init = initialization(city_numbers)
all_path.append(init)
all_ex.append(e(init, length_mat))
t = t0
while t > t_final:
path = search(all_path, length_mat)
ex = e(path, length_mat)
for i in range(inner_iter):
new_path = generate_new(path)
new_ex = e(new_path, length_mat)
if metropolis(ex, new_ex, t):
path = new_path
ex = new_ex
all_path.append(path)
all_ex.append(ex)
t = alpha * t
return all_path, all_ex
⑦ 主程序
先通过Floyd算法求每个节点的最短路径
graph = nx.Graph()
graph.add_nodes_from(range(0, 8))
edges = [(0, 1, 3), (0, 4, 2),
(1, 4, 3), (1, 5, 4), (1, 7, 2),
(2, 3, 5), (2, 4, 1), (2, 5, 4), (2, 6, 3),
(3, 6, 1),
(4, 5, 5),
(5, 6, 2), (5, 7, 4),
(6, 7, 4)]
graph.add_weighted_edges_from(edges)
shortest_length = dict(nx.shortest_path_length(graph, weight='weight'))
length_mat = []
for i in range(0, 8):
i_j = []
for j in range(0, 8):
i_j.append(shortest_length[i][j])
length_mat.append(i_j)
再调用模拟退火算法求最优解，并可视化
all_path, all_ex = sa(3000, pow(10, -1), 0.98, 200, 7, length_mat)
print(search(all_path, length_mat), round(e(search(all_path, length_mat), length_mat)))
iteration = len(all_path)
all_path = np.array(all_path)
all_ex = np.array(all_ex)
plt.xlabel("Iteration")
plt.ylabel("Length")
plt.plot(range(iteration), all_ex)
plt.show()
3、输出结果
[1, 7, 5, 6, 3, 2, 4] 19
————————————————
参考链接：https://blog.csdn.net/weixin_58427214/article/details/125901431
matlab代码示例：
SA_1
clear
close all
clc
D = 10;
X_Min = -20;
X_Max = 20;
L = 500; % 每一个温度下的迭代次数
alpha = 0.998; % T_Next = T_Current*h
S0 = 10; % 变异算子
T0 = 10000; % 初始温度
T = T0;
error_Thre = 1e-8; % 代价函数阈值
T_Thre = 0.001;
K = 0.6; % Metropolis准则系数
% 迭代次数初始化
iter = 1;
iter_Max = log(T_Thre/T0)/log(alpha);
% iter_Max = T0/T_Thre;
% 参数初始化
Current_X = rand(1, D) .* (X_Max - X_Min) + X_Min;
Current_X_Fitness = Fitness(Current_X);
Pre_X = Current_X;
Pre_X_Fitness = Current_X_Fitness;
Best_X = Current_X;
MinFitness_List = zeros(1,floor(iter_Max)+100);
MinFitness_List(1) = Fitness(Best_X);
Fitness_List = zeros(1,floor(iter_Max)+100);
Fitness_List(1) = MinFitness_List(1);
% 与进度条有关的东西
waitbarinter = iter_Max/100;
tic
h = waitbar(0, ['已完成:0% 运算中...用时:', num2str(toc)]);
while (MinFitness_List(iter) > error_Thre) && (T > T_Thre)
iter = iter+1;
% 1.简易模拟退火
T = alpha * T;
% % 2. 经典模拟退火
% T = T0/log10(iter);
% % 3. 快速模拟退火
% T = T0/(iter);
% 在当前温度T下的迭代
% 当前温度下最优解的初始化
Fitness_temp = 0;
X_temp = zeros(1,D);
for i = 1:L
% 随机生成新的解
N_temp = randn(1,D);
P_temp = N_temp/sqrt(sum(N_temp.^2));
Current_X = Pre_X + T*P_temp;
% 边界限制 防止解过大或过小
Current_X = BoundaryLimit(Current_X, X_Min, X_Max);
Current_X_Fitness = Fitness(Current_X);
% 更新当前温度下的最优解
if i == 1
Fitness_temp = Current_X_Fitness;
X_temp = Current_X;
else
if Current_X_Fitness < Fitness_temp
Fitness_temp = Current_X_Fitness;
X_temp = Current_X;
end
end
% 判断当前解与最优解以及上一解的大小关系
if Pre_X_Fitness <= Current_X_Fitness
% 当前解不如上一次解 Metropolis准则，减小局部最优解的可能
P = exp(-1 * (Current_X_Fitness - Pre_X_Fitness) / K / T);
if P > rand
Pre_X = Current_X;
Pre_X_Fitness = Current_X_Fitness;
end
% 当前解比上一次解好
else
Pre_X = Current_X;
Pre_X_Fitness = Current_X_Fitness;
end
end
% 更新全局最优解
if MinFitness_List(iter-1) > Fitness_temp
MinFitness_List(iter) = Fitness_temp;
Best_X = X_temp;
else
MinFitness_List(iter) = MinFitness_List(iter-1);
end
% 将每个温度下的最优解都存储起来
Fitness_List(iter) = Fitness_temp;
% 进度条
if mod(iter, waitbarinter) <1
waitbar(iter / iter_Max, h, ['已完成:' num2str(iter / iter_Max * 100) ...
'% 运算中...用时:', num2str(toc),'/',num2str(toc/iter * iter_Max)])
end
end
close(h)
% 结果文字展示
display(['The best solution obtained by SA is : ', num2str(Best_X)]);
display(['The best optimal value of the objective funciton found by SA is : ', num2str(MinFitness_List(iter))]);
% 结果图像展示
figure
semilogy(MinFitness_List(1:iter),'linewidth',1.2)
hold on
semilogy(Fitness_List(1:iter),'linewidth',1.2)
legend(['最小适应度变化'],['适应度变化'])
title('Convergence Curve')
xlabel('Iteration');
ylabel('Best score obtained so far');
axis tight
grid on
box on
% 目标函数
function result = Fitness(x)
result = sum(x.^2);
end
% 边界限制函数
function result = BoundaryLimit(X, Min, Max)
for i_temp = 1:length(X)
if X(i_temp) > Max
X(i_temp) = Max;
end
if X(i_temp) < Min
X(i_temp) = Min;
end
end
result = X;
end
SA_TSP
clc,clear
sj=[53.7121 15.3046 51.1758 0.0322 46.3253 28.2753 30.3313 6.9348
56.5432 21.4188 10.8198 16.2529 22.7891 23.1045 10.1584 12.4819
20.1050 15.4562 1.9451 0.2057 26.4951 22.1221 31.4847 8.9640
26.2418 18.1760 44.0356 13.5401 28.9836 25.9879 38.4722 20.1731
28.2694 29.0011 32.1910 5.8699 36.4863 29.7284 0.9718 28.1477
8.9586 24.6635 16.5618 23.6143 10.5597 15.1178 50.2111 10.2944
8.1519 9.5325 22.1075 18.5569 0.1215 18.8726 48.2077 16.8889
31.9499 17.6309 0.7732 0.4656 47.4134 23.7783 41.8671 3.5667
43.5474 3.9061 53.3524 26.7256 30.8165 13.4595 27.7133 5.0706
23.9222 7.6306 51.9612 22.8511 12.7938 15.7307 4.9568 8.3669
21.5051 24.0909 15.2548 27.2111 6.2070 5.1442 49.2430 16.7044
17.1168 20.0354 34.1688 22.7571 9.4402 3.9200 11.5812 14.5677
52.1181 0.4088 9.5559 11.4219 24.4509 6.5634 26.7213 28.5667
37.5848 16.8474 35.6619 9.9333 24.4654 3.1644 0.7775 6.9576
14.4703 13.6368 19.8660 15.1224 3.1616 4.2428 18.5245 14.3598
58.6849 27.1485 39.5168 16.9371 56.5089 13.7090 52.5211 15.7957
38.4300 8.4648 51.8181 23.0159 8.9983 23.6440 50.1156 23.7816
13.7909 1.9510 34.0574 23.3960 23.0624 8.4319 19.9857 5.7902
40.8801 14.2978 58.8289 14.5229 18.6635 6.7436 52.8423 27.2880
39.9494 29.5114 47.5099 24.0664 10.1121 27.2662 28.7812 27.6659
8.0831 27.6705 9.1556 14.1304 53.7989 0.2199 33.6490 0.3980
1.3496 16.8359 49.9816 6.0828 19.3635 17.6622 36.9545 23.0265
15.7320 19.5697 11.5118 17.3884 44.0398 16.2635 39.7139 28.4203
6.9909 23.1804 38.3392 19.9950 24.6543 19.6057 36.9980 24.3992
4.1591 3.1853 40.1400 20.3030 23.9876 9.4030 41.1084 27.7149];%将原始数据放入括号即可，不需要任何的改动，这里省略。
x=sj(:,1:2:8);x=x(:);
y=sj(:,2:2:8);y=y(:);
sj=[x y];%得到坐标
d1=[70,40];%起始位置
sj=[d1;sj;d1];%将其实位置加入路径中
sjb=sj;%备份
sj=sj*pi/180;%得到弧度
%距离矩阵d
d=zeros(102);
for i=1:101
for j=i+1:102
temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));
d(i,j)=6370*acos(temp);%距离矩阵
end
end
d=d+d';%对称矩阵
S0=[];Sum=inf;
rand('state',sum(clock));%开始计时
for j=1:1000
S=[1 1+randperm(100),102];%生产一个路径
temp=0;
for i=1:101
temp=temp+d(S(i),S(i+1));%计算当前的路径的距离
end
if temp<Sum %更新
S0=S;Sum=temp;
end
end
%2变换------这里使用的是2路变换，可以使用3路变换，两者取其一。
e=0.1^30;L=100000000;at=0.999;T=1;
%退火过程
for k=1:L
%产生新解
c=2+floor(100*rand(1,2));
c=sort(c);
c1=c(1);c2=c(2);%随机产生两个交换点
%计算代价函数变换值
df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-d(S0(c2),S0(c2+1));
%接受准则
if df<0
S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
Sum=Sum+df;%接受，并记录
elseif exp(-df/T)>rand(1)
S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
Sum=Sum+df;%以概率接受，并记录
end
T=T*at;
if T<e
break;%退出
end
end
% 输出巡航路径及路径长度
S0,Sum
plot(sjb(S0,1),sjb(S0,2),'-bp');%输出结果
319
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
