简介
模拟退火算法求解车间调度问题,评估函数采用时间模拟的形式。
接下来给出各部分的介绍
车间调度问题
本文研究的车间调度问题:有N个工件(W1,W2,…,WN),每个工件都有有M道工序(M1,M2,…,MM)得在相应的机器(M个)上加工,每个工件都得从第一道工序顺序加工至最后一道工序,每个工件的每道工序都有一定的时间,每个工件在同一时间只能被一个机器加工,一个机器同一时间只能加工一个工件,工件的每一道工序不可被打断,求得一个调度方案:每个机器上,工件得加工顺序,使得整个加工流程所用时间最少。
我们的问题其实是一道优化问题,这里给出其数学描述:
记矩阵
T
N
×
M
T_{N×M}
TN×M为工件工序得加工时间,设
S
i
j
S_{ij}
Sij为工件i第j道工序的开始时间。则整个优化问题可以表述为:
m
i
n
(
m
a
x
(
S
i
M
+
T
i
M
)
)
,
i
=
1
,
2
,
…
,
N
.
(
1
)
min(max(S_{iM} + T_{iM})), i = 1,2,…,N.\qquad(1)
min(max(SiM+TiM)),i=1,2,…,N.(1)
s
.
t
.
S
i
1
≥
0
,
i
=
1
,
2
,
…
,
N
.
(
2
)
s.t. \quad S_{i1} ≥ 0 , i = 1,2,…,N.\qquad(2)
s.t.Si1≥0,i=1,2,…,N.(2)
S
i
j
≥
S
i
(
j
−
1
)
+
T
i
(
j
−
1
)
,
i
=
1
,
2
,
…
,
N
,
j
=
2
,
…
,
M
.
(
3
)
S_{ij} ≥ S_{i(j-1)} + T_{i(j-1)}, i = 1,2,…,N, j = 2,…,M.\qquad(3)
Sij≥Si(j−1)+Ti(j−1),i=1,2,…,N,j=2,…,M.(3)
S
j
k
+
T
j
k
≤
S
i
k
o
r
S
j
k
≥
S
i
k
+
T
i
k
,
i
=
1
,
2
,
…
,
N
,
j
=
1
,
…
,
i
−
1
,
i
+
1
,
…
,
N
,
k
=
1
,
2
,
…
,
M
.
(
4
)
S_{jk} + T_{jk} ≤ S_{ik} or S_{jk} ≥ S_{ik} + T_{ik}, i = 1,2,…,N, j = 1,…,i-1,i+1,…,N, k = 1,2,…,M.\qquad(4)
Sjk+Tjk≤SikorSjk≥Sik+Tik,i=1,2,…,N,j=1,…,i−1,i+1,…,N,k=1,2,…,M.(4) (1)式是这个问题的优化目标,即所有工件所有工序都加工完成所有时间得最小值,max(SiM + TiM) , i = 1,2,…,N即为加工结束时间。第一个约束条件((2)式)表明:所有工件的第一道工序可以从时间0开始。第二个约束条件((3)式)表明:每一个工件的某道工序开始都以前一道工序结束为前提。第三个约束条件((4)式)表明:在某个工件的某道工序时间内,其他工件都不得进行该道工序。
作为一个线性约束的最优化问题,在模拟退火的过程中,其实可以用线性规划的算法作为评估函数,然而笔者水平有限,没学会线性规划算法,转而求其次实现一个模拟时间的算法。立个flag:下一篇博客介绍线性规划算法。
模拟退火算法
给出模拟退火算法的官方定义:模拟退火算法来源于固体退火原理,是一种基于概率的算法,将固体加温至充分高,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状,内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。
通俗的来说,就是我们把目标函数的值转化为上述过程中的内能(越是好的解,内能越小),然后我们设置一个循环,且循环过程中有一个变量作为温度(随着循环次数衰减),每次循环都搜索一个新解,如果这个新解的内能比旧解的内能小(退火过程的内能下降),那么我们将这个解作为目前最优解,否则我们以与温度有关的一个指数概率接受这个新解。
该算法的伪代码如下图:
图中有两点值得注意:
- 温度的衰减:可以设置一个慢慢减少的数组,也可以采用其他衰减方法例如: T n = β T n − 1 T_n = βT_{n-1} Tn=βTn−1。
- 算法停止条件:可以是多次未找到新解,也可以是接受的新解能量多次变化很小。
这里展示这一部分对应的代码:
int t = 0x3f3f3f3f;//表示无穷大,其实不是无穷大,不懂得同学可以去搜一下
int flag = 0;
double T = T0;//初始温度
for (int i = 0; ; i++) {
std::cout << i << '\t' << t << std :: endl;
T = T * theta; //当前温度
S.randomExchange(); //产生新解
S.resetForNewRound();//重置工件机器状态
int new_t = S.evaluate();//评估新解
if (new_t < t) {
t = new_t;
flag = 0;
continue;
}//接受新解
else {
double p = std::exp((t - new_t) / T);
if (((std::rand() % 10000) / 10000.0) <= p) {
t = new_t;
flag = 0;
continue;
}//以概率接受新解
else {
S.rollBack();//回退至上一状态
flag++;
if (flag > 50) //如果长时间未接受新解,结束搜索
break;
}//拒绝新解
}
}
其中的S是我们自定义的一个类型的变量,接下来会介绍,并且所有代码在附录的链接里。1
评估函数部分
前面提到,由于不会线性规划,于是采用模拟时间进行的方法来实现评估函数。所谓时间模拟其实就是一个循环,循环的轮数代表时间,每一次循环代表当前时间颗粒(模拟的一个时间单元)的开始到结束,可以理解为图中形式:
这也就是说,我们在每一次循环中都得根据上一个时间颗粒整个系统的末尾状态来得到这个时间颗粒系统的末尾状态。为了更好的描述整个系统,我们实现了一个类simulate来描述我们关注的加工流水线,声明了另外两个类Machine,Workpiece来分别描述机器和工件,它们的声明如下:
class Machine
{//自定义类:机器
int processing;//正在加工调度方案中的第几个零件,-1表示还未加工零件
int left_time;//剩余加工时间
public:
Machine():processing(-1) , left_time(0) { }
int leftTime(); //返回该机器当前加工零件的剩余加工时间
int process(); //返回正在加工调度方案中的第几个零件
void nextTime(); //进入下一时间状态,可能空闲
void nextWP(int time); //调度方案下一个工件就为
void reset(); //重置状态
};
class Workpiece
{//自定义类:工件
int being_processed;//正在被第几个机器加工,-1表示未被加工
int left_time;//剩余加工时间
public:
Workpiece():being_processed(-1) , left_time(0) { }
int beingPro(); //返回正在被加工的机器号
int leftTime(); //返回目前机器上的剩余加工时间
void nextTime(); //进入下一时间状态,可能处于空闲
void nextMachine(int time); //进入下一机器状态
void reset(); //重置
};
#define Max_number_Machine 30
#define Max_number_Workpiece 60
class simulate
{ //自定义类:模拟整个加工流水线
int cost[Max_number_Workpiece][Max_number_Machine]; //时间花费
int management[Max_number_Machine][Max_number_Workpiece]; //调度方案
int Time[Max_number_Machine][Max_number_Workpiece]; //记录每个机器的调度方案中的零件开始加工时间时间
int Wcount, Mcount; //工件,机器数量
int switchx, switchy , switchrow; //记录交换信息,便于回退至上一状态
Machine M[Max_number_Machine]; //机器集合
Workpiece W[Max_number_Workpiece]; //工件集合
public:
simulate(int n, int m, int* p);
void randomReset(); //随机初始化整个调度方案
void rollBack(); //回退至上一状态
void randomExchange(); //随机交换某两个零件在某个机器上的的加工顺序
int evaluate();//模拟一次加工,返回加工时间,并且维护Time数组
void outdata(); //输出调度方案,时间方案
void resetForNewRound();//重置所有状态,准备下一次模拟
void setManagement(int* p);//设置调度方案,调试时可用
};
各个属性和方法的意义见代码的注释。接下来给出最核心的评估函数函数体:
int simulate::evaluate()
{
int rtime = 0; //模拟时间进行
for (rtime = 0; ; rtime++) {
#ifndef _DEBUG_
std::cout << "time : " << rtime << std::endl;
#endif
for (int i = 0; i < Mcount; i++) {
int processing = M[i].process();//当前正在加工的调度方案的第几个工件
if (processing != Wcount - 1) { //不是最后一个工件
int x = management[i][processing + 1];//获取调度方案下一个工件的序号
if (M[i].leftTime() == 0 && W[x].leftTime() == 0 && W[x].beingPro() == i - 1) {
//下一个工件空闲,并且已加工至当前机器
M[i].nextWP(cost[x][i]);//下一个工件就位
W[x].nextMachine(cost[x][i]);//工件进入下一机器
Time[i][processing + 1] = rtime ; // 记录开始时间
}
}
}
//工件以及机器都进入下一时间状态
for (int i = 0; i < Mcount; i++)
M[i].nextTime();
for (int i = 0; i < Wcount; i++)
W[i].nextTime();
//检验是否加结束加工
int flag = 0;
for (int i = 0; i < Wcount; i++) {
if (W[i].leftTime() == 0 && W[i].beingPro() == Mcount - 1)
continue;
else {
flag = 1;
break;
}
}
if (flag == 0) {
rtime++;
break;//结束
}
}
return rtime;
}
稍微解释一下该函数,在每个时间颗粒的开始,我们先对每一个机器进行上一时间颗粒末状态检测,如果该机器加工剩余时间为0,并且调度方案中的下一个工件正在上一个机器加工且剩余时间为0,那么就将这个工件和机器进入下一机器和下一工件初始状态。检测完成后,所有工件和机器进入下一时间颗粒状态(加工一个时间颗粒),然后检测是否所有工件都加工完(其实代码中的方法可以更简单化,只要检测最后一个机器调度方案中的最后一个工件是否加工完就行)。
其他定义代码见附录链接。1
甘特图
甘特图可以非常形象的给出整个调度方案在加工时,每个工件在机器上的时间分配。我们用的python用上面代码的结果进行画图:
import numpy
import matplotlib.pyplot as plt
mylist = []
labels = []
def getColor(mstr):
"对于不同的工件产生不同的RGB颜色"
num = int(mstr)
num = num / 23
x = num * 256 * 256 * 256
r = x % 256
g = ((x - r) / 256 % 256)
b = ((x - r - g * 256) / (256 * 256) % 256)
r = round(float(r / 256), 2)
g = round(float(g / 256), 2)
b = round(float(b / 256), 2)
return (r , g , b)
inputfile = open("time.txt" , 'r' , encoding='utf-8')
while 1:
templine = inputfile.readline()
if templine == '' :
break
templine = templine.strip('\t\n')
templist = templine.split('\t')
templist1 = []
for iterm in templist :
templist1.append(int(iterm))
mylist.append(templist1)
inputfile.close() #以上部分读入调度方案
inputfile = open("management.txt" , 'r' , encoding='utf-8')
while 1:
templine = inputfile.readline()
if templine == '' :
break
templine = templine.strip('\t\n')
templist = templine.split('\t')
labels.append(templist)
inputfile.close() #以上部分输入每个机器上的每个工件加工时间段
y = numpy.array(mylist)
for i in range(0 , len(mylist)):
for j in range(0 , int(len(mylist[i])/2)):
plt.barh(i + 1 , y[i][2*j + 1] - y[i][2*j] , left= y[i][2*j] , color= getColor(labels[i][j]))
plt.xlabel("time")
plt.ylabel("machine")
plt.show() #画出甘特图
写在最后
保持我上次程序的风格,对于结构比较明显的问题,我喜欢将程序结构化,这样带来的效果是非常明显的,编程的思路会非常明确,先实现各个部分,定义好属性和方法,然后通过一个main函数所在的cpp将这些模块联系起来,这对调试也有很大的帮助。这个程序的逻辑结构比较简单就不画图介绍了,另外所有代码都可以在我的github开源库中得到。1
所有代码在这,开源的. ↩︎ ↩︎ ↩︎