赞
踩
通过0-1背包问题看穷举法、贪心算法、启发式算法(JAVA)
java集成Cplex:Cplex下载、IDEA环境搭建、docker部署
windows下cplex12.6.3的下载、安装、IDE编程及相关问题解决
windows下cplex20.1.0启动报错: ilog.odms.ide.opllang.IloOplLangPlugin
模拟退火算法来源于固体退火原理,是一种基于概率的算法,将固体加温至充分高,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状,内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。
模拟退火算法(Simulated Annealing,SA)最早的思想是由N. Metropolis 等人于1953年提出。1983 年,S. Kirkpatrick 等成功地将退火思想引入到组合优化领域。它是基于Monte-Carlo迭代求解策略的一种随机寻优算法,其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发,伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解,即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是一种通用的优化算法,理论上算法具有概率的全局优化性能,目前已在工程中得到了广泛应用,诸如VLSI、生产调度、控制工程、机器学习、神经网络、信号处理等领域。
模拟退火算法有如下特点:
public static void main(String[] args) { // 定义初始温度、温度衰减系数 double initT=1000; double decreaseRate = 0.95; // 生成初始解initX, 及该解对应的函数值initY(评价函数可以直接使用目标函数) initX = generateInitSolve(); initY = targetFun(initX); // 初始解作为当前最优解 bestX = initX; bestY = initY; for(i=0; i<1000; i++){ // 外层for循环 每次循环温度降低 // 计算当前的迭代温度 currT = initT * decreaseRate^i; for(j=0; j<200; j++){ // 在当前温度下循环多次 // 在当前最优解附近产生一个新解,并计算新解的函数值 tmpX = getNewX(bestX, currT); tmpY = targetFun(tmpX); // 评价新解和当前最优解,决定取舍 if(tmpY>bestY){ // 新解更优,则接受新解为当前最优解 bestY = tmpY; bestX = tmpX; }else { // 新解更差,则以一定的概率接受新解为当前最优解 // 根据当前温度和两个解的差值生成接受概率 p = Math.exp(-(bestY - tmpY)/currentTemp); // 生成一个均匀分布的概率 randomP = random(0.0, 1.0); // 比较两个概率,决定是否接受 if(p > randomP){ bestY = tmpY; bestX = tmpX; } } } } // 当前最优解就是全局最优解 System.out.println("最优解:"+bestX, "最优解对应的函数值:"+bestX); }
退火是一种金属热处理工艺,先把金属加热到一个很高的温度,然后以适宜的速度冷却,冷却后金属就可以获得一些之前没有的特性。在算法中也是模拟了这个过程。算法开始时,给一个很高的初始温度initT,每次外层循环使温度以一个固定的衰减速率decreaseRate降低(通常在0.5-0.99之间);但是在每个温度下(currT = initT * decreaseRate^i),内层循环在同一个温度下多次寻找当前最优解附近的新解,并决定接受更好的新解,同时决定要不要接受一个更差的新解。
在启发式算法中,初始解的生成也是关键的一步,特别是对应爬山算法,如果初始解选择有问题就容易陷入局部最优的情况。因为模拟退火算法会以一定的概率接受更差的解,可以很大程度上避免陷入局部最优的情况,通常初始解的生成只要满足约束条件即可。
怎样在当前解的附近生成一个新解,也是模拟退火算法的关键一步。对应新解的生成,在学术上没有统一的规定,具体问题具体分析。
比如:
目标函数: m i n f ( X ) = − x 2 + 3 x + 6 ; 约束条件: x > = − 30 并且 x < = 30 目标函数:min f(X) = -x^2 + 3x + 6 ; 约束条件:x>= -30 并且x<=30 目标函数:minf(X)=−x2+3x+6;约束条件:x>=−30并且x<=30
当前解是x1,当前温度是t1,解的下限xMin=-30,解的上限xMax=30,r1是在区间[0,1]之间的服从N(0,1)分布的随机数,r2是在区间[0,1]之间的均匀分布的随机数。新解生成的方法可以简单写成:
public double getNewX(x1,t1,xMin,xMax,r1,r2){
newX = x1 + t1 * r1;
if(newX < xMin){
newX = r2 * xMin + (1 -r2) * x1;
}
if(newX > xMax){
newX = r2 * xMax + (1 -r2) * x1;
}
return newX;
}
这种方式参考Matlab的内置函数。
给定一组城市以及每个城市之间的距离,从一个城市出发,每个城市经过一次,回到出发的城市,走过的最小路径问题。产生新解的三种方法:
参考:《旅行商问题(TSP)的改进模拟退火算法》苗卉,杨韬 ,2007。
有15个书店,每个书店都有各种书籍,每个书店的运费是不一样的。现在需要买20本书,如何买花费最少。
产生新解的方法(参考):
参考《一种改进的模拟退火算法求解0-1背包问题》 梁国宏,张生,黄辉,何尚录 2007
模拟退火算法中,当遇到更差的解时,还是有一定的概率接受这个更差的解,而不是100%拒绝更差的解,从而克服了陷入局部最优的缺陷。
模拟退火算法的关键点就是如何定义接受一个更差解的概率。这个概率有几个比较重要点,1.这个概率函数的值在[0-1]之间;2.当前解的目标函数值f(A)与新解的目标函数值f(B)的差越大,说明新解和当前解差的越多,那么接受新解的概率就越小,反之越大;2.当前温度和概率的关系是,温度越高(搜索的前期),接受更差解的可能越大,这样可以在更大的范围搜索;温度越低(搜索的后期),接受更差解的可能性越小,这样就在小范围搜索。
根据上面的描述,这个概率函数p应该相关的几个部分参数 , 1.两个解的目标函数值之差 f =|f(B)-f(A)|, p与f成反比; 2. 当前温度Tc和时间t(迭代次数)成反比,通常每次迭代温度变成之前的95%,初始温度为Ti时,则 Tc = Ti * 0.95^t; 3.概率p和当前温度的关系是当前温度越高,概率越大。
根据上面的特点,这个概率函数可以构建成:
p
=
e
−
∣
f
(
A
)
−
f
(
B
)
∣
/
T
c
p=e^{-|f(A)-f(B)|/Tc}
p=e−∣f(A)−f(B)∣/Tc
退出搜索的条件可以是:1.达到指定的搜索次数,比如1000次;2.达到指定的温度,比如0.000001;3.找到的连续最优解M在多次(比如50次)迭代后都没有发生变化。
求解函数y = 11*sin(x) + 7*cos(5*x)在[-3,3]内的最大值
package com.wuxiaolong.algorithm.simulatedAnnealing; public class Constants { // 初始温度 public static final Double INIT_TEMPERATURE = 100.0; // 最大迭代次数 public static final Integer ITERATION_TIMES = 2000; // 每个温度下的迭代次数 public static final Integer CURRENT_TEMP_ITERATION_TIMES = 50; // 温度衰减系数 public static final Double TEMP_DECREASE_RATE = 0.95; // x的上、下界 public static final Double MIN_X = -3.0; public static final Double MAX_X = 3.0; }
package com.wuxiaolong.algorithm.simulatedAnnealing; import java.util.Random; public class Util { /** * 获取随机值 min <= 返回值x <= max * @return */ public static Double random(Double min,Double max) { Random rand = new Random(); double result=0; for(int i=0; i<10; i++){ result = min + (rand.nextDouble() * (max - min)); result = (double) Math.round(result * 100) / 100; // System.out.println(result); } return result; } /** * 随机数 服从N(0,1) * * 若随机变量X服从一个数学期望为μ、方差为σ^2的正态分布,记为X~N(μ,σ^2)。 * * 其概率密度函数为正态分布的期望值μ决定了其位置,其标准差σ决定了分布的幅度。当μ = 0,σ = 1时的正态分布是标准正态分布。 */ public static Double random(){ Random r = new Random(); return r.nextGaussian(); } }
package com.wuxiaolong.algorithm.simulatedAnnealing; import java.util.HashMap; import java.util.Map; import static com.wuxiaolong.algorithm.simulatedAnnealing.Constants.*; /** * 模拟退火算法 * // SA 模拟退火: 求解函数y = 11*sin(x) + 7*cos(5*x)在[-3,3]内的最大值 */ public class SA { public static void main(String[] args) { // 迭代中温度会发生改变,第一次迭代时温度就是T0 Double currentTemp = INIT_TEMPERATURE; // 随机生成一个初始解 并计算当前解的函数值 Double initX = Util.random(-3.0,3.0); Double initY = targetFun(initX); System.out.println("初始解:"+initX +" 初始函数值:"+initY); // 全局最优解 Double bestX = initX; Double bestY = initY; // 模拟退火过程 // 外循环, 我这里采用的是指定最大迭代次数 for(int iter=1; iter<=ITERATION_TIMES; iter++){ //System.out.println("-----------------"+iter+"------------------------"); // 本次迭代(当前温度)最优解 Double currTempBestX = bestX; Double currTempBestY = bestY; // 循环,在每个温度下开始迭代 for(int i=1; i<CURRENT_TEMP_ITERATION_TIMES; i++){ //找一个当前解附近的x, 计算新解的函数值 Double tmpX = getNewX(initX, currentTemp, MIN_X, MAX_X); Double tmpY = targetFun(tmpX); // 如果新解函数值大于当前解的函数值 更新当前解为新解 if(tmpY > currTempBestY) { currTempBestX = tmpX; currTempBestY = tmpY; }else{ //根据Metropolis准则计算一个概率 Double p = Math.exp(-(currTempBestY - tmpY)/currentTemp); // 生成一个随机数和这个概率比较,如果该随机数小于这个概率 Double random = Util.random(0.0, 1.0); if(random < p) { // 更新当前解为新解 currTempBestX = tmpX; currTempBestY = tmpY; } } } // 判断是否要更新找到的最佳的解 如果当前解更好,则对其进行更新 if(currTempBestY > bestY ) { bestY = currTempBestY; // 更新最大的y bestX = currTempBestX; // 更新找到的最好的x } System.out.println("iter:" + iter +" 新解:"+bestX +" 新函数值:"+bestY +" 当前温度"+currentTemp); // 温度下降 currentTemp = TEMP_DECREASE_RATE * currentTemp; } System.out.println("最终:"+bestX+"------------"+bestY); } public static Double targetFun(Double x){ return 11 * Math.sin(x) + 7*Math.cos(5*x); } public static Double getNewX(Double x0,Double T,Double x_lb,Double x_ub){ Double y =Util.random(); // 为变量生成N(0,1)随机数 Double x_new = x0 + y*T; // 根据新解的产生规则计算x_new的值 // 如果这个新解的位置超出了定义域,就对其进行调整 if (x_new < x_lb){ Double r = Util.random(0.0,1.0); x_new = r*x_lb + (1-r)*x0; }else if(x_new > x_ub){ Double r = Util.random(0.0,1.0); x_new = r*x_ub+(1-r)*x0; } // System.out.println("x0:"+ x0 + "---------------x_new:"+x_new); return x_new; } }
由于目前每次运行的结果还不一样,可能陷入了局部最优,还需要优化。
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。