智能计算方法与实现1|禁忌搜索算法|Rastrigin函数|压力容器设计问题

智能计算也有人称之为“软计算”，人们受自然（生物界）规律的启迪，根据其原理，模仿求解问题的算法。从自然界得到启迪，模仿其结构进行发明创造，这就是仿生学。这是我们向自然界学习的一个方面。另一方面，我们还可以利用仿生原理进行设计（包括设计算法），这就是智能计算的思想。
这方面的内容很多，如人工神经网络技术、禁忌搜索、遗传算法、模拟退火算法和群集智能技术（蚁群、粒子群等）等。

在数学建模竞赛中，有许多用常规方法难以在规定时间内获得较好解的优化问题，例如大规模整数线性规划问题、约束优化问题、旅行商问题、车辆路径问题等
针对这些问题，我们在这一讲开始介绍能够有效求解上述问题的智能计算方法，主要考虑这些方法在求解优化问题中的应用，重点介绍五种常用算法的应用：
禁忌搜索算法
模拟退火算法
遗传算法
蚁群算法
粒子群算法

智能计算方法概述

优化问题与智能计算方法

优化问题是指在满足一定条件下，在众多方案或参数值中寻找最优方案或参数值，以使得某个或多个功能指标达到最优或使系统的某些性能指标达到最大值或最小值
优化方法是以数学为基础，用于求解各种优化问题的应用技术，通常我们通过对实际问题，进行合理假设，建立适合的优化模型求解。
在数学建模竞赛中的许多赛题，最终都可以转化为某个优化问题求解。

由实际问题得到的优化问题模型，一般都是比较复杂的问题可能有大量的决策变量，或约束条件。
面对这些大型的优化问题，尤其是组合优化问题，传统的优化方法（牛顿法、单纯形法等）需要遍历整个解空间，无法在有限时间内获得比较满意的解。
鉴于实际工程问题的复杂性、非线性、约束性以及建模困难等特点，要求我们研究更有效地求解优化问题的计算技术和方法。

受到人类智能、生物群体社会性或自然现象规律的启发人们发明了许多智能优化方法来解决上述复杂优化问题。主要包括：
模拟人类记忆过程的禁忌搜索算法；
源于固体物质退火过程的模拟退火算法
模仿自然界生物进化机制的遗传算法；
模拟蚂蚁集体寻径行为的蚁群算法；
模拟鸟群群体行为的粒子群算法等

禁忌搜索算法简介

禁忌搜索算法TS（TabuSearch）是由美国科罗拉多州大学的 FredGlover教授在1986年左右提出来的，是一个用来跳出局部最优的搜寻方法。
禁忌搜索是一种亚启发式随机搜索算法，它从一个初始可行解出发，选择一系列的特定搜索方向（移动）作为试探，选择实现让特定的目标函数值变化最多的移动。为了避免陷入局部最优解，TS搜索中采用了一种灵活的记忆技术，对已经进行优化过程进行记录和选择，指导下一步的搜索方向

TS是人工智能的一种体现，是局部领域搜索的一种扩展
禁忌搜索是在领域搜索的基础上，通过设置禁忌表来禁忌一些已经历的操作，并利用视准则来奖励一些优良状态，其中涉及邻域、禁忌表、禁忌长度、候选解、蔑视准则等影响禁忌搜索算法性能的关键因素
迄今为止，TS算法在组合优化等计算机领域取得了很大的成功，近年来又在函数全局优化方面得到较多的研究，并大有发展的趋势。

模拟退火算法简介

金属退火的原理
金属退火是将金属加热到一定温度，保持足够时间，然后以适宜速度冷却（通常是缓慢冷却，有时是控制冷却）的一种金属热处理工艺。
模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态内能减为最小。

如上图，处在低温状态时，固体中分子具有的内能很低，在原本的位置上做小范围的振动。
若是将固体加热到一定温度，分子内能将会增加，热运动加剧，分子排列的无序度增加。
此时再将温度缓缓降低在每个温度都达到平衡态（即准静态过程），分子具有的能量逐渐降低，最终回归到有序排列的状态，分子内能也跟着降到最低，
2. 模拟退火算法机制
模拟退火算法（SimulatedAnnealing，SA）最早的思想是由N. Metropolis等人于1953年提出。1983年，S.Kirkpatrick等成功地将退火思想引入到组合优化领域。它是基于Monte-Carlo 送代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性
![[Pasted image 20240827143707.png]]

遗传算法简介

