《算法设计与分析》模拟退火算法解决作业车间调度问题

一、作业车间调度问题描述

作业车间调度问题(Job Shop Scheduling, JSP)是最经典的几个NP-hard问题之一。其应用领域极其广泛,涉及航母调度,机场飞机调度,港口码头货船调度,汽车加工流水线等。

JSP问题描述:一个加工系统有M台机器,要求加工N个作业,其中,作业i包含工序数为Li。令,则L为任务集的总工序数。其中,各工序的加工时间已确定,并且每个作业必须按照工序的先后顺序加工。调度的任务是安排所有作业的加工调度排序,约束条件被满足的同时,使性能指标得到优化。

作业车间调度需要考虑如下约束:

Cons1:每道工序在指定的机器上加工,且必须在其前一道工序加工完成后才能开始加工;

Cons2:某一时刻1台机器只能加工1个作业;

Cons3:每个作业只能在1台机器上加工1次;

Cons4:各作业的工序顺序和加工时间已知,不随加工排序的改变而改变。

  • 作业车间调度问题的数学模型

在本课程的综合设计与实现环节中,我们将作业车间调度问题的优化目标设为最大完工时间最小:令(i,j)表示作业i的第j个工序。SijTij分别表示(i,j)的加工起始时刻和加工时间。Zijk表示(i,j)是否在第k台机器上加工:如果(i,j)在第k台机器上加工,Zijk=1;否则,Zijk=0。Ck为第k台机器的完工时间,则问题的数学模型如下:

(1)

(2)

(3)

(4)

公式(1)为目标函数,使最迟完工的机器尽早完成,即加工时间最短;公式(2)表示1个作业只能在加工完成前一道工序后才可以加工后一道工序;公式(3)表示1个作业的第1道工序的起始加工时刻大于或等于0;公式(4)表示在1台机床上不会同时加工1个以上的作业。

三、问题实例

下面给出作业车间调度问题的一个实例,其中每个工序上标注有一对数值(m,p),其中,m表示当前工序必须在第m台机器上进行加工,p表示第m台机器加工当前工序所需要的加工时间。(注:机器和作业的编号从0开始)

  1. jop0=[(0,3),(1,2),(2,2)]
  2. jop1=[(0,2),(2,1),(1,4)]
  3. jop2=[(1,4),(2,3)]

在这个例子中,作业jop0有3道工序:它的第1道工序上标注有(0,3),其表示第1道工序必须在第0台机器上进行加工,且需要3个单位的加工时间;它的第2道工序上标注有(1,2),其表示第2道工序必须在第1台机器上进行加工,且需要2个单位的加工时间;余下的同理。总的来说,这个实例中共有8道工序。

该问题的一个可行解是L=8道工序开始时间的一个排列,且满足问题的约束。下图给出了一个可行解(注:该解不是最优解)的示例:

 

在上图中,我们能看到每个作业的工序按照问题给定的顺序进行了加工,且相互之间没有时间区间重叠。这个可行解的结果是12,即三个作业均被加工完成的时间。

四、算法设计思想

JSP是典型NP-hard问题之一,首先我想解释一下什么是NP-hard问题

在了解NP-hard问题之前,必须了解一个概念叫做多项式时间:在计算复杂度理论中,指的是一个问题的计算时间不大于问题大小的多项式倍数,通俗些理解,我觉得就是一定规模的问题,总有一个时间范围内可以将它解决,这个范围就是多项式时间,然后NP是指非确定性多项式(non-deterministic polynomial,缩写NP)。所谓的非确定性是指,可用一定数量的运算去解决多项式时间内可解决的问题。NP问题就像是如果给你一个确定的例子,你可以很容易的验证他是不是正确,但是如果让你去找这个最优解确实无穷无尽,不太容易求解的。

像这一类在有限个可行解的集合中找出最优解的一类优化问题称为组合优化问题,它是运筹学中的一个重要分支。所研究的问题涉及信息技术、经济管理、工业工程、交通运输、通讯网络等诸多领域。组合优化算法是一类在离散状态下求极值的问题。所以对于这样的组合优化问题,就出现了各种算法作为解决方案。

常用的非数值算法有遗传算法、模拟退火算法和神经网络算法。

