CN103294929B - 一种自适应的直接搜索模拟退火方法 - Google Patents

Abstract

本发明公开了一种自适应的直接搜索模拟退火（ADSA）算法，从初始温度计算、迭代终止条件和新工作点产生机理等三方面对标准直接搜索模拟退火（DSA）算法进行了改进。初始温度结合问题规模和难易程度按自适应方法计算并增加最大迭代次数终止条件，以提高算法运行效率；通过对工作点集中最佳工作点变化的监测，引入新工作点搜索范围随迭代温度动态变化的自适应方法，并将标准DSA算法中的新工作点按均匀分布随机产生方式改为按变尺度柯西分布随机产生，以增强算法跳出局部极值能力。仿真结果表明，该算法能明显改进收敛解的精度和搜索效率，有效提高算法的稳定性和成功率。

Description

一种自适应的直接搜索模拟退火方法

技术领域

本发明涉及一种自适应的直接搜索模拟退火算法。

背景技术

文献1“Kirkpatrick S,Gelatt C D,Vecchi M P.Optimization by simulatedannealing.Science 1983,220(4598):671-680.”公开了一种模拟退火(SA)算法。这种算法不具备记忆功能。

文献2“Ali M M,A,Viitanen S.A direct search variant of the simulatedannealing algorithm for optimization involving continuousvariables.Computers&Operations Research.2002,29:87-102.”公开了一种直接搜索模拟退火(Direct Search Simulated Annealing,以下简称为“标准DSA”)算法，这种算法设有一个工作点集，按随机概率选择基于工作点集以受控方式产生新工作点或基于整个搜索空间按均匀分布随机方式产生新工作点，新工作点按概率取代工作点集中的最差工作点，具有记忆功能。

文献3“申玮，模拟退火算法的改进及其在水环境工程中的应用，硕士学位论文，浙江大学，2004.”公开了一种每次迭代按随机或正交设计方式生成包含当前工作点的工作点集、在工作点集中的最佳工作点邻域内随机产生新工作点的DSA算法。这种算法需动态生成工作点集，具有记忆功能。

文献4“Erdal O,Sonmez F O.Optimum design of composite laminates formaximum buckling load capacity using simulated annealing.Computers&Structures.2005,71:45-52.”公开了一种按随机概率选择基于工作点集以受控方式产生新工作点或在工作点集中随机选取一工作点进行随机扰动产生新工作点的DSA算法。这种算法没考虑在整个搜索空间以随机方式产生新工作点，具有记忆功能。

文献5“Akbulut M,Sonmez F O.Optimum design of composite laminates forminimum thickness.Computers&Structures.2008,86:1974-1982.”公开了一种在工作点集中随机选取一工作点进行随机扰动产生新工作点的DSA算法。这种算法没考虑基于工作点集以受控方式或在整个搜索空间以随机方式的新工作点产生方式，具有记忆功能。

文献6“Akbulut M,Sonmez F O.Design optimization of laminatedcomposites using a new variant ofsimulated annealing.Computers&Structures.2011,89:1712-1724.”公开了一种在工作点集中随机选取一工作点进行随机扰动产生新工作点、新工作点按概率随机取代工作点集中最差的n+1个工作点之一的DSA算法(n为设计变量个数)。这种算法提出了新的温度下降系数、新工作点取代机制，具有记忆功能。

以上文献的典型特点是：1)DSA算法因存在工作点集而具有记忆功能；2)标准DSA算法具有较好的全局和局部协调的搜索能力，但对高维优化问题，算法的稳定性差、成功率低，且易陷入局部极值；3)改进的DSA算法发展了基于工作点集的新工作点产生机制，对算法求解高维优化问题的性能表现有一定改善，但因舍弃在整个搜索空间的随机产生新工作点方式，弱化了算法的全局搜索能力，对算法易陷入局部极值问题没有改进。

发明内容

本发明的目的在于提供了一种自适应的直接搜索模拟退火(ADSA)算法，它实现了算法跳出局部极值进行全局搜索、加快搜索效率、提高算法稳定性和成功率的要求。