遗传算法（GeneticAlgorithm，简称GA）起源于1975年，Holland教授对生物系统所进行的计算机模拟研究，是一种随机全局搜索优化方法，它模拟了自然选择和遗传中发生的复制、交叉（crossover）和变异（mutation）等现象，从任一初始种群（Population）出发，通过随机选择、交叉和变异操作，产生一群更适合环境的个体，使群体进化到搜索空间中越来越好的区域，这样一代一代不断繁衍进化，最后收敛到一群最适应环境的个体（lndividual），从而求得问题的优质解
![[Pasted image 20240827143854.png]]

蚁群算法简介

蚁群算法（AntClonyOptimization，Aco）是一种群智能算法，它是由一群无智能或有轻微智能的个体（Agent）通过相互协作而表现出智能行为，从而为求解复杂问题提供了一个新的可能性。
蚁群算法最早是由意大利学者ColorniA.，DorigoM.等于1991年提出。经过20多年的发展，蚁群算法在理论以及应用研究上已经得到巨大的进步

蚁群算法是一种仿生学算法，是由自然界中蚂蚁觅食的行为而启发的。在自然界中，蚂蚁觅食过程中，蚁群总能够按照寻找到一条从蚁巢和食物源的最优路径。下图显示了这样一个觅食的过程
![[Pasted image 20240827144135.png]]

粒子群算法简介

1995年，美国学者Kennedy和Eberhart共同提出了粒子群算法
其基本思想源于对鸟类群体行为进行建模与仿真的研究结果的启发。
粒子群优化算法（PSO：Particleswarmoptimization）是一种进化计算技术，源于对鸟群捕食的行为研究，它的核心思想是利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得问题的可行解，
PSO的优势：在于简单容易实现并且没有许多参数的调节。目前已被广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域。
![[Pasted image 20240827144719.png]]

禁忌搜索算法原理

TS（tabusearch）是LocalSearch（LS）的扩展，是一种全局逐步寻优的全局性邻域搜索算法
TS模仿人类的记忆功能，在搜索过程中标记已经找到的局部最优解及求解过程，并于之后的搜索中避开它们
算法通过禁忌策略实现记忆功能，通过破禁准则继承LS的强局部搜索能力。使得TS一方面具备高局部搜索能力，同时又能防正算法在优化中陷入局部最优

TS主要构成要素

评价函数(Evaluation Function)
邻域移动(Move Operator)
禁忌表（TabuTable）
邻居选择策略(Neighbor Selection Strategy)
破禁准则（AspirationCriterion，又称藐视准则）
停止规则(Stop Criterion)

评价函数（EvaluationFunction）：用来评价邻域中的
邻居、判断其优劣的衡量指标。
邻域移动（MoveOperator）：进行解转移的关键，又
称“算子”，影响整个算法的搜索速度
禁忌表（TabuTable）：记录被禁正的变化，以防出现搜索循环、陷
入局部最优。
对其设计中最关键的因素是禁忌对象（禁忌表限定的对象）和禁忌步长（对象的禁忌在多少次迭代后失效）。
禁忌表是禁忌搜索算法的核心，禁忌表的对象、步长及更新策略在很大程度上影响着搜索速度和解的质量。若禁忌对象不准确或者步长过小，算法避免陷入局部最优的能力会大打折扣；若禁忌表步长过大，搜索区域将会限制，好的解就可能被跳过
邻居选择策略(Neighbor Selection Strategy)：选择最佳邻域移动的规则。目前最广泛采用的是“最优解优先策略”及“第一个改进解优先策略”。
破禁准则（AspirationCriterion，视准则）：对于禁忌表的适度放松。当某个被禁忌的移动可得到优于未被禁忌的移动得到的最优邻域解和历史所得到的最优解时，算法应接受该移动，不受禁忌表的限制
停止规则（StopCriterion）：禁忌搜索中停止规则的设计多种多样，如最大送代数、算法运行时间、给定数目的送代内不能改进解或组合策略等等。

TS算法性能

通过一系列关于TSP的实验得到以下结论：
禁忌策略大大加强了算法的搜索能力
问题规模较小时，禁忌搜索能得到最优解；
问题规模较大时，禁忌搜索能在规定时间内输出满意解。禁忌对象的选择对算法效果存在较大影响

TS基本步骤

