基于模拟退火算法的PFSP研究（最优化方法实验附github源码）
基于模拟退火算法的PFSP研究
摘 要: 置换流水车间调度问题（PFSP）是典型的NP 问题,其研究的是假设有若干工件和机器，不同工件在不同机器上加工时间不同，每个机器加工工件的顺序相同，每个工件加工的机器顺序是固定的，求解一个加工工件的顺序使得总完工时间最小。近年来随着智能优化算法的出现和发展，解决车间生产调度问题也有了新的思路。本文采用的是模拟退火算法，模拟退火算法（Simulated Annealing, SA）是一种以蒙特卡罗迭代策略为基础寻找全局最优解的技术。它基于固体物质的退火过程,在温度由高变低的过程中,通过合理设置参数,并根据接受准则在相应的解空间中寻找目标函数的最优解。由于这种算法在组合应用问题中应用很广泛,因此适合用来解决PFSP。该算法将总工时降低到一个比较趋近全局最优解的水平，易于跳出局部最优解，效果良好。
关键词: 模拟退火算法；置换流水车间调度问题；PFSP
1 引言
传统的FSP问题是针对n个工件在m台机器上流水加工的过程进行研究。在有限的资源约束下，n个工件在m台机器上进行加工，单个工件需多项操作才可完成，且每个工件在每台机器上加工顺序相同，同时假定每个工件在每台机器上只加工一次，每台机器在某一时刻一次只能够加工一个工件。已知单个工件在每台机器上的加工时间和准备时间，目的是确定使某项生产指标最优的调度方案。如果每台机器上加工的各工件顺序相同，则称为置换流水车间调度问题(permutation flow-shop scheduling problem,PFSP) 。虽然置换流水车间调度问题形式上比较简单，但3台机器以上的置换流水车间调度问题已经被证明属于NP-Hard问题[1]。到目前为止还没有一个具有多项式计算复杂性的全局优化算法。为了能够在合理的时间内快速得到问题的最优解，开发研究高速有效的优化技术显得尤为重要[2]。
在PFSP 问题的研究中，大部分研究文献都将最小化最大完工时间作为性能衡量指标。假设\(p_{i,j}\)表示工件i在机器j上的加工时间（假设各个工件的加工准备时间为零或者已经包含在加工时间内）；\(\pi=(j_1,j_2,\dots , j_n)\)表示工件集合的一个排序；\(C(j_i,k)\)表示工件\(j_{i}\)在机器上的加工完成时间。数学描述如下：
首先，Makespan可被定义为：
\(C(j_i,1)=C(j_{i-1},1)+p_{j_i,k},k=2,\dots,n;\quad\quad(4)\)
\(C(j_i,k)=C(j_1,k-1)+p_{j_i,k},k=2,\dots,m;\quad\quad(5)\)
\(C_{j_i,k}=max\{ C(j_{i-1},k),C(j_i,k-1) \}+p_{j_i,k},i=2,\dots,n;k=2,\dots,m;\quad\quad(6)\)
本文采用的是模拟退火算法。其中的步进规则是随机交换某两个零件的加工顺序，在计算目标函数时模拟时间的进行，先确定每台机器将要加工的工件，然后再处理机器和工件的剩余加工时间，得到加工结束时间。最后对不同用例适当调整温度以及其他参数可以分别使其在较短时间内求出最优解。
本文后续部分组织如下。第2节详细陈述使用的方法，第3节报告实验结果，第4节对讨论本文的方法并总结全文。
2 算法设计
2.1 模拟退火算法思路简介
模拟退火算法（Simulated Annealing, SA）是一种以蒙特卡罗迭代策略为基础寻找全局最优解的技术。它基于固体物质的退火过程,在温度由高变低的过程中,通过合理设置冷却进度表中的参数,并根据接受准则在相应的解空间中以一定的步进规则寻找目标函数的最优解。在这个问题中的步进规则是随机交换某两个零件的加工顺序，在计算目标函数时模拟时间的进行，先确定每台机器将要加工的工件，然后再处理机器和工件的剩余加工时间，得到加工结束时间。在温度下降的过程中，若超过50次未接受新解或者超过100次未发生目标函数的变化则得到最优解。
2.2 算法详述
2.2.1 模拟退火算法
考虑一个优化函数为\(f:x\rightarrow R^+\)，其中\(x\in S\)，它表示优化问题的一个可行解，\(R^+=\{ y|y\in R,y>0 \}\)，S表示函数的定义域。\(N(x)\subseteq S\)表示x的一个邻域集合。
首先给定一个初始温度\(T_0\)和该优化问题的一个初始解x(0)，并由x(0)生成下一个解\(x'\in N(x(0))\)，是否接受x'作为一个新解依赖于下面概率：
\(P(x(0)\rightarrow x')=\left\{\begin{matrix} 1 & 若f(x')<f(x(0))\\ e^{-\frac{f(x')-f(x(0))}{T_0} } & 其他 \end{matrix}\right.\quad \quad (7)\)
即，如果生成的解x'的函数值比前一个解的函数值更小，则接受x(1)=x'作为一个新解。否则以概率$e^{-\frac{f(x')-f(x(0))}{T_0}}$接受x'作为一个新解。
对于某一个温度$T_i$和该优化问题的一个解x(k)，可以生成x'。接受x'作为下一个新解x(k+1)的概率为：
\(P(x(k)\rightarrow x')=\left\{\begin{matrix} 1 & 若f(x')<f(x(k))\\ e^{-\frac{f(x')-f(x(k))}{T_0} } & 其他 \end{matrix}\right.\quad \quad (8)\)
在温度$T_i$下，经过很多次的转移之后，降低温度$T_i$，得到$T_{i+1}
\(P_i(T_i)=\frac{e^{-\frac{f(x_i)}{T} }}{ {\textstyle \sum_{j\in S}^{}}e^{-\frac{f(x_i)}{T_i} } } \quad \quad (9)\)
当温度T降为0时，$x_i$的分布为：
\(P_o^*=\left\{\begin{matrix} \frac{1}{S_{min}} & 若S_i\in S_{min}\\ 0 & 其他 \end{matrix}\right.\quad \quad (10)\)
并且
\(\sum_{x_i\in S_{min}}P_i^*=1 \quad \quad (11)\)
2.2.2 Metropolis准则
Metropolis准则是最常用的接受准则。在一定的温度下，状态由当前状态i转换为状态j，他们对应的目标函数为f(i)和f(j)。若\(f(i)\ge f(j)\)，则接受新状态。若f(i)<f(j)，则状态转换以如下概率被接受：
\(q=e^{\frac{f(i)-f(j)}{KT} }\quad \quad (12)\)
同时，在区间[0,1)间随机生成一个数\(\lambda\)。Metropolis准则表明，若\(q>\lambda\)则接受新状态j，否则舍弃它。
2.2.3 目标函数与步进规则
在算法过程中，目标函数标志当前状态效果的优劣，在本题中目标函数就是加工结束时间\(C(j_n,m)\)，其中\(j_n\)为序列中第n个工件，m为第m台机器，即序列中第n个工件在第m台机器中加工结束时间。
模拟退火算法中重要的一环就是当前状态如何转化到下一状态。这个过程实现的规则又称步进规则。在本题，步进规则遵循随机生成原则：随机选取序列中的两个工件交换顺序。
2.2.4 初始状态
其他实验表明[5]，对于模拟退火算法来说，一个好的初始状态有很重要的意义，它可以大幅提升算法的优化结果并减少优化时间。很多模拟退火的运用实例中都用其它优化算法得到一个较优的初始解，再用模拟退火算法进一步优化，实现模拟退火和其他算法的结合使用。
本题对于初始解的生成采用了随机策略，即随机生成一个工件序列。
2.2.5 算法伪代码
本文所使用的模拟退火算法不同于平常的模拟退火算法，是以超过Max_management_Maintain次未接受新解或者超过Max_f_Unchanged次未发生目标函数的变化作为算法结束条件的，基于本题的模拟退火算法伪代码如下：
1. 初始化参数T、冷却因子theta、连续未接受新解计数器count、目标函数连续重复次数count1、当前解目标函数值t、未接受新解次数上限Max_management_Maintain和目标函数值未变化上限Max_f_Unchanged
2. 生成初始解$S_{init}$，令当前解S=$S_{init}$，其目标函数值为f(S)
3. for t←1 to ∞ do
4. T←T*theta
5. S步进
6. 重置工件和机器状态
7. new_t←f(S)
8. if new_t < t then
9. t←new_t
10. 重置两计数器
11. else
12. 生成0到1之间的随机浮点数pp
13. if pp<$e^{\frac{t-new\_t}{T} }$ then
14. 接受新解且更新两计数器
15. else
16. S回退到上一个状态
17. 更新两计数器
18. If count>Max_management_Maintain || count1>Max_f_Unchanged then break
2.2.6 算法核心部分的复杂度分析
若设置了初始温度T、最小温度\(T_{min}\)和内循环次数b，n为工件数量，m为机器数量，theta为冷却因子，每次计算目标函数的时间复杂度为O(mntmax)，因为计算目标函数的复杂度是以时间为循环因子的，而时间可以大致认为是小于tmaxmn的（此处的\(t_{max}\)为某工件在某机器加工的最大加工时间）。所以总的计算复杂度为：
\(O(t_{max}*m*n*log_{theta}\frac{T_{min}}{T} *b)\)
但是由于本文采用的方法不同于普通的模拟退火算法，所以计算时间复杂度十分复杂。但仍可以给出粗略的估计，分析如下：
先将参数Max_management_Maintain记为M，当出现M次未接受新解的概率为0.99时，\(e^{\frac{t-new\_t}{T_{min}} }\)约为\(1-0.99^{\frac{1}{M}}\)时，此时\(e^{\frac{t-new\_t}{T_{min}} }\)约为ln(\(1-0.99^{\frac{1}{M}}\))，假设最后的每次步进是最大幅度的，则\(new\_t-t>2*t_{max}\)，这里的\(t_{min}\)是某工件在某机器加工的最小加工时间，所以\(T_{min}>2*t_{min}/(-ln(1-0.99^{\frac{1}{M} }))\).接下来我们对该式子进行一定的放缩：
3 实验
3.1 实验设置
模拟退火主程序是在java的jdk11.0.12的环境下运行的，运行时间随参数变化可以短到毫秒，也可以长至半分钟。主要参数有温度T、降温因子theta、未接受新解次数上限Max_management_Maintain和目标函数值未变化上限Max_f_Unchanged。
3.2 实验结果
见表1。
表1 实验结果
| instance | 调度方案 | 加工工时 |
|---|---|---|
| 0 | [7, 4, 2, 10, 1, 8, 6, 0, 3, 5, 9] | 7038 |
| 1 | [4, 1, 2, 0, 5, 3] | 6923 |
| 2 | [6, 2, 3, 10, 7, 8, 1, 0, 4, 5, 9] | 6309 |
| 3 | [10, 5, 11, 13, 12, 9, 4, 8, 1, 3, 2, 7, 6, 0] | 8244 |
| 4 | [3, 11, 12, 5, 8, 6, 0, 13, 9, 10, 14, 1, 2, 15, 7, 4] | 9373 |
| 5 | [4, 1, 3, 0, 2, 7, 5, 9, 8, 6] | 7720 |
| 6 | [1, 13, 8, 3, 9, 15, 11, 18, 7, 10, 19, 12, 4, 2, 14, 0, 17, 16, 6, 5] | 1431 |
| 7 | [5, 14, 8, 4, 12, 0, 19, 13, 16, 7, 18, 2, 15, 11, 17, 9, 10, 3, 1, 6] | 1973 |
| 8 | [13, 1, 0, 6, 5, 18, 2, 10, 3, 9, 7, 16, 15, 4, 12, 8, 19, 14, 11, 17] | 1111 |
| 9 | [12, 17, 16, 1, 14, 11, 4, 7, 8, 18, 6, 10, 2, 3, 13, 9, 15, 5, 0, 19] | 1935 |
| 10 | [12, 13, 39, 49, 34, 8, 35, 9, 1, 2, 41, 26, 47, 43, 24, 45, 37, 18, 44, 36, 33, 3, 6, 4, 17, 22, 28, 21, 20, 38, 10, 29, 40, 16, 7, 19, 25, 31, 11, 14, 42, 46, 5, 32, 48, 0, 27, 23, 15, 30] | 3280 |
这些调度方案和对应的加工工时都是在运行了多次程序多次出现的最优解。所以在一般情况下运行程序得到的结果可能和最优解偏差0~0.2%。
3.2.1 用例角度的结果分析
对于instance0，如果我们只有条件超过50次未接受新解的话，它会在两个目标函数相等的地方来回跳变以至于难以跳出循环，所以我们添加了新的条件，以目标函数值超过100次未改变为条件来跳出循环。
对instance6和instance7可以将theta改为0.99999。对instance8和instance9可以将theta改为0.999，这样运行得更快而且更容易收敛到最优解。对其他用例theta采用默认值0.9999可以得到较好的效果。
对instance6、instance7和instance9可以将Max_management_Maintain和Max_f_Unchanged全部设置为250，因为这些用例中的有很多加工时间较少的工件，这样可以在温度较低时进行微调，找到更加趋近于最优解的解。
图1 instance6目标函数变化图
3.2.2 参数角度的结果分析
在实验中初始温度一般为10000，theta一般为0.9999。若将温度T或降温因子theta从合适值开始降低，运行时间会有所降低，所得到的结果方差也会增大，对某些用例将降温因子theta从合适值开始提高不一定能提高解的质量，比如将instance1中的theta改为0.9999反而难以找到最优解。
一般将Max_management_Maintain和Max_f_Unchanged都设置为100会得到较好的结果。将Max_management_Maintain和Max_f_Unchanged升高，这样可以在温度较低时进行微调，找到更加趋近于最优解的解。
初始温度高则在前几次步进中目标函数变化会比较快，见图2，增大温度适用于解空间较大的情况。
图2 不同温度下的目标函数变化
4 总结
本文采用了不同于平常的模拟退火算法的方式对置换流水车间调度问题的一些实例进行了寻找最优解。步进规则是随机交换某两个零件的加工顺序，在计算目标函数时模拟时间的进行，先确定每台机器将要加工的工件，然后再处理机器和工件的剩余加工时间，得到加工结束时间。本文所用算法是以超过50次未接受新解或者超过100次未发生目标函数的变化作为算法结束条件的。
笔者对每个用例进行了多次实验和调参，我的感悟是步进规则要尽量温和，不然可能难以收敛到最优解。参数温度和冷却因子并不是越大越好，关键是要找出最合适的参数。在实验中可以作出目标函数变化图来分析和调参。
然而本文进行的参数改动仍然较少，可能并未对每个用例找到一个兼顾时间和得到最优解的概率的最合适的参数。本文没有对调度结果图象化，可以考虑用python作出甘特图。
详细代码见https://github.com/Tracker1701/SA_PFSP
参考文献:
[1] Garey M R,Johnson D S.Computers and intractability: a guide to the theory of NP-completeness[M].San Francisco:Freeman,1979．
[2] 盛晓华,叶春明.蝙蝠算法在PFSP调度问题中的应用研究[J].工业工程,2013,16(01):119-124.
[3] 于承敏,郑丽萍,张民.粒子群算法解决置换流水车间调度问题方法综述[J].机械设计与制造,2012(08):84-86.DOI:10.19356/j.cnki.1001-3997.2012.08.031.
[4] 庞峰. 模拟退火算法的原理及算法在优化问题上的应用[D].吉林大学,2006.
[5] 陈华根,李丽华,许惠平. 模拟退火定位算法研究[J]. 同济大学学报（自然科学版）,2005,33(9):1240-1243. DOI:10.3321/j.issn:0253-374X.2005.09.019.