本发明是这样来实现的，改变标准DSA算法的初始温度计算方法、迭代终止条件设置及新工作点产生方式，即初始温度计算结合问题规模和难易程度按自适应方法处理；迭代终止条件增加最大迭代次数限制；通过对工作点集中最佳工作点变化的监测，引入新工作点搜索范围随迭代温度动态变化的自适应方法，并将标准DSA算法的新工作点按均匀分布随机产生方式改为按变尺度柯西分布随机产生。它包括下述步骤：

(a)给定算法的收敛精度参数ε₁、ε₂，概率参数q，温度下降系数α_max、α_min，最大迭代次数V_max，函数值变化参数Δf，区间参数[e,f]。

(b)在搜索空间S内按均匀分布随机产生20n个工作点(n为设计变量个数，每个工作点对应一个设计变量向量X＝(x₁，x₂，…，x_n)^T)，如果20n个工作点中有可行解，则初始温度T₀＝500n；否则，初始温度T₀＝2000n。从20n个工作点中选出N＝7(n+1)个评价函数值最小的工作点储存在工作点集A中，并标记A中的最大和最小评价函数值为F_max和F_min，计算dF＝F_max-F_min，初始化Flag1(当A中无可行解时，Flag1＝1；当A中既有可行解又有不可行解时，Flag1＝2；当A中全是可行解时，Flag1＝3)。

(c)初始化迭代次数k＝0、A中最佳工作点未显著更新的迭代次数Δk＝0，确定初始温度下降系数α₀＝(α_min+α_max)/2和Markov链长度L₀＝10n+int(10n(1-e^-dF))，其中int(*)为取整函数。

(d)若k＜1，进入步骤(f)；否则，进入步骤(e)。

(e)k＝k+1；计算Markov链长度L_k、温度T_k：

T_k＝a_kT_k-1

(2)

其中：

为温度T_k-1时A中最大与最小评价函数值；L′_k-1为温度T_k-1时评价函数的实际计算次数。

(f)进入子循环，初始化循环次数j＝0；初始化A中最佳工作点更新记录FLAG＝0(FLAG＝0表示A中最佳工作点未更新，FLAG＝1表示A中最佳工作点已更新)；初始化Flag1变更记录AF＝0(当Flag1＝1变更为Flag1＝2时，AF＝1；其它情况，AF＝0)。

(g)j＝j+1；在设定区间[e,f]内产生一个随机整数

(h)若则在(0,1)内产生一个随机数w。当w≥q时，以式(3)按受控方式产生新工作点X；否则，以式(4)在整个搜索空间S内按变尺度柯西分布产生新工作点X。

X＝2G-X_n+1

(3)

其中：

X_i(i＝1,…,n+1)由A中记为X₁的最佳工作点与A中余下N-1个工作点中随机获取的n个工作点组成；X^U和X^L分别为X的上下限；符号“ο”表示各向量对应位置分量依次相乘；abs(*)为取其绝对值函数；

r为在区间(0,1)内按均匀分布产生的一个随机数；y是以位置参数为0、尺度参数为的柯西分布随机产生的n维单位列向量，按下式计算：

若按式(5)产生新工作点X：

X＝BestX+bIο(X^U-X^L)οy (5)

其中：

BestX为A中最佳工作点；

y是按标准正态分布随机产生的n维单位列向量；m_k为在区间(0,5)内按均匀分布随机产生的一个随机数；α_i为在区间[0,1]内按均匀分布产生一个随机数；β为自适应概率。

(i)若X∈S，进入步骤(j)；否则，按(6)式修正X的各分量：

其中：

x_i为X的第i个分量，和分别为分量x_i的上下限；mod(*)为取余数函数。

(j)计算X的评价函数值F(X)，按如下规则更新FLAG及计算X被接受的概率P_a：

(j.1)当Flag1＝1时，若X是可行解，则FLAG＝1，P_a＝1；若X是不可行解，则分别按式(7)和(8)计算FLAG和P_a的值。

(j.2)当Flag1＝2时，若X是可行解，则FLAG按式(7)计算，P_a＝1；若X是不可行解，则FLAG＝0，P_a按式(8)计算。

(j.3)当Flag1＝3时，若X是可行解，则分别按式(7)和(8)计算FLAG和P_a的值；若X是不可行解，则FLAG＝0，P_a＝0。

(k)在[0,1)内随机产生一个随机数P_r，若P_a＞P_r，则接受X并用X取代A中最差工作点；否则，直接进入步骤(m)。