初始化
利用贪婪算法等局部搜索算法生成一个初始解，清空禁忌表，设置禁忌长度
邻域搜索产生候选解
根据步骤1产生初始解，通过搜索算子（searchoperators），如relocation、exchange、2-opt等，产生候选解（candidatesolution），并计算各个候选解的适应值（即解对应的目标函数值）
选择最好的候选解
从步骤2产生的所有候选解中选出适应值最好的候选解，将其与当前最好解（即搜索算法开始到现在找到的最好解）进行比较。如果优于当前最好解，那么就不考虑其是否被禁忌，用这个最好的候选解来更新当前最好解，并且作为下一个迭代的当前解，然后将对应的操作加入禁忌表；如果不优于当前最好解，就从所有候选解中选出不在禁忌状态下的最好解作为新的当前解，然后将对应操作加入禁忌表。
判断终止条件
若满足终止条件，则立即停止并输出当前最好解：否则继续搜索。一般终止条件为是否到达一定的迭代次数或者达到了一个时间限制

禁忌搜索算法涉及：
编码解码（Encoding and decoding)
搜索算子（searchoperators）
邻域(Neighborhood)
禁忌表（Tabu list）
禁忌长度（Tabutenure）
候选解（Candidatesolution）
藐视准则（Aspirationcriterion）等关键组成部分

TS算例-TSP问题

![[Pasted image 20240827150507.png]]

%%%%禁忌搜索算法解决TSP问题
%%初始化
clear;       %清除所有变量
close all;   %清图
clc;         %清屏
figure(1);   %过程图设置%%数据预处理
[C, N, D] = DataO; %C-城市坐标N-城市个数D-城市间距离矩阵%%设置禁忌表与禁忌长度Tabu = zeros(N);   %禁忌表
TabuL = round((N*(N-1)/2)^0.5); %禁忌长度%%初始解
S0 = randperm(N);  %随机产生初始解
bestsofar = S0;    %当前最佳解
BestL = Inf;       %当前最佳解距离%%禁忌搜索主循环 
p = 1;
Gmax = 100;         %最大迭代次数
ArrBestL = zeros(Gmax, 1); 
while p <= Gmax%   ALong(p) = fitnessFunction(D, S0);  %当前解适配值[BestCa, BestCaNum, CaNum] = genNeighbor(N, S0, D); %生成邻域解%选择最好解，更新禁忌表[Tabu, S0, BestL] = selectBest(N, CaNum, BestL, BestCa, Tabu, TabuL); if fitnessFunction(D, bestsofar) > BestLbestsofar = S0; %更新到目前为止的最优解 endArrBestL(p) = BestL; %记录第p次迭代得到的最优解 plotSolution(C, bestsofar, BestL, p); %绘制求解过程图p = p + 1; %下一次迭代 
end%%记录最终结果，绘制适应度变化曲线
BestShortcut = bestsofar;    %最佳路线
theMinDistance = BestL;      %最佳路线长度
figure(2);
plotFitCurve(ArrBestL);      %绘制适应度函数值的选代变化

