0 引言

车间调度问题广泛存在于钢铁生产、冶金、炼油、化工、制药等生产领域^[1-2],在某些工艺约束下,要求工件的加工不允许等待,即在当前机器上完成加工后必须立即进入下一台机器加工,这样的约束称为零等待约束,相应的流水车间调度问题称为零等待流水车间调度问题(no-wait flow shop scheduling problem,NWFSP)。

NWFSP已被证明在复杂度理论上是NP-难的,由于该问题的复杂性和在工业实际中的广泛存在性,NWFSP引起了众多研究者的关注。较早的工作是采用启发式规则求解^{[1, 3-4]},而近年来更多的研究者采用智能方法进行求解。文献^[5]提出了混合模拟退火的遗传算法和变邻域搜索两种方法进行求解;文献^[6]提出了下降搜索方法和禁忌算法;文献^[7-9]设计了微粒群算法求解该问题并取得了显著效果;其他的方法还有离散差分进化算法^[10]、离散自组织迁移算法^[11]、迭代贪婪算法^[12]、蚁群算法^[13]等。上述方法以最大完工时间为最小化目标进行求解,而对总流程时间的最小化,文献^[14-15]提出了一种离散和声搜索算法与一种混合和声搜索算法;文献^[16]通过嵌入局部搜索技巧,提出了一种新型微粒群算法求解;文献^[17-19]提出了几种构造性启发式规则,取得了不错的效果。

上述研究都是对于单一目标进行调度优化,而针对多目标NWFSP(multi-objective NWFSP,MNWFSP)的研究目前还十分有限。文献^[20]以最大完成时间和最大推迟完成时间为目标,提出了一些混合算法;文献^[21]以加权平均完成时间和加权平均拖期为目标,提出了一种混合多目标免疫算法;文献^[22]引入三种局部搜索方法,以最大完工时间和最大延迟时间为目标,提出了多目标混合差分进化算法;针对同样的多目标,文献^[23]提出了离散编码的离散差分进化算法进行求解。

针对多目标优化问题,目前主要有先验法和后验法两种多目标处理方法。先验法需要决策者在优化的初始阶段即给出多目标相关的必要信息,如各目标的权重或各目标的优先级,然后采用单目标方法求解;后验法为决策者提供有效的解集合,通常表现为Pareto最优解或非支配解集。本研究考虑MNWFSP,提出一种离散差分进化算法进行Pareto求解。

作为一种智能优化方法,差分进化算法自提出后便不断得到改进,同时,也逐渐被应用到其他各种领域中^[24-25]。在车间调度领域,差分进化算法由于其简单高效的并行搜索机制而被广泛地应用^{[10, 22-23]}。对多目标调度而言,由于各个目标相互制约,Pereto前沿的分布将非常不规则,从而导致采用单一的邻域搜索结构难以取得满意效果。因此,本研究考虑引入变邻域搜索方法,在差分进化操作完成后,算法在当前非支配解的基础上进行变邻域搜索,从而提出一种多目标变邻域搜索差分进化混合算法MDEVNS,并将其用于求解以最大完工时间和总流程时间为目标的MNWFSP。

1 MNWFSP 1.1 NWFSP问题描述

NWFSP研究具有n个工件和m个机器的流水加工过程,每个工件j(j=1,2,…,n)必须在第k(k=1,2,…,m)台机器上加工,且加工次序必须是机器1,机器2,直到机器m。工件j在机器k上的加工过程记为工序O_jk,相应的加工时间已知为p(j,k)。零等待约束存在于该流水车间中,因此,每个工件一旦开始加工,则必须不停歇地完成所有工序的加工,每个工序(O_jk)的完成时刻必须等于下一工序(O_j,k+1)的开始时刻。显然,零等待的约束很可能导致工件的开始加工时间后延,从而导致整个生产周期的增大。

为使问题描述严格,还需做以下假设:(1) 任意时刻,一台机器只能加工至多一个工件,且一个工件只能在至多一台机器上加工;(2) 所有工件和机器在零时刻是可用的;(3) 工件不允许拆分,必须整体加工;(4) 加工时间已知且大于零;(5) 工件的准备时间、在机器间的转移时间均忽略不计。

1.2 MNWFSP 模型

零等待约束必然要求各台机器上的工件先后次序是一致的,即意味着是置换流水车间。因此,若确定一个工件的排列次序,则每台机器按该次序加工,可得到一个可行调度。显然,工件排列不同,则最大完工时间也可能不同。由于总流程时间也是一个很重要的调度指标,这里以总流程时间和最大完工时间为目标,寻找Pareto最优解集。