这次的作业车间调度问题(Job-shop scheduling problem, JSP)中,得满足约束条件:(1) 每个工件使用每台机器不多于1次;(2) 每个工件利用每台机器的顺序可以不同;(3) 每个工件的工序必须依次加工,后工序不能先于前工序;(4) 任何工件没有抢先加工的优先权,应服从任何生产顺序;(5) 工件加工过程中没有新工件加入,也不临时取消工件的加工。作业车间调度问题需要在有效时间内寻找到最小的加工时间。

鉴于这样的需求分析,我选择的是采用模拟退火算法。模拟退火算法现在运用于数学建模和竞赛较多,以简单有效的搜索方式既避免了数值算法的高计算量,又避免了局部搜索算法快速收敛于局部最优解的缺点。模拟退火算法,需要首先由暴力算法生成一个可行的工序序列,作为退火算法的初始解,构建出解空间,然后在解的空间内找寻命题的最优解。

首先我想介绍一下退火算法,采用模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法来源于固体退火原理。大自然在缓慢降温(亦即,退火)时,可“找到”最低能量状态:结晶。如下图4-1所示,首先(左图)物体处于非晶体状态。我们将固体加温至充分高(中图),再让其徐徐冷却,也就退火(右图)。加温时,固体内部粒子随温升变为无序状,内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小(此时物体以晶体形态呈现)。似乎,大自然知道慢工出细活:缓缓降温,使得物体分子在每一温度时,能够有足够时间找到安顿位置,则逐渐地,到最后可得到最低能态,系统最安稳。

 

图4-1

 

 

模拟退火其实也是一种贪心算法,但是它的搜索过程引入了随机因素。在迭代更新可行解时,以一定的概率来接受一个比当前解要差的解,因此有可能会跳出这个局部的最优解,达到全局的最优解。以下图为例,假定初始解为左边蓝色点A,模拟退火算法会快速搜索到局部最优解B,但在搜索到局部最优解后,不是就此结束,而是会以一定的概率接受到左边的移动。也许经过几次这样的不是局部最优的移动后会到达全局最优点D,于是就跳出了局部最小值。

 

 

图4-2

    

物理退火过程:根据Metropolis(蒙特卡洛)准则,在温度T,分子停留在状态r满足Boltzman概率分布(r是指能量状态,离散化的数据)

 