(m)若FLAG＝1或j＝L_k,更新A中工作点的评价函数值及F_max、F_min、dF、Flag1，更新AF，按式(9)更新Δk，转入步骤(n)；否则，更新A中F_max，返回步骤(g)。

(n)迭代终止判断：如Flag1＞1且T_k≤ε₁,dF≤ε₂、或k≥V_max，则终止计算；否则，进入步骤(e)。

本发明将初始温度设置为与设计变量的个数及问题的难易度相关联，且利用计算初始温度时所产生的工作点来初始化工作点集A，可明显提高算法的通用性及计算效率。按步骤(h)设置的新工作点产生机理既可实现在算法运行前期进行大范围粗搜索，又可实现在后期进行局部精细搜索，而式(5)表示的自适应寻优方法还可以提前跳出局部极值点，最终使得算法的全局搜索能力、稳定性和搜索效率大大提高。按步骤(n)设置的终止条件可以很好的平衡算法运行时间和解的计算精度之间的矛盾。本发明的评价函数值计算可以是基于可行规则方法、惩罚函数法等多种有约束优化问题的约束处理方法。

本发明的技术效果是：本发明通过改变标准DSA算法的初始温度计算方法、迭代终止条件设置及新工作点产生方式，实现算法跳出局部极值进行全局搜索、加快搜索效率、提高算法稳定性和成功率的要求。

附图说明

图1是按可行规则法处理约束函数的自适应DSA算法流程图。

具体实施方式

实施例1：采用本发明的自适应DSA算法计算文献“罗辞勇,陈民铀,韩力,适应性粒子群寻优算法II,控制与决策,2009,24(6):859-863.”中的Rosenbrock函数(n＝30)：

s.t.-50≤x_i≤50,i＝1,2,…,n

该函数全局极小值为minf(X)＝f(1,1,…,1)＝0。

本实施例是一个无约束优化问题，该问题的评价函数F(X)和目标函数f(X)是一致的，工作点集A中的工作点均为可行解，因此Flag1＝3，且被接受的新工作点均在搜索空间内，也是可行解。

利用本发明的算法求解该实施例的主要步骤如下：如图1所示；

Step1：算法参数设定为：ε₁＝ε₂＝1e-4，q＝0.5，V_max＝5.5×10⁴，

α_min＝0.8，α_max＝0.99，β＝0.1，Δf＝0.1，e＝1，f＝5。

Step2：计算初始温度T₀＝500n＝15000(本例n＝30)并初始化工作点集A，标记A中的最大和最小评价函数值为F_max和F_min，计算dF＝F_max-F_min，Flag1＝3；

Step3：k＝0；Δk＝0；计算α₀＝0.895；计算Markov链长度L₀。

Step4：若k＜1，进入Step6；否则，进入Step5。

Step5：k＝k＋1；计算Markov链长度L_k和温度T_k。

Step6：进入子循环，j＝0；FLAG＝0；AF＝0。

Step7：j＝j＋1；在设定区间[1,5]内产生一个随机整数

Step8：若则在(0,1)内产生一个随机数w。当w≥0.5时，以式(3)按受控方式产生新工作点X；否则，以式(4)在整个搜索空间S内按变尺度柯西分布产生新工作点X。

若按式(5)产生新工作点X。

Step9：若X∈S，进入Step10；否则，按(6)式修正X的各分量。

Step10：计算X的评价函数值F(X)，按如下规则更新FLAG及计算X被接受的概率P_a：

当Flag1＝3时，X是可行解，按式(7)和(8)计算FLAG和P_a。

Step11：在[0,1)内随机产生一个随机数P_r，若P_a＞P_r，则接受X并用X取代A中最差工作点；否则，进入Step12。

Step12：若FLAG＝1或j＝L_k,更新A中工作点的评价函数值及F_max、F_min、dF，Flag1＝3，AF＝0，按式(9)更新Δk，转入Step13；否则，更新A中F_max，返回Step7。

Step13：迭代终止判断：如Flag1＞1且T_k≤ε₁,dF≤ε₂、或k≥V_max，则终止计算；否则，进入Step5。

采用本发明的自适应DSA算法对该实施例重复运行200次，统计的测试结果如表1所示。表中“APSO-Ⅱ算法”是文献“罗辞勇,陈民铀,韩力,适应性粒子群寻优算法II,控制与决策,2009,24(6):859-863.”的适应性粒子群寻优算法Ⅱ。