下面给出Pareto 优化的几个重要概念。考虑多目标优化模型:

$Min\left\{ {{f}_{1}}(x) \right.,{{f}_{2}}(x),\ldots ,\left. {{f}_{g}}(x) \right\},subjected\text{ }to\text{ }x\in X,$

(1)

其中: f₁(x),f₂(x),…,f_g(x)是需要最小化的g个目标; X是可行解集。

(1) Pareto 支配:当且仅当以下两个条件成立时,称为解a(a∈X)支配解b(b∈X),记为a$\prec $b:

${{f}_{i}}(a)\le {{f}_{i}}(b),\forall i\in \left\{ 1,2,\ldots ,g \right\}$

(2)

${{f}_{j}}(a)＜{{f}_{j}}(b),\exists j\in \left\{ 1,2,\ldots ,g \right\}$

(3)

(2) Pareto 最优解:对解x∈X,若不存在x′∈X且x′$\prec $x,则称x是Pareto 最优解,也称非支配解。

(3) 合并$\uplus $:将两个解集X₁和X₂合并,得到的解集X仅保留X₁和X₂集合的非支配解,记为X=X₁$\uplus $X₂。

对MNWFSP ,其可行解采用工件排列π={π₁,π₂,…,π_n}表示,其中π_j表示工件号(j=1,2,…,n)。用d(π_j-1,π_j)表示第一台机器上的工件π_j-1和π_j开始时刻的最小延迟,C(π_j,m)表示π_j在机器m上的完成时刻,f₁表示最大完工时间目标,f₂表示总流程时间目标,则d(π_j-1,π_j)可通过式(4)计算:

$d({{\pi }_{j-1}},{{\pi }_{j}})=p({{\pi }_{j-1}},1)+max\left[ 0,\underset{2\le k\le m}{\mathop{max}}\,\left\{ \sum\limits_{h=2}^{k}{p({{\pi }_{j-1}},h)}-\sum\limits_{h=1}^{k-1}{p({{\pi }_{j}},h)} \right\} \right]$

(4)

在此基础上可计算

$C({{\pi }_{j}},m)=\sum\limits_{i=2}^{j}{d({{\pi }_{i-1}},{{\pi }_{i}})}+\sum\limits_{k=1}^{m}{p({{\pi }_{j}},k)},$

(5)

于是目标函数可计算

${{f}_{1}}(\pi )=C({{\pi }_{n}},m)$

(6)

${{f}_{2}}(\pi )=\sum\limits_{j=1}^{n}{C({{\pi }_{j}},m)}$

(7)

令所有的排列π构成集合Π,则问题为

$Min\left\{ {{f}_{1}}(\pi ),{{f}_{2}}(\pi ) \right\},subjected\text{ }to~\pi \in \Pi $

(8)

为说明MNWFSP 的目标函数计算,下面给出一个4个工件3台机器的例子。所有的加工时间p(j,k) 构成时间矩阵$P=\left( \begin{matrix} 2 & 3 & 5 \\ 5 & 1 & 3 \\ 4 & 5 & 2 \\ 1 & 3 & 2 \\ \end{matrix} \right)$,若工件排列为π={1,2,3,4},则调度甘特图如图 1所示,目标函数计算过程如下:

$\begin{align} & d\left( 1,2 \right)=p\left( 1,1 \right)+max\text{ }\left[ 0 \right.,max\left\{ p\left( 1,2 \right)-p\left( 2,1 \right),p\left( 1,2 \right)+p\left( 1,3 \right)- \right. \\ & \left. \left. \left( p\left( 2,1 \right)+p\left( 2,2 \right) \right) \right\} \right]=4, \\ & d\left( 2,3 \right)=p\left( 2,1 \right)+max\text{ }\left[ 0,max\left\{ p\left( 2,2 \right) \right.-p\left( 3,1 \right),p\left( 2,2 \right)+p\left( 2,3 \right)- \right. \\ & \left. \left. \left( p\left( 3,1 \right)+p\left( 3,2 \right) \right) \right\} \right]=5, \\ & d\left( 3,4 \right)=p\left( 3,1 \right)+max\left[ 0,maxp\left( 3,2 \right)-p\left( 4,1 \right),p\left( 3,2 \right)+p\left( 3,3 \right)- \right. \\ & \left. \left. \left( p\left( 4,1 \right)+p\left( 4,2 \right) \right) \right\} \right]=8, \\ & C\left( 1,3 \right)=\sum\limits_{k=1}^{3}{p\left( 1,k \right)}=10,C\left( 2,3 \right)=d\left( 1,2 \right)+\sum\limits_{k=1}^{3}{p\left( 2k \right)}=13, \\ & C\left( 3,3 \right)=d\left( 1,2 \right)+d\left( 2,3 \right)+\sum\limits_{k=1}^{3}{p\left( 3,k \right)}=20, \\ & {{f}_{1}}\left( \pi \right)=C\left( 4,3 \right)=d\left( 1,2 \right)+d\left( 2,3 \right)+d\left( 3,4 \right)+\sum\limits_{k=1}^{3}{p\left( 4,k \right)}=23, \\ & {{f}_{2}}\left( \pi \right)=\sum\limits_{j=1}^{4}{C\left( {{\pi }_{j}},3 \right)}=66 \\ \end{align}$