function [BestCa, BestCaNum, CaNum] = genNeighbor(N, S0, D)Ca = 200;      %候选集的个数(全部领域解个数)
CaNum = zeros(Ca, N);   %候选解集合
i = 1;
A = zeros(Ca, 2);   %解中交换的城市矩阵%%求领域解中交换的城市矩阵 
while i <= CaM = N*rand(1,2;M = ceil(M); if M(1) ~= M(2)A(i, 1）= max(M(1),M(2));A(i, 2）= min(M(1),M(2));if i == 1 isa = 0; elsefor j = 1 : i-1if A(i, 1) == A(j,1) && A(i,2)==A(j,2)isa = 1; break;elseisa=0;end endend if ~ isa i = i+1;else endelse end 
endfunction[Tabu, S0, BestL] = selectBest(N, CaNum, BestL0, BestCa, Tabu0, Tabu)Tabu = TabuO; BestL = BestLO;BestCaNum = size(BestCa,1);%%候选解选择和禁忌表更新if BestCa(1,2) < BestL %满足藐视准则 BestL = BestCa(1,2);S0 = CaNum(BestCa(1,1),:); for m = 1:Nfor n = 1:Nif Tabu(m,n) ~= 0%更新禁忌表，表中元素禁忌长度-1 Tabu(m,n)=Tabu(m,n)-1;end end end%更新禁忌表，增加新的最优解入表Tabu(BestCa(1,3), BestCa(1,4)) = TabuL;else %否则for i = 1:BestCaNumif Tabu(BestCa(i,3),BestCa(i,4)) == 0S0 = CaNum(BestCa(i,1),:); for m = 1:Nfor n = 1:Nif Tabu(m, n) ~= 0%更新禁忌表，表中元素禁忌长度-1 Tabu(m,n) = Tabu(m,n)-1;end end end%更新禁忌表，增加新的最优解入表 Tabu(BestCa(i,3), BestCa(i,4)) = TabuL; break;end endend 
end

%%适应度函数
function F = fitnessFunction(D, s)DistanV = 0; n = size(s, 2); for i = 1:(n-1)DistanV = DistanV + D(s(i), s(i+1));endDistanV = DistanV + D(s(n), s(1)); F = DistanV;
end

禁忌搜索算法应用-无约束优化

Rastrigin函数

$f(x_{1},x_{2})=20+x_{1}^{2}+x_{2}^{2}-10(\cos 2\pi x_{1}+\cos 2\pi x_{2})$
$-5\le x_{i}\le 5,\ i=1,2$
![[Pasted image 20240827215629.png]]
![[Pasted image 20240827215645.png]]

$-1\le x_{i}\le 1,\ i=1,2$
![[Pasted image 20240827215718.png]]
![[Pasted image 20240827215728.png]]

%%%%禁忌搜索算法求解Rastrigin函数最小值
%%初始化
clear;     %清除所有变量
close all;    %清图
clc;            %清屏
plotRastrigin();   %绘制Rastrigin函数图像%%设置禁忌表与禁忌长度 
narNum = 2;
TabuL = randi([5 11], 1, 1);    %禁忌长度
Tabu = zeros(TabuL, narNum);    %禁忌表%%初始解
xu = 1；%解的上界 
xl = -1; %解的下界
x0 = rand(1, narNum) *（xu-xl) + xl;   %随机产生初始解
bestsofar.key=x0;     %当前最佳解
bestsofar.value = fitnessFunction(x0);
xnow.key = x0;  %当前解
xnow.value = bestsofar.value;
w = 1;     %自适应权重系数

%%适应度函数
function F = fitnessFunction(x)F = 20 + x(1).*x(1) + x(2).*x(2) - 10*c0s(2*pi.*x(1)) - 10*c0s(2*pi.*x(2)); 
end

%%禁忌搜索主循环 
p = 1;
Gmax = 20;
ArrBestF = zeros(Gmax,1);
while p <= Gmax w = w*0.998;candidate = genNeighbor(xnow, w, xu, xl); %生成邻域解，选出候选解%选择最好解，更新禁忌表[Tabu, xnow_new, bestsofar_new] = selectBest(candidate, xnow, bestsofar, Tabu, Tabu); xnow = xnow_new;bestsofar = bestsofar_new; %更新到目前为止的最优解ArrBestF(p) = bestsofar.value; %记录第p次迭代得到的最优解 plotSolution(bestsofar, p); %绘制求解过程图p = p+1;%下一次迭代 
end

function[candidate] = genNeighbor(xnow, w, xu, xl) Ca = 5;x_near = zeros(Ca, 2);x_near_value = zeros(Ca, 1);
%%产生邻域解 for i = 1:Cax_temp = xnow.key;x1 = x_temp(1); x2 = x_temp(2);x_near(i, 1) = x1 + (2*rand-1)*w*(xu-xl);%边界条件处理 - 边界吸收if x_near(i,1) < xl x_near(i,1) = xl; endif x_near(i,1) > xu x_near(i,1) = xu; endx_near(i,2) = x2 + (2*rand-1)*w*(xu-xl);%边界条件处理-边界吸收if x_near(i,2) < xl x_near(i,2) = xl; endif x_near(i,2) > xu x_near(i,2) = xu; end%计算邻域解点的函数值x_near_value(i) = fitnessFunction(x_near(i, :));end%%最优邻域解为候选解temp = find(x_near_value == min(x_near_value)); temp = temp(1);candidate.key = x_near(temp, :);candidate.value = x_near_value(temp); end

function[Tabu, xnow_new, bestsofar_new] = selectBest(candidate, xnow, bestsofar, TabuO, Tabu) Tabu = Tabu0;
bestsofar_new = bestsofar;
%%候选解和当前解的评价函数差 
delta1 = candidate.value - xnow.value;%%候选解和目前最优解的评价函数差 
delta2 = candidate.value - bestsofar.value;%%候选解并没有改进解，把候选解赋给下一次迭代的当前解 
if delta1 >= 0xnow_new.key = candidate.key;xnow_new.value = candidate.value;%更新禁忌表Tabu = [Tabu; xnow_new.key]; if size(Tabu, 1) > TabuLTabu(1,:) = []; end%%如果相对于当前解有改进，则应与目前最优解比较 elseif delta2 < 0 %候选解比目前最优解优，把改进解赋给下一次迭代的当前解 xnow_new.key = candidate.key;xnow_new.value = candidate.value; %更新禁忌表Tabu = [Tabu;xnow_new.key]; if size(Tabu,1) > TabuLTabu(1,:) = []; end%把改进解赋给下一次选代的目前最优解%包含藐视准则bestsofar_new.key = candidate.key;bestsofar_new.value = candidate.value; else%判断改进解时候在禁忌表里[M,~] = size(Tabu); r = 0;for m = 1:Mif candidate.key(1) == Tabu(m,1)..&candidate.key(2) == Tabu(m,1)r=1; end end
if r == 0%改进解不在禁忌表里，把改进解赋给下一次迭代的当前解xnow_new.key = candidate.key;xnow_new.value = candidate.value;%更新禁忌表Tabu = [Tabu;xnow_new.key]; if size(Tabu,1) > TabuLTabu(1,:) = []; endelse%%%如果改进解在禁忌表里，用当前解重新产生邻域解xnow_new = xnow;end endend 
end

%%记录最终结果，绘制适应度变化曲线
% bestsofar 最优解
% BestF最优函数值 
figure(3);
%绘制适应度函数值的选代变化 
plotFitCurve(ArrBestF);

![[Pasted image 20240827222143.png]]

禁忌搜索算法应用-约束优化

压力容器设计问题

![[Pasted image 20240827222317.png]]

目标：使压力容器制作（配对、成型和焊接）成本最低。
压力器两端都由封盖封住，头部一端的封盖为半球状。L是不考虑头部的圆柱体部分的截面长度，R是圆柱体的内壁半径，T和T，分别表示圆柱体的壁厚和头部的壁厚，这四个变量是压力容器设计问题的4个优化变量。对应的自标和约束如下：
$min\ f(x)=0.6224x_{1}x_{3}x_{4}+1.7781x_{2}x_{3}^{2}+3.1661x_{1}^{2}x_{4}+19.84x_{1}^{2}x_{3}$
$s.t.\qquad g_{1}(x)=-x_{1}+0.0193x_{3}\le 0$
$g_{2}(x)=-x_{2}+0.00954x_{3}\le 0$
$g_{3}(x)=-\pi x_{3}^{2}-4\pi \frac{x_{3}^{3}}{3}+129600\le 0$
$g_{4}(x)=x_{4}-240\le 0$
$x=[x_{1},x_{2},x_{3},x_{4}]=[T_{s},T_{h},R,L];$
$100\ge x_{i}\ge 0,i=1,2$
$100\ge x_{i}\ge 10,i=3,4$

%%%%禁忌搜索算法求解压力容器设计问题
%%初始化
clear;     %清除所有变量
close all;   %清图
clc;     %清屏%%设置禁忌表与禁忌长度 
narNum = 4;
TabuL = randi([5 11], 1, 1);   %禁忌长度
Tabu = zeros(TabuL, narNum);   %禁忌表%%初始解
xl = [0, 0, 10, 10];  %解的下界
xu = [100, 100, 100, 100];  %解的上界 
x0 = rand(1, narNum).*(xu = xl) + xl;  %随机产生初始解
bestsofar.key = x0;  %当前最佳解
bestsofar.value = fitnessFunction(x0);
xnow.key = x0;  %当前解
xnow.value = bestsofar.value;
w = 1;    %自适应权重系数

%%适应度函数
function F = fitnessFunction(x)x1 = x(1); %Ts x2 = x(2); %Th x3 = x(3); %R x4 = x(4); %L%%约束条件判断 g1 = -x1 + 0.0193*x3;g2 = -x2 + 0.00954*x3;g3 = -pi*x3^2 - 4*pi*x3^3/3 + 1296000； g4 = x4 - 240;
%如果满足约束条件则计算适应度值if(g1 <= 0 && g2 <= 0 && g3 <= 0 && g4 <= 0)F = 0.6224*x1*x3*x4 + 1.7781*x2*x3^2 + 3.1661*x1^2*x4 + 19.84*x1^2*x3;else %否则适应度值无效 F = inf;end 
end

%%禁忌搜索主循环 
p = 1;
Gmax = 200;
ArrBestF = zeros(Gmax, 1);
while p <= Gmax w = w*0.998;candidate = genNeighbor（xnow, w, xu, xl); %生成邻域解，选出候选解[Tabu, xnow_new, bestsofar_new] = selectBest(candidate, xnow, bestsofar, Tabu, Tabu);
%选择最好解，更新禁忌表 xnow = xnow_new;bestsofar = bestsofar_new; %更新到目前为止的最优解ArrBestF(p) = bestsofar.value; %记录第p次迭代得到的最优解p = p+1；%下一次迭代 
end

![[Pasted image 20240827224633.png]]