表1Rosenbrock函数(n＝30)测试结果

表1数据说明，本发明的自适应DSA算法能有效处理复杂的高维优化问题。

实施例2：

0≤x_i≤1(i＝1,2,…,n)

该函数全局极大值为本实施例取n＝30。

该实施例需要将最大化等式约束问题转换为最小化不等式问题，且采用文献“Kalyanmoy D.An efficient constraint handling method for geneticalgorithms.ComputerMethods inAppliedMechanics andEngineering.2000,186:311-338.”提出的可行规则法处理不等式约束。

首先，通过引入松弛系数δ＝10^-3将等式约束转换为两个不等式约束：

则转换为最小化不等式问题后的测试函数为：

0≤x_i≤1(i＝1,2,…,n)

基于可行规则方法，优化问题的评价函数F(X)按式(14)计算，

其中F为可行域；f(X)为目标函数；f_max为工作点集A中最差可行解(即目标函数值最大的可行解)的目标函数值，当A中无可行解时，设定f_max＝0；C(X)为不可行度，按式(15)计算：

l为不等式约束的个数，本实施例为2。当X为可行解时，C(X)＝0；否则C(X)＞0。

采用可行规则方法，在每次工作点集A更新后，计算A中工作点的目标函数值f(X)和不可行度C(X)，确定A中最差可行解的目标函数值f_max，然后按式(14)计算A中各工作点的评价函数值F(X)；并根据A中各工作点的C(X)值确定Flag1。

利用本发明的自适应DSA算法求解该实施例的主要步骤如下：

Step1：算法参数设定为：ε₁＝ε₂＝0.05，q＝0.5，V_max＝5.5×10⁴，α_min＝0.8，α_max＝0.99，β＝0.1，Δf＝0.1，e＝1，f＝5。

Step2：计算初始温度T₀，并按式(14)初始化工作点集A，初始化Flag1，标记A中的最大和最小评价函数值为F_max和F_min，计算dF＝F_max-F_min。

Step3：k＝0；Δk＝0；计算α₀＝0.895；计算Markov链长度L₀。

Step4：若k＜1，进入Step6；否则，进入Step5。

Step5：k＝k+1；计算Markov链长度L_k和温度T_k。

Step6：进入子循环，j＝0；FLAG＝0；AF＝0。

Step7：j＝j+1；在设定区间[1,5]内产生一个随机整数

Step8：若则在(0,1)内产生一个随机数w。当w≥0.5时，以式(3)按受控方式产生新工作点X；否则，以式(4)在整个搜索空间S内按变尺度柯西分布产生新工作点X。

若按式(5)产生新工作点X。

Step9：若X∈S，进入Step10；否则，按(6)式修正X的各分量。

Step10：计算X的评价函数值F(X)，按如下规则更新FLAG及计算X被接受的概率P_a：

(a)当Flag1＝1时，若X是可行解，则FLAG＝1，P_a＝1；若X是不可行解，则分别按式(7)和式(8)计算FLAG和P_a的值。

(b)当Flag1＝2时，若X是可行解，则FLAG按式(7)计算，P_a＝1；若X是不可行解，则FLAG＝0，P_a按式(8)计算。

(c)当Flag1＝3时，若X是可行解，则分别按式(7)和式(8)计算FLAG和P_a的值；若X是不可行解，则FLAG＝0，P_a＝0。

Step11：在[0,1)内随机产生一个随机数P_r，若P_a＞P_r，则接受X并用X取代A中最差工作点；否则，直接进入Step12。

Step12：若FLAG＝1或j＝L_k,按式(14)更新A中各工作点的评价函数值及F_max、F_min、dF、Flag1，更新AF，按式(9)更新Δk，转入Step13；否则，更新A中F_max，返回Step7。

Step13：迭代终止判断：如Flag1＞1且T_k≤ε₁,dF≤ε₂、或k≥V_max，则终止计算；否则，进入Step5。

采用标准DSA算法和本发明的自适应DSA算法分别对该实施例重复运行100次，统计的测试结果如表2所示。

表2有约束优化问题的测试结果

(注：表2中的数据为最小化结果；表中出现的最好解(-1.0081)好于理论极值(-1)的主要原因是将等式约束转换为不等式约束时的松弛系数δ＝10^-3引起的)