`E表示分子能量的一个随机变量,E(r)表示状态r的能量,KB>0是Boltzman常数,Z(T)是概率分布的标准化因子,因为概率是一个[0,1]的数值,这就是能量因子的作用,KBT是一个整体,一个温度

物理退火的过程,同一个温度T,选定两个能量E1<E2,有

 

模拟退火算法的基本思想就是:在一定温度下,搜索从一个状态随即地变化到另一个状态,然后这些状态的最小值就是这个温度下的局部最优值,随着温度的不断下降直到最低温度,搜索过程以概率1停留到最优解。

用固体退火模拟组合优化问题,组合优化问题的解可以看做粒子的状态,最优解就是固体能量最低的状态,组合优化问题设定的初始温度便是固体的退火过程,Metropolis抽样过程,就是固体的等温降温的过程,组合优化问题中的控制参数的下降,就是固体的冷却过程,也就是温度T下降:由初始解i和控制参数初值t开始,对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代,并逐步衰减t值,算法终止时的当前解即为所得近似最优解,这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制,包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。

算法先以搜寻空间内一个任意点作起始:每一步先选择一个“邻居”,然后再计算从现有位置到达“邻居”的概率。模拟退火算法是一种通用的优化算法,其物理退火过程由加温过程、等温过程、冷却过程这三部分组成,即:

初始高温 => 温度缓慢下降=> 终止在低温 (这时能量函数达到最小,目标函数最小),对于模拟退火的要求:(1)初始温度足够高,粒子能足够自由运动(2)降温过程足够慢,在每一个温度下都能找到一个最低的状态(3)终止温度足够低,使概率接近于1,虽然可能不是理想状态,但是得达到可接受的范围。

Metropolis准备就是以概率接受新状态,固体在恒定温度下达到热平衡的过程可以用MonteCarlo方法(计算机随机模拟的方法)加以模拟,需大量采样才能得到精确的值。

模拟退火算法也是贪心算法,但是在其过程中引入了随机因素,以一定的概率接受一个比当前解要差的解,并且这个概率随着时间的推移而逐渐降低。

若在温度T,当前状态i—>j

若Ej<Ei,则接受j为当前状态,否则,若概率p大于[0,1)区间的随机数,则仍接受状态j为当前状态,若不成立,则保留状态i为当前状态

在高温下,可接受状态能量差较大的新状态,在低温下,只接受与当前状态能量差较小的状态。

  

模拟退火算法的模型:

模拟退火算法可以分解为解空间、目标函数和初始解三部分。

模拟退火的基本思想:

(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步。

 

流程图如图4-3所示。

 

 

图4-3.算法流程图

 

 模拟退火算法新解的产生和接受可分为如下四个步骤:

  第一步是由一个产生函数从当前解产生一个位于解空间的新解;为便于后续的计算和接受,减少算法耗时,通常选择由当前新解经过简单地变换即可产生新解的方法,如对构成新解的全部或部分元素进行置换、互换等,注意到产生新解的变换方法决定了当前新解的邻域结构,因而对冷却进度表的选取有一定的影响。

  第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生,所以目标函数差的计算最好按增量计算。事实表明,对大多数应用而言,这是计算目标函数差的最快方法。

  第三步是判断新解是否被接受,判断的依据是一个接受准则,最常用的接受准则是Metropo1is准则: 若Δt′<0则接受S′作为新的当前解S,否则以概率exp(-Δt′/T)接受S′作为新的当前解S。

  第四步是当新解被确定接受时,用新解代替当前解,这只需将当前解中对应于产生新解时的变换部分予以实现,同时修正目标函数值即可。此时,当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时,则在原当前解的基础上继续下一轮试验。

模拟退火算法与初始值无关,算法求得的解与初始解状态S(是算法迭代的起点)无关;模拟退火算法具有渐近收敛性,已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法;模拟退火算法具有并行性。

如用模拟退火算法,则首先应该解决编号问题。本程序采用以每个工序为基本变量,对每个工序设置一个工序号进行编码,比如上文的举例,则可以将每个工序编号为12345678,表示jop0的第一个工序、第二个工序、第三个工序、……jop2的第二个工序,这种编码方式的好处在于可以保证工件的加工顺序无误,并且方便快捷。如图4-4所示。

图4-4.编号示意图

 

在已知初始解或者说是第一个解已经生成的情况下,就可以对这个已知的调度进行模拟退火过程,搜索其邻域的解,查看是否有最优解。其中,在SA()函数中主体由外循环和内循环组成。外循环即主要更新参数t,模拟退火过程,由Δ控制设置外循环迭代次数,设置终止温度的阈值;内循环则是用于决定在各温度下产生候选解的数目。内部循环则是采用了前面提到的Metropolis抽样稳定准则,寻找在一定温度下的最优值,在查找过程中,如果找到更优值,也就是比当前已经产生的作业车间调度的解的时间更小的解的话,则直接更新当前最优值,并且记录当前程序的甘特图,如果找到比当前更差的解或者说比当前程序差的解,即生成的作业车间调度时间更长,也以一定概率接受该解,但是,发生这种情况的概率会越来越小,趋近于0。如图4-5所示。

 

图4-5.内外循环示意图

 

五、伪代码

 

退火算法伪代码:

/*

* J(y):在状态y时的评价函数值

* Y(i):表示当前状态

* Y(i+1):表示新的状态

* r: 用于控制降温的快慢

* T: 系统的温度,系统初始应该要处于一个高温的状态

* T_min :温度的下限,若温度T达到T_min,则停止搜索

*/

while( T > T_min )

{

dE = J( Y(i+1) ) - J( Y(i) ) ;

 

if ( dE <=0 ) //表达移动后得到更优解,则总是接受移动

Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动

else

{

// 函数exp( dE/T )的取值范围是(0,1) ,dE/T越大,则exp( dE/T )也大

if ( exp( dE/T ) > random( 0 , 1 ) )

Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动

}

T = r * T ; //降温退火 ,0<r<1 。r越大,降温越慢;r越小,降温越快

/*

* 若r过大,则搜索到全局最优解的可能会较高,但搜索的过程也就较长。若r过小,则搜索的过程会很快,但最终可能会达到一个局部最优值

*/

i ++ ;

}

 

编号程序伪代码:

 

//tmp为编号数组

For i=0; i<=cnt; i++

    tmp[i]=i; //标号为i

初始解伪代码:

 

//tmp[]为所有工序号的全排列

//cnt为所有工序的总和

//machin数组为当前机器上次在工作的时间

//wktime为本任务上次执行的时间

//wktime为本任务上次执行的时间

For I = 1;i<=cnt ;i++

  id = a[tmp[i]].id;//当前任务所属标号

//先判断当前jop执行到了第几个工序

//如果现在排列的工序的前一工序没有执行,直接跳出函数,进行下一次

        if(num[ id ] != a[ tmp[i] ].num)

            return;

//若没有跳出,则

num[id]++;

if(matime[a[ tmp[i] ].machine] < wktime[ a[tmp[i] ].id])//机器在工作

for(int kk=1; kk<=a[ tmp[i] ].time ; kk++)

   if(vis[a[tmp[i]].machine][wktime[a[tmp[i]].id]+kk]==0)//状态数组未被写过,则赋值                   vis[ a[ tmp[i] ].machine ][wktime[a[tmp[i]].id]+kk]=a[tmp[i]].id+1;

            }

            matime[a[ tmp[i] ].machine] = wktime[ a[tmp[i] ].id] = wktime[a[ tmp[i] ].id] +a[ tmp[i] ].time; //机器本次工作的时间记录

        }

        else//机器是空闲的

        {

            for(int kk=1; kk<=a[tmp[i]].time; kk++)

                vis[ a[ tmp[i] ].machine ][ matime[a[ tmp[i] ].machine ] + kk ] = id+1;

            matime[ a[ tmp[i] ].machine ] = wktime[ a[tmp[i] ].id ] =

                                                matime[a[ tmp[i] ].machine] + a[ tmp[i] ].time;//机器本次工作的时间记录

//所有序列组合判断结束,使用temp记录本次结果

temp=-1;

    for(int i=0; i<m; i++)

    {

        temp=max(matime[i],temp);//总时间

}

 

搜索领域解空间伪代码:

 

//使用随机数选则[1,cnt]之间的两个工序

int x = (int)(knum * (rand() / (RAND_MAX + 1.0)));

    int y = (int)(knum * (rand() / (RAND_MAX + 1.0)));

    while(x==y)

    {

        x = (int)(knum * (rand() / (RAND_MAX + 1.0)));

        y = (int)(knum * (rand() / (RAND_MAX + 1.0)));

}

//交换这两个工序

swap(tmp[x],tmp[y]);

//判断这个序列是否合法,合法则执行出结果temp

//不合法则直接跳出

for(int i=0; i<cnt; i++)

    {

        id = a[tmp[i]].id;//当前任务所属标号

        if(num[ id ] != a[ tmp[i] ].num)

            return 0;

        num[id]++;

        if(matime[a[ tmp[i] ].machine] < wktime[ a[tmp[i] ].id])//machine is working

        {

            for(int kk=1; kk<=a[ tmp[i] ].time ; kk++)

            {

                if(vis[a[tmp[i]].machine][wktime[a[tmp[i]].id]+kk]==0)

                    vis[ a[ tmp[i] ].machine ][wktime[a[tmp[i]].id]+kk]=a[tmp[i]].id+1;

            }

            matime[a[ tmp[i] ].machine] = wktime[ a[tmp[i] ].id] = wktime[a[ tmp[i] ].id]

                                          +a[ tmp[i] ].time;

        }

        else//mathine is waiting

        {

            for(int kk=1; kk<=a[tmp[i]].time; kk++)

            {

                vis[ a[ tmp[i] ].machine ][ matime[a[ tmp[i] ].machine ] + kk ] = id+1;

            }

            matime[ a[ tmp[i] ].machine ] = wktime[ a[tmp[i] ].id ] =

                                                matime[a[ tmp[i] ].machine] + a[ tmp[i] ].time;

        }

    }

 

作业车间调度的退火算法伪代码:

//初始温度T

double t = T;

    srand(time(NULL));

    int CurCost = temp;

    int NewCost = temp;

    int P_L = 0;

    int P_F = 0;

    while(1) //外循环,主要更新参数t,模拟退火过程

    {

        for(int i=0; i<OLOOP; i++)//内循环,寻找在一定温度下的最优值

        {

//产生一个可行解

            while(!GetNext())

            {

                ;

            }

            NewCost = temp;//为新解赋值

            double dE = NewCost - CurCost;

            if(dE<0)//找到更优解,更新

            {

                CurCost = NewCost;

                ans = temp;

                //更新甘特图数组

                for(int i=0; i<3; i++)

                {

                    for(int j=0; j<=ans; j++)

                    {

                        gt[i][j]=vis[i][j];//保存最优态

                    }

                }

                P_L = 0;

                P_F = 0;

            }

            else

            {

                double rd = rand() / (RAND_MAX +

                                      1.0);    //如果找到比当前更差的解,以一定概率接受该解,并且这个概率会越来越小

                if(exp(dE/t)>rd&&exp(dE/t)<1)

                {

                    CurCost = NewCost;

                }

                P_L++;

            }

//内循环超限

            if(P_L > LIMIT)

            {

                P_F++;

                break;

            }

        }

        //内循环条件判断,择机跳出退火过程

        if(P_F > OLOOP || t < EPS)

            break;

        t *= DELTA;

}

 

六、程序代码及运行结果

程序代码:

#include <bits/stdc++.h>

#define mme(i,j) memset(i,j,sizeof(i))//新申请的内存做初始化工作

#define ll long long int

#define maxs 202020

#define NUM 110

#define T 3000//系统的初始温度,需要达到一个高温的状态,保证粒子可以足够自由运动,放在组合优化中,就是我们的解的覆盖性尽可能全

#define EPS 1e-8

#define DELTA 0.98//温度衰减率,控制降温的快慢

//若r过大,则搜索到全局最优解的可能会较高,但搜索的过程也就较长。若r过小,则搜索的过程会很快,但最终可能会达到一个局部最优值

//所以一般的选取就是0.95—-0.99

#define LIMIT 10000//概率选择上限

#define OLOOP 1000//内循环

#define ILOOP 15000//外循环

 

using namespace std;

 

int tmp[20];//所有工序号的全排列

int n,m;

int vis[10][20],gt[15][20];

int k[10];

int num[6];

int matime[6];//数组为当前机器上次在工作的时间

int wktime[6];//为本任务上次执行的时间

int cnt=0;

int ans = INT_MAX;

int knum = 0;

int temp;

int flag = 0;

 

struct work

{

    int id;//工作号

    int machine;

    int time;

    int num;//工序号

} a[100],b[100];

 

void Init()//定义初始解空间,解决编号问题,需要通过赋值写入的方式

{

    mme(matime,0);

    mme(num,0);

    mme(wktime,0);

    mme(vis,0);

    int id;

    for(int i=0; i<cnt; i++)

        //cnt为工序的总和

        //wktime为本任务上次执行的时间

        //wktime为本任务上次执行的时间

    {

        id = a[tmp[i]].id;//当前任务所属标号

        //如果现在排列的工序的前一个工序没有执行,直接跳出函数

        if(num[ id ] != a[ tmp[i] ].num)

            return;

            //如果符合,则

        num[id]++;//表示可进入目前工序的加工,用来下一次判断

        if(matime[a[ tmp[i] ].machine] < wktime[ a[tmp[i] ].id])//机器在工作

        {

            for(int kk=1; kk<=a[ tmp[i] ].time ; kk++)

            {

                if(vis[a[tmp[i]].machine][wktime[a[tmp[i]].id]+kk]==0)//状态数组未被写过,则赋值

                    vis[ a[ tmp[i] ].machine ][wktime[a[tmp[i]].id]+kk]=a[tmp[i]].id+1;

            }

            matime[a[ tmp[i] ].machine] = wktime[ a[tmp[i] ].id] = wktime[a[ tmp[i] ].id]

                                          +a[ tmp[i] ].time;// //机器本次工作的时间记录

        }

        else//机器是空闲的

        {

            for(int kk=1; kk<=a[tmp[i]].time; kk++)

            {

                vis[ a[ tmp[i] ].machine ][ matime[a[ tmp[i] ].machine ] + kk ] = id+1;

            }

            matime[ a[ tmp[i] ].machine ] = wktime[ a[tmp[i] ].id ] =

                                                matime[a[ tmp[i] ].machine] + a[ tmp[i] ].time;

        }

    }

//所有序列组合判断结束,使用temp记录本次结果

    temp=-1;

    for(int i=0; i<m; i++)

    {

        temp=max(matime[i],temp);//总时间

    }

    if(temp<ans)

    {

        ans=temp;

        mme(gt,0);

        for(int i=0; i<3; i++)

        {

            for(int j=0; j<=ans; j++)

            {

                gt[i][j]=vis[i][j];//保存最优态

            }

        }

        flag = 1;

    }

}

//搜索邻域解空间,这也就是刚才讲到的第一步,由当前解构造出新解,交换工序,之后就得在判断是不是合法的序列

int GetNext()

{//使用随机数选择[1,cnt]之间的两个工序,

    int x = (int)(knum * (rand() / (RAND_MAX + 1.0)));

    int y = (int)(knum * (rand() / (RAND_MAX + 1.0)));

    while(x==y)

    {

        x = (int)(knum * (rand() / (RAND_MAX + 1.0)));

        y = (int)(knum * (rand() / (RAND_MAX + 1.0)));

    }

     swap(tmp[x],tmp[y]);

    mme(matime,0);

    mme(num,0);

    mme(wktime,0);

    mme(vis,0);

    //判断这个序列是否合法,合法则执行出结果temp,合法直接退出

    int id;

    for(int i=0; i<cnt; i++)

    {

        id = a[tmp[i]].id;//当前任务所属标号

        if(num[ id ] != a[ tmp[i] ].num)

            return 0;

        num[id]++;

        if(matime[a[ tmp[i] ].machine] < wktime[ a[tmp[i] ].id])//machine is working

        {

            for(int kk=1; kk<=a[ tmp[i] ].time ; kk++)

            {

                if(vis[a[tmp[i]].machine][wktime[a[tmp[i]].id]+kk]==0)

                    vis[ a[ tmp[i] ].machine ][wktime[a[tmp[i]].id]+kk]=a[tmp[i]].id+1;

            }

            matime[a[ tmp[i] ].machine] = wktime[ a[tmp[i] ].id] = wktime[a[ tmp[i] ].id]

                                          +a[ tmp[i] ].time;

        }

        else//mathine is waiting

        {

            for(int kk=1; kk<=a[tmp[i]].time; kk++)

            {

                vis[ a[ tmp[i] ].machine ][ matime[a[ tmp[i] ].machine ] + kk ] = id+1;

            }

            matime[ a[ tmp[i] ].machine ] = wktime[ a[tmp[i] ].id ] =

                                                matime[a[ tmp[i] ].machine] + a[ tmp[i] ].time;

        }

    }

    temp=-1;

    for(int i=0; i<m; i++)

    {

        temp=max(matime[i],temp);//总时间

    }

    return temp;

}

//退火部分

int SA()

{

    double t = T;//初始温度T

    srand(time(NULL));

    int CurCost = temp;

    int NewCost = temp;

    int P_L = 0;

    int P_F = 0;

    while(1)//外循环,主要更新参数t,

        for(int i=0; i<OLOOP; i++)//内循环,寻找一定温度的最优值,内循环则是用于决定在各温度下产生候选解的数目。

        {

//产生一个可行解

            while(!GetNext())

            {

                ;

            }

            NewCost = temp;//为新解赋值

            double dE = NewCost - CurCost;//公式中的变化量

            if(dE<0)//表达移动后得到更优解,则总是接受移动,如果找到更优解,直接更新//计算增量的好处

            {

                CurCost = NewCost;

                ans = temp;

                //更新甘特图

                for(int i=0; i<3; i++)

                {

                    for(int j=0; j<=ans; j++)

                    {

                        gt[i][j]=vis[i][j];//保存最优态

                    }

                }

                P_L = 0;

                P_F = 0;

            }

            else

            {

                double rd = rand() / (RAND_MAX +

                                      1.0);    //如果找到比当前更差的解,以一定概率接受该解,并且这个概率会越来越小

                if(exp(dE/t)>rd&&exp(dE/t)<1)

                {

                    CurCost = NewCost;

                }

                P_L++;

            }

            //内循环超限

            if(P_L > LIMIT)

            {

                P_F++;

                break;

            }

        }

        //内循环条件判断,择机跳出退火过程

        if(P_F > OLOOP || t < EPS)

            break;

        t *= DELTA;

    }

}

 

int main()

{

    while(~scanf("%d%d",&n,&m))

    {

        for(int i=0; i<n; i++)

        {

            scanf("%d",&k[i]);

            knum+=k[i];

            for(int j=0; j<k[i]; j++)

            {

                scanf("%d%d",&a[cnt].machine,&a[cnt].time);

                a[cnt].id=i;//属于哪个工作

                a[cnt].num=j;//第几个工序

                cnt++;

            }

        }

        for(int i=0; i<=cnt; i++)

            tmp[i]=i;

        do

        {

            Init();

            if(flag)break;

        }

        while(next_permutation(tmp,tmp+cnt));//做全排列

        SA();

        printf("The answer is %d\n",ans);

        for(int i=0; i<3; i++)

        {

            for(int j=1; j<=ans; j++)

            {

                printf("%3d",gt[i][j]);

            }

            puts("");

        }

    }

}

 

 

运行结果:

 

甘特图:

 

 

七、总结

通过这一次的作业车间调度问题,让我体会到NP-hard问题的解题难度,在刚拿到这个题目理解这个NP-hard的时候,参考了著名的推销员收获问题,假设一个推销员需要从香港出发,经过广州,北京,上海,…,等 n 个城市, 最后返回香港。 任意两个城市之间都有飞机直达,但票价不等。假设公司只给报销 C 元钱,问是否存在一个行程安排,使得他能遍历所有城市,而且总的路费小于 C?仔细分析此问题和作业车间调度问题的异同点,在知网上搜索了解决作业车间调度问题的方法,论文大部分是采用的遗传算法,但是在变异部分的优化方法,有点想不通,然后所以就想到了竞赛经常用到的模拟退火算法,杭电上的求最小外接球半径,均分数据的题,这个算法的难度会比遗传算法简单点。它能以简单有效的搜索方式既避免了数值算法的高计算量,又避免了局部搜索算法快速收敛于局部最优解的缺点。在代码实现的过程中,我也遇到了各种产生的错误,所以又在知网上和CSDN上搜索讲义,反复修改代码,直到最终的运行出正确的结果,在这其中收获了许多。

首先是在分析NP-hard问题上,在分析问题的过程中,搜索了很多算法,比如分支限界法,遗传算法,模拟退火算法,像遗传算法就是之前未接触过的,了解了这些算法的基本定义,以及区分出了一个基本的难易度,再落回分析jsp问题上,刚拿到这个问题的时候感觉很无力,只能暴力出解空间的具体值,无从下手,最后将拿到大的问题逐步分解,使原本看起来难解的一个题目变成一个可以解决的易解问题,再重新逐步推演到原本的大的问题。

然后就是代码实现的问题,优化代码的问题的过程中,遇到了各种各样的问题,尤其是改bug的过程中,有时候越改越出现的连锁的问题,所以就需要耐心和更多的知识去一起更改。

这一次综合设计作业车间调度的程序,仍有一些需要优化的地方,比如可以从当前的工序中挑出最长的作业序列,然后再在其领域进行搜索,或许会得到一个运行时间更好的程序,或时记录每次的序列当序列存在时即采用已经出现的解决值,均可对程序进行进一步的优化。而且对于我本身来说,还需要更加认真的学习多一些的算法,比如下次争取使用遗传算法解决NP-hard算法,实际掌握新的算法。

 

  • 15
    点赞
  • 113
    收藏
    觉得还不错? 一键收藏
  • 20
    评论
流水车间调度问题是指在一个车间中有多个工序需要完成,每个工序需要在一台特定机器上完成,不同工序之间存在先后顺序的限制,同时每个机器只能同时处理一个工序。流水车间调度问题是一个NP难问题,难以通过精确求解算法得到最优解。模拟退火算法是一种常用的启发式算法,被广泛应用于求解NP难问题,包括流水车间调度问题模拟退火算法的基本思想是通过模拟物理系统的退火过程来寻找问题的全局最优解。在流水车间调度问题中,可以将不同的工序看作组成系统的粒子,每个粒子的状态对应着一种调度方案。算法从一个随机的初始状态开始,通过一系列随机变换(如交换两个工序的位置)来寻找更优的解,并以一定概率接受劣解,以避免陷入局部最优解。 具体地,模拟退火算法通过定义一个能量函数来衡量当前调度方案的质量,通常采用流水车间调度问题中的加权完工时间(weighted completion time)作为能量函数。算法从一个随机的初始状态开始,通过一系列随机变换(如交换两个工序的位置)来寻找更优的解,并以一定概率接受劣解,以避免陷入局部最优解。随着算法迭代的进行,接受劣解的概率逐渐减小,直到最终收敛于一个较优的解。 需要注意的是,模拟退火算法的求解结果并不能保证是最优解,但可以得到一个较优的近似解。同时,算法的求解时间也比较长,需要通过参数调优和合理的实现方式来提高求解效率。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 20
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值