2 MDEVNS

差分进化算法种群内的个体称为目标个体,在每次迭代中,算法最初通过变异操作产生变异个体,再通过交叉操作产生试验个体,最后通过选择操作对目标个体更新。本研究采用离散的工件排列编码这一直接的编码方式,又由于本研究需要进行多目标寻优,因此差分进化的变异、交叉、选择机制需要进行相应设计。此外，在算法中混合了多目标变邻域搜索以增强算法的多目标搜索性能。

2.1 种群初始化

作为性能优异的启发式规则,Nawaz-Enscore-Ham(NEH)规则被广泛应用在流水车间调度问题中。其思想为首先按工件在机器上的总加工时间递减的顺序得到n个工件的排列,然后依次将该排列中的第k(k=2,…,n) 个工件插入到前面的最优位置,直到所有工件调度完成。该规则是针对最大完工时间目标提出来的,而针对总流程时间目标,WANG L等对NEH第一步中的工件按照递增排序,得到了更优的NEH_-WPT规则^[26]。根据初步的计算试验,将该步骤中工件随机排列,得到的调度解多样性更好,且质量并无明显差异,本研究将该方法记为规则NEH_-RAN^[27]。

MDEVNS算法在每一次的迭代中始终维持大小为ps的目标种群,记为PL。在算法进入迭代操作之前,需初始化种群中的每个个体。本研究设计以下多样化种群初始方式:一个初始个体采用以f₁ 为目标的NEH产生,一个初始个体采用以f₂ 为目标的NEH_-WPT产生,其他ps-2个个体采用以f₃ 为目标的NEH_-RAN产生。其中,NEH_-RAN过程中,目标函数f₃做如下设计:

${{f}_{3}}=\frac{k}{ps-1}{{f}_{1}}+\frac{ps-k-1}{ps-1}{{f}_{2}},k=1,2,\ldots ,ps-2$

(9)

这样,初始种群的目标函数能较均匀地分布在f₁和f₂的二维平面,从而使得初始种群具有较好性能的同时兼具多样性。

由于算法进行多目标搜索,因此算法始终维持一个非支配解集合,记为 NDS=｛S₁,S₂,…,S_nb｝,nb为当前非支配解集合的大小。在初始化种群之后,使用PL的非支配集来初始化NDS。为配合后续的局部搜索方法,NDS中的每个解都初始化标记为“未搜索”,而在后续搜索中,该标记可能改变为“已搜索”。

2.2 差分进化操作

设第t次迭代种群中的第i个目标个体、变异个体和试验个体分别用π_i^t,V_i^t和U_i^t表示。

变异个体通过插入操作产生:

