|
模拟退火算法
1982年,KirkPatrick将退火思想引入组合优化领域，提出一种解大规模组合优化问题的算法，对NP完全组合优化问题尤其有效。这源于固体的退火过程，即先将温度加到很高,再缓慢降温(即退火)，使达到能量最低点。如果急速降温(即为淬火)则不能达到最低点.。
模拟退火算法是一种能应用到求最小值问题或基本先前的更新的学习过程（随机或决定性的）。在此过程中，每一步更新过程的长度都与相应的参数成正比，这些参数扮演着温度的角色。然后，与金属退火原理相类似，在开始阶段为了更快地最小化或学习，温度被升得很高，然后才（慢慢）降温以求稳定。
模拟退火算法是一种用于求解大规模优化问题的随机搜索算法，它以优化问题求解过程与物理系统退火过程之间的相似性为基础；优化的目标函数相当于金属的内能；优化问题的自变量组合状态空间相当于金属的内能状态空间；问题的求解过程就是找一个组合状态，使目标函数值最小。利用Metropolis准则并适当地控制温度的下降过程实现模拟退火，从而达到在多项式时间内求解全局优化问题的目标。
模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis准则，粒子在温度T时趋于平衡的概率为e-ΔE/(kT)，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。
模拟退火算法的模型
模拟退火算法可以分解为解空间、目标函数和初始解三部分。
模拟退火的基本思想:
(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L
(2) 对k=1，……，L做第(3)至第6步：
(3) 产生新解S′
(4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数
(5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)接受S′作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解，结束程序。
终止条件通常取为连续若干个新解都没有被接受时终止算法。
(7) T逐渐减少，且T->0，然后转第2步。
算法对应动态演示图：
模拟退火算法新解的产生和接受可分为如下四个步骤：
第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropo1is准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。
第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。
模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。
模拟退火算法的简单应用
作为模拟退火算法应用，讨论货郎担问题(Travelling Salesman Problem，简记为TSP)：设有n个城市，用数码1,…,n代表。城市i和城市j之间的距离为d(i，j) i, j=1,…,n．TSP问题是要找遍访每个域市恰好一次的一条回路，且其路径总长度为最短.。
求解TSP的模拟退火算法模型可描述如下：
解空间解空间S是遍访每个城市恰好一次的所有回路，是{1，……，n}的所有循环排列的集合，S中的成员记为(w1,w2 ,……，wn)，并记wn+1= w1。初始解可选为(1，……，n)
目标函数此时的目标函数即为访问所有城市的路径总长度或称为代价函数：
我们要求此代价函数的最小值。
新解的产生 随机产生1和n之间的两相异数k和m，若k，则将
(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn)
变为：
(w1, w2 ,…，wm , wm-1 ,…，wk+1 , wk ,…，wn).
如果是k>m，则将
(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn)
变为：
(wm, wm-1 ,…，w1 , wm+1 ,…，wk-1 ,wn , wn-1 ,…，wk).
上述变换方法可简单说成是“逆转中间或者逆转两端”。
也可以采用其他的变换方法，有些变换有独特的优越性，有时也将它们交替使用，得到一种更好方法。
代价函数差 设将(w1, w2 ,……，wn)变换为(u1, u2 ,……，un), 则代价函数差为：
根据上述分析，可写出用模拟退火算法求解TSP问题的伪程序：
Procedure TSPSA:
begin
init-of-T; { T为初始温度}
S={1，……，n}; {S为初始值}
termination=false;
while termination=false
begin
for i=1 to L do
begin
generate(S′form S); { 从当前回路S产生新回路S′}
Δt:=f(S′))-f(S);{f(S)为路径总长}
IF(Δt<0) OR (EXP(-Δt/T)>Random-of-[0,1])
S=S′;
IF the-halt-condition-is-TRUE THEN
termination=true;
End;
T_lower;
End;
End
模拟退火算法的应用很广泛，可以较高的效率求解最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)、调度问题(Scheduling Problem)等等。
模拟退火算法的参数控制问题
模拟退火算法的应用很广泛，可以求解NP完全问题，但其参数难以控制，其主要问题有以下三点：
(1) 温度T的初始值设置问题。
温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。实际应用过程中，初始温度一般需要依据实验结果进行若干次调整。
(2) 退火速度问题。
模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索(退火)是相当必要的，但这需要计算时间。实际应用中，要针对具体问题的性质和特征设置合理的退火平衡条件。
(3) 温度管理问题。
温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中，由于必须考虑计算复杂度的切实可行性等问题，常采用如下所示的降温方式：
T(t+1)＝k×T(t)
式中k为正的略小于1.00的常数，t为降温的次数。
模拟退火算法的主要思想
就函数最小值问题来说，模拟退火的主要思想是：在搜索区间（二维平面中）随机游走（即随机选择点），再以Metropolis抽样准则，使随机游走逐渐收敛于局部最优解。而温度即是Metropolis算法中的一个重要控制参数，可以认为这个参数的大小控制了随时过程向局部或全局最优解移动的快慢。
冷却参数表、领域结构和新解产生器、接受准则和随机数产生器（即Metropolis算法）一起构成算法的三大支柱。
l 重点抽样与Metroplis算法：
Metropolis是一种有效的重点抽样法，其算法为：系统从能量一个状态变化到另一个状态时，相应的能量从E1变化到E2，概率为p = exp[ - (E2- E1)/kT ]。如果E2 < E1，系统接收此状态，否则，以一个随机的概率接收此或丢弃此状态。这种经常一定次数的迭代，系统会逐渐趋于一引稳定的分布状态。
重点抽样时，新状态下如果向下则接受（局部最优），若向上（全局搜索），以一定机率接受。模拟退火方法从某个初始解出发,经过大量解的变换后,可以求得给定控制参数值时组合优化问题的相对最优解。然后减小控制参数T的值，重复执行Metropolis算法，就可以在控制参数T趋于零时，最终求得组合优化问题的整体最优解。控制参数的值必须缓慢衰减。
其中温度是一个Metropolis的重要控制参数，模拟退火可视为递减控制参数什时Metroplis算法的迭代。开始T值大，可能接受较差的恶化解，随着T的减小，只能接受较好的恶化解，最后在T趋于0时，就不再接受任何恶化解了。
在无限高温时，系统立即均匀分布，接受所有提出的变换。T的衰减越小，T到达终点的时间越长；但可使马可夫链越小，到达准平衡分布的时间越短，
l 参数的选择：
我们称调整模拟退火法的一系列重要参数为冷却进度表。它控制参数T的初值及其衰减函数，对应的MARKOV链长度和停止条件，非常重要。
一个冷却进度表应当规定下述参数：
1． 控制参数t的初值t0;
2． 控制参数t的衰减函数；
3．马尔可夫链的长度Lk。（即每一次随机游走过程，要迭代多少次，才能趋于一个准平衡分布，即一个局部收敛解位置）
4． 结束条件的选择
有效的冷却进度表判据：
一．算法的收敛：主要取决于衰减函数和马可夫链的长度及停止准则的选择
二．算法的实验性能：最终解的质量和CPU的时间
参数的选择：
一）控制参数初值T0的选取
一般要求初始值t0的值要充分大，即一开始即处于高温状态，且Metropolis的接收率约为1。
二）衰减函数的选取
衰减函数用于控制温度的退火速度，一个常用的函数为：T(n + 1) = K*T(n)，其中K是一个非常接近于1的常数。
三）马可夫链长度L的选取
原则是：在衰减参数T的衰减函数已选定的前提下，L应选得在控制参数的每一取值上都能恢复准平衡。
四）终止条件
有很多种终止条件的选择，各种不同的条件对算法的性能和解的质量有很大影响，本文只介绍一个常用的终止条件。即上一个最优解与最新的一个最优解的之差小于某个容差，即可停止此次马尔可夫链的迭代。
使用模拟退火法求函数f(x,y) = 5sin(xy) + x2 + y2的最小值
解：根据题意，我们设计冷却表进度表为：
即初始温度为100
衰减参数为0.95
马可夫链长度为10000
Metropolis的步长为0.02
结束条件为根据上一个最优解与最新的一个最优解的之差小于某个容差。
模拟退火算法在自动排课系统中的应用
using System;
using System.Collections.Generic;
using System.Text;
using System.Xml;
using System.IO;
using System.Runtime.Serialization;
using System.Runtime.Serialization.Formatters.Binary;
using System.Windows;
using System.Threading;
using System.Windows.Forms;
namespace WinPaike
{
public class CoursePriority
{
public static int[] Priority = new int[ClassUnit.CourseCount] { 5, 5, 4, 4, 2, 1 };
}
public class PaiKe
{
int cc = 0;//当前要变异的班级序号
const int MaxCourseWay = 10;
List<string> ResolveList = new List<string>(); //保存最后结果
int[] next;
#region 线程操作所用
public TextBox TxtBox;
public Form Frm;
public bool IsChongTu = true;
public ToolStripProgressBar Bar;
delegate void SetNumber(int value, ToolStripProgressBar bar);
delegate void SetText(string text, Control txtbox);
SetText addtext ;//= new SetText(AddControlText);
SetText settext ;//= new SetText(SetControlText);
SetNumber setnumber;// = new SetNumber(SetControlNumber);
SetNumber setmaxvalue ;//= new SetNumber(SetMaxValue);
void AddControlText(string text, Control txtbox)
{
txtbox.Text += text;
}
void SetControlText(string text, Control txtbox)
{
txtbox.Text = text;
}
void SetControlNumber(int value, ToolStripProgressBar bar)
{
bar.Value = value;
}
void SetMaxValue(int value, ToolStripProgressBar bar)
{
bar.Maximum = value;
}
#endregion
Random rnd = new Random();
ClassUnit tClassUnit;
List<ClassUnit> ClassList = new List<ClassUnit>();
List<Course> CourseList = new List<Course>();
//－－－－－－－－－－－－－－－惩罚值列表－－－－－－－－－－－－－－
//冲突
const int ChongTu = 200;
//一天有2节相同的课
const int OneDayForTowCourse = 150;
//连续2天有2节相同课
const int TowDayForTowCourse = 40;
//课程优先级单位值
const int CoursePr = 15;
public PaiKe()
{
addtext = new SetText(AddControlText);
settext = new SetText(SetControlText);
setnumber = new SetNumber(SetControlNumber);
setmaxvalue = new SetNumber(SetMaxValue);
}
//初始化
void Init()
{
ResolveList.Clear();//结果清零
CourseList = CommonClass.GetCourseListFromDB("Course.dat");
Course tcourse = new Course("自习", -1, "NO", 1);
tcourse.ID = -1;
CourseList.Add(tcourse);
List<CourseInClass> courseinclasslist = new List<CourseInClass>();
courseinclasslist = CommonClass.GetCourseInClassListFromDB("courseinclass.dat");
ClassList = CommonClass.GetClassListFromDB("Class.dat");
int[] txulie = new int[ClassList.Count];
foreach (CourseInClass courseinclass in courseinclasslist)
{
for (int i = 0; i < courseinclass.Count; i++)
{
int classid = courseinclass.ClassID - 1;
if (txulie[classid] < ClassUnit.CourseCount * ClassUnit.WeekDay)
{
ClassList[classid].XuLie[txulie[classid]++] = courseinclass.CourseID;
}
}
}
//将未添满的课设置为自习课
for (int i = 0; i < ClassList.Count; i++)
{
for (int j = txulie[i]; j < ClassUnit.CourseCount * ClassUnit.WeekDay; j++)
{
ClassList[i].XuLie[j] = -1;
}
}
//ClassList.RemoveAt(ClassList.Count-1);
next = new int[ClassUnit.WeekDay * ClassUnit.CourseCount];
Assign(next,ClassList[0].XuLie);
}
//得到单个班级的课程表
string GetCourseTable(ClassUnit cu)
{
string re = " ";
for (int i = 0; i < ClassUnit.CourseCount; i++)
{
re += string.Format("第{0}节 ", i + 1);
}
re += string.Format("\r\n－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－\r\n");
for (int i = 0; i < ClassUnit.WeekDay * ClassUnit.CourseCount; i++)
{
if ((i) % (ClassUnit.CourseCount) == 0)
{
re += string.Format("周{0}: ", i / ClassUnit.CourseCount + 1);
}
int tcourseid = cu.XuLie[i];
if (tcourseid > 0)
{
re += string.Format(CourseList[(cu.XuLie[i] - 1)].Name.ToString() + " ");
}
else
{
re += string.Format("自习" + " ");
}
if ((i + 1) % ClassUnit.CourseCount == 0)
{
re += string.Format("\r\n");
}
}
return re;
}
//得到所有班级的课程表
string GetCourseTableList()
{
string re = "";
for (int i = 0; i < ClassList.Count; i++)
{
re += string.Format("班级{0}:\r\n", ClassList[i].Name);
re += GetCourseTable(ClassList[i]);
}
return re;
}
//复制数组
void Assign(int[] c, int[] n)
{
for (int i = 0; i < c.Length; i++)
c[i] = n[i];
}
//产生下一组解（班级课表变异）
void CreateNext(int[] c, int[] n)
{
for (int i = 0; i < c.Length; i++)
n[i] = c[i];
int i1 = (int)(rnd.Next(ClassUnit.WeekDay * ClassUnit.CourseCount));
int i2 = (int)(rnd.Next(ClassUnit.WeekDay * ClassUnit.CourseCount));
int aux = n[i1];
n[i1] = n[i2];
n[i2] = aux;
}
//计算单个班级的除冲突外的惩罚代价
int ComputeClassValue(int[] cu)
{
int Sum = 0;
for (int i = 0; i < ClassUnit.WeekDay; i++)
{
int row = i * ClassUnit.CourseCount;
for (int j = row; j < row + ClassUnit.CourseCount; j++)
{
//一天2节
for (int k = j + 1; k < row + ClassUnit.CourseCount; k++)
{
if (cu[j] == cu[k])
{
Sum += OneDayForTowCourse;
}
}
//2天2节
for (int k = (j / ClassUnit.CourseCount + 1) * ClassUnit.CourseCount; k < (j / ClassUnit.CourseCount + 1) * ClassUnit.CourseCount + ClassUnit.CourseCount && k < ClassUnit.WeekDay * ClassUnit.CourseCount; k++)
{
if (cu[j] == cu[k] && j != k)
{
Sum += TowDayForTowCourse;
}
}
//4天2节
for (int k = (j / ClassUnit.CourseCount + 3) * ClassUnit.CourseCount; k < (j / ClassUnit.CourseCount + 3) * ClassUnit.CourseCount + ClassUnit.CourseCount && k < ClassUnit.WeekDay * ClassUnit.CourseCount; k++)
{
if (cu[j] == cu[k] && k != j)
{
Sum += TowDayForTowCourse;
}
}
int t1 = CoursePriority.Priority[j % ClassUnit.CourseCount];
if (cu[j] > 0)
{
int t2 = CourseList[cu[j] - 1].Priority;
if (t1 < t2)
{
Sum += CoursePr * (t2 - t1);// * (Math.Max((j % 4 + 1 - cu[j]), 0));
}
}
}
}
return Sum;
}
//计算总的排课方案代价
int ComputeValue()
{
int value;
value = ComputeClassValue(next);
for (int i = 0; i < ClassList.Count; i++)
{
if (i == cc) { continue; }
value += ComputeClassValue(ClassList[i].XuLie); ;
}
value += ComputeConflictValue();
return value;
}
//输出可行的排课方案
string GetResolve()
{
string txt = "";
txt += GetCourseTableList();
return txt;
}
//计算冲突所产生的惩罚代价
int ComputeConflictValue()
{
int Sum = 0;
int t1, t2;
int tcourseid1, tcourseid2;
for (int i = 0; i < ClassUnit.WeekDay * ClassUnit.CourseCount; i++)
{
for (int j = 0; j < ClassList.Count; j++)
{
if (j == cc)
{
tcourseid1 = next[i];
}
else
{
tcourseid1 = ClassList[j].XuLie[i];
}
if (tcourseid1==-1)
{
continue;
}
t1 = CourseList[tcourseid1-1].TeacherID;
for (int k = j + 1; k < ClassList.Count; k++)
{
if (k == cc)
{
tcourseid2 = next[i];
}
else
{
tcourseid2 = ClassList[k].XuLie[i];
}
if (tcourseid2 == -1)
{
continue;
}
t2 = CourseList[tcourseid2 - 1].TeacherID;
//int t1 = CourseList[ClassList[j].XuLie[i] - 1].TeacherID;
//int t2 = CourseList[ClassList[k].XuLie[i] - 1].TeacherID;
//if (j == cc) { t1 = CourseList[next[i] - 1].TeacherID; }
//if (k == cc) { t2 = CourseList[next[i] - 1].TeacherID; }
if (t1 == t2)
{
Sum += ChongTu;
}
}
}
}
if (Sum == 0)
{
IsChongTu = false;
}
else
{
IsChongTu = true;
}
return Sum;
}
//计算结束,输出排课方案
void Finished()
{
SetText add = new SetText(AddControlText);
if (ResolveList.Count == 0)
{
MessageBox.Show("对不起，没有找到可行解，请重新尝试！");
}
else
{
MessageBox.Show("成功完成！");
}
int ct = 0;
foreach (string ss in ResolveList)
{
ct++;
TxtBox.Invoke(add, new object[2] { string.Format("方案{0}\r\n", ct), TxtBox });
TxtBox.Invoke(add, new object[2] { ss, TxtBox });
}
}
//主函数，开始排课
public void Run()
{
int iteration = -1;//迭代次数
double proba;//产生的0～1间的随机数
double alpha = 0.999;//降温速率
double temperature = 4000.0 + (100 * ClassList.Count) * (100 * ClassList.Count);//温度
double temperature_copy = temperature;//温度副本
double epsilon = 0.002;//终止温度
double delta;//
int i;
double Value = 0;//当前代价
#region 待优化(思想：将每个班的惩罚代价进行比较，然后尽量选取代价大的进行变异)
/*
double[] XuanZhongGaiLv = new double[ClassList.Count];//每个班级的变异概率，总和＝1
for (i = 0; i < XuanZhongGaiLv.Length; i++)
{
XuanZhongGaiLv[i] = 1 / XuanZhongGaiLv.Length;
}
*/
#endregion
//初始化班级，课程等信息
Init();
Value = ComputeValue();
Value += (Value + 20000);
int battc = 0;
TxtBox.Invoke(setmaxvalue, new object[2] {(int)(temperature), Bar });
//while the temperature didnt reach epsilon
while (temperature > epsilon)
{
battc++;
if (battc > 100)
{
//Thread.Sleep(100);
TxtBox.Invoke(setnumber, new object[2] { (int)(temperature_copy-temperature), Bar });
battc = 0;
}
iteration++;
cc = cc % ClassList.Count;
//cc = rnd.Next(ClassList.Count);
CreateNext(ClassList[cc].XuLie, next);
delta = ComputeValue() - Value;
if (delta == 0 && ResolveList.Count <MaxCourseWay && !IsChongTu)
{
ResolveList.Add(GetResolve());
}
#region 如果新生成的解更优，则替换以前的，否则以接受概率接受新解
//如果新生成的解更优，则替换以前的
if (delta < 0)
{
Assign(ClassList[cc].XuLie, next);
Value = delta + Value;
if (!IsChongTu)
{
ResolveList.Clear();
//输出代价
TxtBox.Invoke(settext, new object[2] { "自动排课" + "------当前最优值:(" + Value.ToString() + ")", Frm });
ResolveList.Add(GetResolve());
}
}
else
{
proba = rnd.Next();
if (proba < Math.Exp(-delta / temperature))
{
Assign(ClassList[cc].XuLie, next);
Value = delta + Value;
}
}
#endregion
//降温
temperature *= alpha;
cc++;
}
Finished();
}
}
}
|