表2数据显示，与标准DSA算法相比，本发明自适应DSA算法的平均解精度提高了190.8％、评价函数平均计算次数降低了98.7％(即计算速度提高了98.7％)、平均解的标准差(即算法稳定性)提高了97.3％。验证了本发明的自适应DSA算法处理高维约束优化问题的有效性。

Claims (2)

1.一种自适应的直接搜索模拟退火方法，其特征在于：(1)改变初始温度的计算方法，即初始温度结合问题规模和难易程度按自适应方法计算，且工作点集A的初始化在计算初始温度阶段完成；(2)改变迭代终止条件，即在标准DSA算法基础上增加最大迭代次数的限制；(3)改变新工作点产生方式，即通过对工作点集中最佳工作点变化的监测，引入新工作点搜索范围随迭代温度动态变化的自适应方法，并将标准DSA算法的新工作点按均匀分布随机产生方式改为按变尺度柯西分布随机产生；所述方法包括如下步骤：

(1)给定相关参数；按自适应方法计算初始温度T₀，即在搜索空间S内随机产生一组工作点，如果该组工作点中有可行解，则T₀＝500n，其中，n为设计变量个数，否则，T₀＝2000n；从该组工作点中选取N＝7(n+1)个评价函数值最小的工作点储存在工作点集A中，计算A中最大和最小评价函数值之差dF，完成A的初始化；

(2)初始化迭代次数k＝0、A中最佳工作点未显著更新的迭代次数Δk＝0，计算初始温度下降系数α₀和Markov链长度L₀；

(3)若k＜1，进入步骤(5)，否则，进入步骤(4)；

(4)k＝k+1，计算Markov链长度L_k、温度T_k；

(5)进入子循环，初始化循环次数j＝0；

(6)j＝j+1，在设定区间内产生一个随机整数

(7)若则按概率以受控方式或在整个搜索空间S内按变尺度柯西分布产生新工作点；否则，在A中最佳工作点附近按搜索范围随迭代温度T_k动态变化的自适应方法产生新工作点，

若按下式产生新工作点X：

X＝BestX+bIο(X^U-X^L)οy

其中：X中每个分量X_i(i＝1，…，n+1)由A中记为X₁的最佳工作点与A中余下N-1个工作点中随机获取的n个工作点组成；X^U和X^L分别为X的上下限；符号“ο”表示各向量对应位置分量依次相乘；BestX为A中最佳工作点；

$b=\left\{\begin{array}{c}1e-4,T_k<1\\ \&lsqb;\arctan \left(T_k\right)-0.75\&rsqb;/m_k,T_k\&GreaterEqual;1and(\stackrel{\&OverBar;}{\mathrm{\&Delta;}K}\&le;\mathrm{\&Delta;}k<15or\mathrm{\&Delta;}k\&GreaterEqual;75)\\ \&lsqb;\arctan \left(T_k\right)-0.75\&rsqb;\&times;{1.1}^{\mathrm{\&Delta;}k-15}/m_k,T_k\&GreaterEqual;1and15\&le;\mathrm{\&Delta;}k<75\end{array}\right.$

$I_i=\left\{\begin{array}{c}1,a_i<\mathrm{\&beta;}\\ 0,a_i\&GreaterEqual;\mathrm{\&beta;}\end{array}\right.,(i=1,2,\mathrm{...},n)$

y是按标准正态分布随机产生的n维单位列向量；m_k为在区间(0,5)内按均匀分布随机产生的一个随机数；a_i为在区间[0,1]内按均匀分布产生一个随机数；β为自适应概率；

新工作点若在S之外，则修正在S之外的新工作点分量；

(8)计算新工作点的评价函数值，按概率接受新工作点并取代A中最差工作点；

(9)若新工作点优于A中最佳工作点或j＝L_k，则更新A中工作点的评价函数值及dF、Δk，进入步骤(10)；否则，返回步骤(6)；

(10)迭代终止判断：若A中含有可行解且dF和温度T_k达到设定精度，或迭代次数k超过设定的最大迭代次数，则终止计算；否则，返回步骤(4)。

2.如权利要求1所述的一种自适应的直接搜索模拟退火方法，其特征在于，所述的评价函数值的计算方法是基于可行性规则法和惩罚函数法之一的多种约束处理方法。