${{V}_{i}}^{t}=\left\{ \begin{matrix} insert\text{ }({{S}_{r1}}), & if rand ＜{{p}_{m}}, \\ insert\text{ }\left( \pi _{r2}^{t-1} \right), & otherwise, \\ \end{matrix} \right.$

(10)

其中: insert (·)表示一次随机插入移动;p_m 是变异率,表示变异个体通过NDS 产生的概率;r1是[1,nb ]间的随机整数;r2是[1,ps ]间的随机整数且r≠i; rand 是(0,1)范围内的随机数。

试验个体通过交叉操作产生:

${{U}_{i}}^{t}=\left\{ \begin{matrix} crossover\text{ }({{V}_{i}}^{t},\pi _{i}^{t-1}), & if\text{ }rand\text{ }{{p}_{c}}, \\ {{V}_{i}}^{t}, & otherwise \\ \end{matrix} \right.$

(11)

其中: crossover (·)表示对两个个体进行交叉操作,并随机返回一个个体; p_c 是交叉率,表示试验个体通过交叉产生的概率;交叉操作选用部分匹配交叉。

选择操作根据试验个体和目标个体的支配关系,并结合当前NDS 进行设计,方式为:

$\pi _{i}^{t}=\left\{ \begin{matrix} U_{i}^{t}, & if U_{i}^{t}\prec \pi _{i}^{t-1}, \\ \pi _{i}^{t-1}, & if \pi _{i}^{t-1}\prec U_{i}^{t}, \\ \Phi , & otherwise, \\ \end{matrix} \right.$

(12)

其中,Φ表示根据以下情况决定目标个体:

情况1: U_i^t被NDS 中的某个解支配,若π^t-1_i被NDS 中的某个解支配,则Φ=U_i^t;若π^t-1_i不被NDS 中的任何解支配,则Φ=π^t-1_i。

情况2: U_i^t不被NDS 中的任何解支配,若U_i^t和NDS 中的某个解相同,则计算U_i^t和π^t-1_i的拥挤距离,选择拥挤距离较大者作为Φ;若U_i^t不同于NDS 中任何解,则Φ=U_i^t,并且NDS=NDS $\uplus $U_i^t。

注意,式(12)中U_i^t$\prec $π^t-1_i的情况下,需要执行NDS=NDS $\uplus $U_i^t,且若U_i^t加入了NDS,应立即将其标记为“未搜索”。上述差分进化操作可保持目标个体不断地向Pareto最优解进化,且种群多样性亦能较好地维持。

2.3 多目标变邻域搜索

在给出多目标变邻域搜索(multi-objective variable neighborhood search,MVNS)方法之前,这里先给出对插入邻域的局部搜索过程,记为基于插入的Pareto局部搜索(insert-based Pareto local search,IPLS)。其步骤为:

Step 1 产生一个随机的排列πr=｛πr₁,πr₂,…,πr_n｝,令i=0,j=1;

Step 2 找出x中产品πr_j的位置,将x中产品πr_j插入到x的其他n-1个位置上,可得n-1个解,通过计算目标值可得这n-1个解的局部非支配集合 (local non-dominated solution set),记为LNDS,若存在x′∈ LNDS且x′$\prec $x,则 LNDS=LNDS＼x′,x=x′并置i=1;否则,i=i+1。置j=(j+1) ％ n;

Step 3 NDS=NDS$\uplus $LNDS,若NDS中有新解加入,则将其标记为“未搜索”;

Step 4 若i ＜ n,则转 Step 2;否则,NDS=NDS$\uplus $x,若x加入 NDS中,则标记为“已搜索”。

将IPLS的邻域结构改为交换邻域,则可得到基于交换的Pareto局部搜索(swap-based Pareto local search,SPLS)。考虑采用多个邻域结构有利于搜索多目标最优,因此,本研究设计以下多目标变邻域搜索,通过交替使用IPLS和SPLS,不断搜索解空间,直到不再有更优解出现,其步骤为:

Step 1 对x执行 IPLS;

Step 2 对x执行 SPLS,若x未更新,则结束;否则,转 Step 3;

Step 3 对x执行 IPLS,若x未更新,则结束;否则,转 Step 2。

上述多目标变邻域搜索在算法的每次迭代中对x执行一次。解x的选择方法为:若当前 NDS中存在“未搜索”的解,则随机选择一个“未搜索”的解为x;若不存在“未搜索”的解,则随机选择一个 NDS中的解,对该解进行d次随机插入移动,所得解作为x。

2.4 MDEVNS算法流程

MDEVNS算法流程如图 3所示。算法首先利用NEH、NEH_-WPT、NEH_-RAN的多样化种群初始方法产生初始种群PL,并初始化非支配解集NDS,然后进入迭代。在每一代中,算法首先执行差分进化的变异、交叉和选择操作,然后对NDS中的某个个体执行多目标变邻域搜索MVNS。

3 计算研究 3.1 算法运行设置

为了验证MDEVNS算法的效果,选取了各种规模下的Taillard算例进行测试。该算例集包含工件数为20、50、100、200、500,机器数为5、10、20的情形,共包括120个算例,本研究对每种规模选取1个算例,得到Ta10、Ta20、…、Ta120共12个算例,基本涵盖了实际工业应用上可能遇到的问题规模。涉及的算法均采用C++语言编程实现,运行环境是4 GB内存的Intel(R) Core(TM) i7-2600 CPU @3.06 GHz PC。MDEVNS算法4个参数设置为:种群大小ps=15,变异概率p_m=0.4,交叉概率p _c=0.8,扰动规模d =5。算法运行终止条件为CPU时间50 mn ms。

对多目标优化算法而言,有多种性能评价准则,这里使用Coello和Cortés给出的Inverted Generational Distance (IGD),IGD定义如下:

假设P^*是已知的某个参考解集,A是某算法找到的非支配解集。解x∈A和解y∈P^*的归一化欧氏距离为

$d\left( x,y \right)=\sum\limits_{i=1}^{2}{{{\left( \frac{{{f}_{i}}(x)-{{f}_{i}}(y)}{f_{i}^{\max }-f_{i}^{\min }} \right)}^{2}},}$

(13)

其中: f_i(·)是解的第i个目标函数值; f_i ^max、 f_i ^min分别是P^*中该目标函数的最大值和最小值。

在此基础上,可计算A的 IGD值为

$IGD(A,{{P}^{*}})=\frac{1}{\left| {{P}^{*}} \right|}\sum\limits_{y\in {{P}^{*}}}{\underset{x\in A\text{ }}{\mathop{min}}\,}~d\left( x,y \right),$

(14)

式中: IGD是A到P^*的所有解的最短距离的平均,该指标可反映A的多样性和收敛性。

3.2 算法性能比较

本研究将MDEVNS算法和以下两种较先进的群智能算法进行比较:多目标多初始点模拟退火(bi-objective multi-start simulated annealing,BMSA)算法^[28]和非支配排序遗传算法(non-dominated sorting genetic algorithm,NSGA-Ⅱ)^[29]。每种参与比较的算法对考虑的每个算例均运行10次,每次运行可得到1个非支配解集。

完成计算试验后,需要对结果进行处理以比较算法间的性能差异。为了计算方便,对每个算例,将各个算法10次运行得到的非支配解集合并,作为算法获得的非支配解集 A,并将所有运行结果的非支配解集合并,作为参考解集P^*。以此为基础,计算IGD值,并按照算例列于表 1,表中加粗字体表示最小的IGD值。

分析表 1可知,对工件数为20的较小规模算例,3个算法获得的结果都等于或非常接近0(仅NSGA-Ⅱ算法获得了0.03的非零值),这意味着各个算法的非支配解集和参考解集在收敛性和多样性上相同或非常接近。从数据上看,似乎3个算法的求解效果在该问题规模上并无明显差别,然而,这一现象也很可能由算例难度不大的因素引起。对工件数为50、100、200、500的中大规模算例,MDEVNS算法获得的IGD值要明显小于BMSA和NSGA-Ⅱ算法,例如,对Ta100,MDEVNS算法获得的IGD值为0.01,而BMSA和NSGA-Ⅱ算法的IGD值分别为0.93和3.06,可见,MDEVNS算法对中大规模算例的优势逐渐明显。总体上看,MDEVNS算法的效果明显优于BMSA和NSGA-Ⅱ算法。表 1还从IGD值的角度说明了BMSA和NSGA-Ⅱ算法的性能差异。对算例Ta10-Ta50,两者的IGD值差异较小,而对算例Ta60-Ta120,BMSA算法获得的IGD值要明显小于NSGA-Ⅱ算法,可见,随着算例规模的增大,相比于NSGA-Ⅱ算法,BMSA算法的性能优势越来越明显。

为从直观的角度说明算法间的性能差异,本研究选取Ta50、Ta90两个算例,将各个算法随机运行1次,并将获得的非支配解集分别示于图 4和图 5中。分析图 4和图 5可知,MDEVNS算法获得的非支配解在分布上具有较好的多样性,特别是解的质量要明显优于另两种算法。

MDEVNS算法性能表现较好的原因为:一方面,差分进化算法的初始化、变异个体、试验个体、目标个体等操作设计适当,可以维持种群的多样性和算法搜索的广度;另一方面,算法中嵌入的多目标变邻域搜索过程可持续搜索邻域结构内的解,具有较好的搜索深度。正是通过这两方面的设计,MDEVNS算法在开发性和探索性上具有较好的平衡,因此获得了优于另两种算法的求解效果。

4 总结

针对存在零等待约束的流水车间调度问题,以最大完工时间和总流程时间为双目标,给出了多目标调度模型,设计了一种多目标变邻域搜索差分进化混合算法进行Pareto求解。提出的算法采用基于启发式规则的多样性种群初始化策略,并采用适用于多目标优化的差分进化操作。特别地,该算法还混合了一种多目标变邻域搜索方法,搜索调度解的插入和交换邻域,以增强算法的多目标寻优能力。对各种规模下的算例试验表明,提出的算法能较好地解决零等待流水车间多目标调度问题,获得的非支配解集在多样性和性能上均有较好的表现。本研究仅求解了最大完工时间和总流程时间的双目标情形,在后续研究中,可考虑其他优化目标,如具有交货期的滞后目标,也可考虑将双目标推广到三目标的情形。