一种改进离散磷虾群的复杂产品装配调度算法

庄存波， 熊辉， 刘检华， 唐承统

(北京理工大学 机械与车辆学院， 北京 100081)

0 引言

复杂产品是指客户需求复杂、产品组成复杂、产品技术复杂、制造过程复杂、项目管理复杂的一类产品，如导弹、卫星、火箭、飞机等[1]。装配是复杂产品生产的最后环节，也是最为重要的环节之一，其直接关系到产品的质量、寿命、性能和可靠性[2]。生产调度是指在一定约束条件下，把有限资源在时间上分配给若干个任务，以满足或优化一个或多个性能指标的过程[3]。生产调度是产品装配过程的关键环节，也是装配过程管理与控制的核心问题之一。复杂产品装配多以流程为主组织生产，其装配过程是典型的工作流。在航天器等复杂产品的装配过程中，通常包括产品总装和组件部装两个阶段，每个阶段均需绘制总装工艺流程图和部装工艺流程图。总装装配和部装装配多数是分布在不同的装配车间，且二者的装配过程相似，故本文仅研究单个阶段的装配调度问题。复杂产品装配调度问题(总装或部装)可简化为一个以最大完工时间最小为调度目标、针对关键路径进行排产的多阶段混合流水车间调度问题(HFSP)。HFSP为流水车间调度问题(FSP)与并行机调度问题的集合，也被称为柔性FSP(FFSP)。与一般FSP相比， HFSP加大了机器的可选择性，扩大了可行解的搜索范围，是更加复杂的NP-hard问题。传统HFSP按并行机类型可分为相同并行机HFSP[4]、均匀并行机HFSP[5]及不相关并行机HFSP[6]共3种，按阶段数量可分为2阶段HFSP、3阶段HFSP及多阶段HFSP[7].

自1954年Johnson[8]研究具有两台机器、工件加工顺序已知的HFSP开始，国内外学者对HFSP进行了研究，并取得了显著成果。目前，HFSP求解算法主要包括精确算法、启发式算法和元启发式算法。精确算法中分支定界(B&B)算法使用最为广泛，如Brah等[9]利用B&B算法求解HFSP时提出了一种对工件排列和机器分配进行枚举的分支策略，Neron等[10]为了增强B&B算法搜索效率提出了一种基于推理的全局操作方法。精确算法在理论上能够获取问题的最优解，但面对大规模调度问题时，其求解时间过长，故只适用于求解小规模调度问题。Gupta[11]提出了一种简单的启发式算法来求解第1阶段含1台机器、第2阶段含2台相同机器的2阶段HFSP. Riane等[12]提出了2种新的简单启发式算法来求解3阶段相同并行机调度问题。对于多阶段HFSP，Ruiz等[13]和Pan等[14]在种群初始化阶段利用Nawaz-Enscore-Ham(NEH)启发式算法[15]生成少数优质个体以提高种群质量，从而优化算法寻优性能。启发式算法虽然能在较短时间内获得可行解，但难以保证解的质量。为了获得更优质的解，国内外许多学者在近些年更多采用元启发式算法来求解HFSP，如遗传算法(GA)[16-17]、禁忌搜索(TS)[18-19]、模拟退火(SA)[20]、人工免疫系统(AIS)[21-22]、迭代局部搜索(ILS)[23]、粒子群优化(PSO)[24-25]、人工蜂群(ABC)算法[26-27]和蚁群优化(ACO)[28-29]等群体智能优化算法(SIOA)，这些算法多用于求解多阶段相同并行机HFSP. 2012年，Gandomi等[30]通过模拟磷虾种群觅食行为提出了一种新颖的仿生SIOA，即磷虾群(KH)算法。在函数优化领域，KH算法及其改进算法在许多标准测试问题上被证明比一些熟知的元启发式算法更好或非常具有竞争性[31]。KH算法具有较强局部探索能力，且控制参数少，易于实现，非常适用于并行计算。然而，KH算法无法一直很好地实现全局搜索，易陷入局部最优，且KH算法更多依赖于随机运动，易导致收敛速度缓慢[32-33]。目前，KH算法已在特征识别、人工神经网络训练、比例—积分—微分(PID)控制、逆几何设计、基于地图的网络路径优化和结构优化等问题上得到广泛研究与应用[34]。如Kowalski等[35]将KH算法应用于训练人工神经网络，并通过仿真结果和算法比较得出KH算法在精度和时间上比其他元启发式算法和传统方法更有效。Yaghoobi等[36]提出了一种改进的混沌KH(ICKH)算法，将其应用于PID控制器参数的确定，通过对比发现ICKH算法比其他优化算法性能更优。KH算法在复杂产品装配车间调度问题和HFSP上的应用研究尚无文献可查。

本文以复杂产品装配车间为研究对象，采用基于排列的编码方法和基于启发式规则的解码方法，提出一种改进的离散KH(IDKH)算法，通过实验设计方法分析不同参数设置对算法性能的影响，最后通过实例和算法比较验证所提算法的高效性和鲁棒性。

1 复杂产品装配调度模型

1.1 问题描述

1.2 数学模型构建

基于上述分析，构建复杂产品装配调度问题的数学模型，具体描述如下：

minCmax，

(1)

(2)

k∈{1,2,…,s-1},

(3)

(4)

zghk+zhgk≤1,g,h∈J,g≠h,

(5)

(6)

xjkm,zghk∈{0,1},

(7)

式中：Cmax为最大完工时间；g、h为产品编号；m∈AT，AT为装配班组集合，AT={1,2,…,L}，L为装配班组总数；pjkm为产品j在工序k、由班组m进行装配操作所需要的时间；tjk为产品j的第k道工序开工时间；U为足够大的正数；xjkm=1为产品j的第k道工序被分派到装配班组m上，否则xjkm=0；zghk=1为在第k道工序上产品g先于产品h被装配，否则zghk=0. (1)式表示调度性能指标，即工期最小；(2)式表示最后一个产品在最后一道工序s的完工时间为Cmax；(3)式表示同一个产品只有在前一道工序装配完毕后才能开始下一道工序；(4)式表示每个产品在各个工序阶段只能由一个装配班组进行装配操作；(5)式表示产品g和产品h间存在先于、后于、同时3种顺序关系；(6)式表示同一个装配班组在一个时刻只能装配一个产品，即前一个产品装配完成后才能开始下一个产品的装配；(7)式表示变量的取值为0或1.

2 标准KH算法介绍

KH算法在优化过程中的目标是磷虾与食物之间距离最短及磷虾种群密度最大[30]。磷虾个体i(i∈{1,…,NP}，NP为种群数量)在t时刻的位置向量Xi由3个因素决定：1)受其他磷虾个体诱导所产生的运动向量Ni；2)觅食运动向量Fi；3)随机的物理扩散运动向量Di. 用拉格朗日模型来描述磷虾个体的移动位置向量：

(8)

1)Ni由局部种群密度(即局部影响)、目标种群密度(即目标影响)和排斥种群密度(即排斥影响)3个因素决定。磷虾个体i受其他磷虾个体诱导所产生的运动向量可表示为

(9)

2)Fi由食物位置和磷虾个体i的过往觅食经验2个因素决定。磷虾个体i的觅食运动向量可表示为

(10)

3)Di为一个随机过程，对于磷虾个体i，有

Di=Dmaxδ,

(11)

式中：Dmax为最大扩散速度，取值范围为[0.002,0.010]，一般取0.005；δ为随机方向向量，每项都是位于[-1,1]区间的随机数。

4)种群位置更新[30]：对于磷虾个体i，第gen+1次迭代后的位置向量可表示为

(12)

式中：Δt为速度向量的比例因子，是非常重要的常数之一，需根据优化问题仔细设定。Δt完全由搜索空间决定，可简单表示为

(13)

式中：c为一个常数，取值范围为[0，2];UBv和LBv为变量v的上界和下界。

5)遗传算子：受进化算法的启发，交叉算子和变异算子两种遗传繁殖机制被加入到KH算法当中，经过仿真实验发现，只将交叉算子应用到KH算法中综合性能最优。关于标准KH算法的详细介绍见参考文献[30].

3 IDKH设计

3.1 编码与解码

目前，求解HFSP的算法大多采用基于排列的方式对群体中个体进行编码，即取所有工件序号排列作为一个个体。群体中的每一个个体都对应一个可行调度解，个体长度为工件数量n，工件号在排列中位置代表其在第一个阶段加工顺序。工件排序和后续机器分配由解码规则决定[38-39]。

机器分配部分一般采用最先空闲机器规则[40]，但主要面向相同并行机调度问题。由于复杂产品装配调度问题属于不相关并行机调度问题，因此，本文基于该规则，提出如下改进策略：

1)对于每个工件j，选择完工时间最短的机器作为工件j的加工机器，完工时间为am+pjkm，其中，am为工件j最早允许加工时间且am=max (rm,Cj,k-1)，rm为机器m释放时间，Cj,k-1为工件j在上一阶段的完工时间；

2)对于完工时间相同的机器，基于效率最大化原则，优先选择加工时间pjkm最小的机器；

3)对于完工时间和加工时间都相同的机器，为了使调度结果更加紧凑，选择am-rm最小的机器。

工件排序部分采用如下改进策略：

1)第一阶段按照编码方式确定工件加工顺序，后续阶段基于先到先加工方式进行排序；

2)对于前一阶段完成时间相同的工件，选择当前阶段加工时间短的工件优先加工。对于复杂产品装配调度问题，可用产品序号代替工件序号，用装配代替加工，用装配班组代替机器。

3.2 种群初始化

为了保证种群的多样性以获得全局最优解，采用随机方式生成种群个体的位置向量X=(X1,X2,…,XNP)，个体i位置向量为Xi，变量数为工件数n，上下界范围可自定义。然而，标准KH算法只适用于连续优化问题的求解，并不适用于离散优化问题的求解，因此需将个体位置向量Xi转化为工件序号的排列πi后才能够被解码。

为此，本文提出转化策略，对于每一个磷虾个体i，其对应的位置向量为

Xi=(Xi1,Xi2,…,Xin)，

(14)

3.3 局部搜索

针对标准KH算法收敛速度较为缓慢的问题，本文采用局部搜索策略来加快算法收敛速度，并增强算法局部开采能力。目前，常用的局部搜索策略包括：1)交换策略：从排列编码中随机选择2个位置的元素进行交换；2)插入策略：从排列编码中随机选择2个位置，将第1个位置中的元素放到第2个位置上。

为了不影响算法运算速度，在每一次迭代过程中，只对目前出现的最优解执行局部搜索策略。本文采用一种修改的变邻域局部搜索策略，具体实现步骤如下：

1)对目前已知的最优解Xbest执行插入操作得到Xbest′.

2)对Xbest′执行混合邻域局部搜索策略(即在每一次局部搜索过程中，从交换策略和插入策略中随机选择一种，选择概率均为50%)得到其邻域X，如果邻域X适应度值比Xbest′适应度值低，则用X替代Xbest′进行下一次混合邻域搜索。该步骤迭代次数为n(n+1).

3)比较Xbest′与Xbest的适应度值，即比较f(Xbest′)与f(Xbest)，若f(Xbest′)

3.4 重启操作

针对KH算法全局搜索性能有限而导致的容易陷入局部最优问题，本文提出了一种重启策略。具体运行流程为：1)当种群最优解不更新代数Age达到设定常数Limit时，算法将执行重启操作；2)保留百分比为η的最优个体，剩下(1-η)·NP个个体则被随机生成的种群个体替代。重启操作伪代码如图2所示。

3.5 IDKH算法流程

基于上述分析，构建如图3所示的求解复杂产品装配调度问题的IDKH算法运算流程。

4 仿真结果与算法比较

4.1 参数讨论

采用MATLAB 2013平台对IDKH算法进行编程，并在性能为Inter(R) Core(TM) i7-3770 CPU 和3.40 GHz RAM的计算机上运行IDKH算法，操作系

SetLimit,η;

IfAge>=Limit

Keep the current bestXbest

Randomly generate the left (1-η)*NPindividuals

While not termination

Repeat

Evaluate the fitness of the krill

Perform the three motions

Implement the crossover operator

Update the krill position

一切都和过往的经验不一样了。在国家队时，邹习惯谨慎地谈论目标。第一场职业赛前新闻发布会，一个泰国记者看到他笑眯眯的，感到吃惊。而现在，他会训练自己不怒自威的样子，摆一张臭脸，不光给媒体看，更多给对手看。

Evaluate the krill based on its new position

Until a number ofNPreplications

Sort all the krill and find the current bestXbest′

Iff(Xbest′)

Xbest=Xbest′

END If

END while

END If

END

图2 重启操作的伪代码

Fig.2 Pseudo code of restart operation

统为Windows 7.

IDKH算法中，种群个数NP、控制时间间隔C、重启操作常数Limit和最优种群保留比例η等4个参数会对其性能产生影响。各参数水平取值如表1所示。

为了测试IDKH算法的性能，本文随机生成一系列计算实例，产品(工件)数量J∈{10,15,20}，

表1 IDKH算法参数水平表Tab.1 Combinations of parameter values of IDKH

工序(阶段)数S=5，装配(加工)工时为取值范围[3, 40]的离散随机整数分布，每阶段的并行机数均为3. 每种组合随机生成3个实例，共计NI=3×3=9个实例。为每个实例选择规模为L16(44)的正交实验，IDKH算法在每种参数组合下均独立运行20次，设定终止条件为最大迭代次数200，同时限制每次运行时间不超过20 s. 计算得到每个实例正交结果如表2所示。

表2 IDKH算法正交表Tab.2 Orthogonal array and response values

由表2可知，组合[3 2 4 1]的不同参数组合的平均值最小，是较为稳定的、获得最优参数值的组合，此时NP=80，C=1.0，Limit=100，η=0.1. 根据表2，种群个数NP、控制时间间隔C、重启操作常数Limit和最优种群保留比例η的极差和重要程度如表3所示。

表3 IDKH算法各参数极差表Tab.3 Statistics of response values

由表3可知，Limit和η极差最大，其次为NP，C最小，这表明重启操作对IDKH算法的性能影响最大。从数值上看，4个参数极差均不超过0.50%，说明IDKH算法鲁棒性较强。不同参数在不同水平下对算法性能的影响趋势图如图4所示。

由图4(a)可知，NP越大，算法性能越优，原因在于NP决定了群体在搜索空间的覆盖范围。当NP较小时，算法全局探索性能较差，易早熟收敛而陷入局部最优；当NP较大时，群体能够覆盖更多搜索空间，更易获得高质量的解。NP较大时，计算量大、耗时长，同时变优幅度随着NP增加不断变小(水平3和水平4差距不大，仅为0.03)，故NP不宜过大，NP=80为最佳。C决定了算法寻优过程中在整个搜索空间的移动幅度，C过大易导致收敛速度过快而陷入局部最优，C过小易导致收敛速度缓慢从而影响IDKH算法性能，由图4(b)可知，C=1.0时IDKH算法整体性能最优。对于重启操作，Limit决定了重启时机。当Limit较小时，IDKH算法还未达到局部最优就被迫重启，不仅重启效果一般，且会影响IDKH算法性能。η决定了重启时保留最佳个体比例，当保留比例较大时，IDKH算法易再次陷入局部最优，导致重启效果不明显。由图4(c)和图4(d)可知，当Limit=100，η=0.1时，重启效果最佳。综上所述，组合[3 2 4 1]为最佳参数组合，即NP=80，C=1.0，Limit=100，η=0.1.

4.2 初始实验

与离散KH(DKH)算法相比，IDKH算法主要是在局部搜索和重启操作两个方面进行了改进。

本文实验中，考虑对比IDKH算法及其3种变型：DKH算法、包含重启操作的DKH(DKH-RS)算法、包含局部搜索的DKH(DKH-LS)算法。为了验证DKH算法、DKH-RS算法、DKH-LS算法和IDKH算法求解HFSP的性能，采用4.1节中的9个实例对算法性能进行比较，选择最优参数组合NP=80，C=1.0，Limit=100，η=0.1，每个实例均运行20次，每次运行时间不超过20 s. 如图5所示为IDKH算法不同变型下的平均值、最优值和标准差变化趋势图。

通过DKH算法和DKH-LS算法对比以及DKH-RS算法和IDKH算法对比发现，加入局部搜索操作后，算法的平均值、最小值和标准差都会得到优化，说明局部搜索策略能够全面提高算法寻优性能和鲁棒性。通过DKH算法和DKH-RS算法对比以及DKH-LS算法和IDKH算法对比发现，加入重启操作后，算法的平均值和标准差得到了优化，即算法平均寻优性能和鲁棒性得到了增强 但最小值略有下降，原因在于算法还未达到局部最优就被迫重启了，然而大多数情况下并不会出现这种情况，如DKH-LS算法仅在一个实例上比IDKH算法获得的最小值小，其他均一致。综合来看，IDKH算法性能最优。

4.3 计算结果

4.3.1 文献实例验证

目前，求解HFSP比较常用的智能算法为GA[13,41]，近几年涌现出了一些新颖算法，如万有引力搜索算法(GSA)[42]、差分进化(DE)算法[43]、分布式估计算法(EDA)[38,44-45]。本文利用参考文献[38,46]所提3个HFSP实例L1、L2、L3，对比IDKH算法和DKH算法与GA、GSA、DE算法、EDA在求解实例时的性能差异。其中：DE算法数据来源于参考文献[43]，评价次数为10 000；GA、GSA、EDA数据均来源于参考文献[45]，采用运行时间10 s为终止条件。对于IDKH算法，同样设定终止条件为10 s，参数组合为NP=80，C=1.0，Limit=100，η=0.1，运算结果如表4所示。另外，由于智能算法均存在一些随机因子， 运行10次的结果有时并不能准确反映出算法性能(如算法稳定性)。因此，为了能够更加准确地验证所提算法性能，本文对DKH算法和IDKH算法运行70次，根据70次运行结果及其平均值得出10次运行结果。对于GA、GSA、DE算法、EDA(包括EDA_W算法[38]和EDA_I算法[45])，则采用参考文献[45]所提供的10次运算结果，得出其平均值。

表4 3个实例在不同智能算法中的结果对比

由表4可知， IDKH算法和EDA_I算法结果要明显好于其他算法，且都很稳定。二者求解性能相差不大，但IDKH算法略优于EDA_I算法，基本上每次都能获得目前已知的最优解，如：对于L2和L3，每次都能获得目前已知的最优解297和13；对于L2，在运行70次过程中，仅有一次没有获得已知的最优解21，运行70次得到的平均值为21.014 3. 其次为GSA和DKH算法，再次为GA和EDA_W算法，最后为DE算法，从而可知IDKH算法在求解HFSP时是高效且稳定的。

4.3.2 工程实例验证

以某航天器装配为例，在调度开始时刻，需要完成10个产品的装配，其装配工艺流程如图6所示。装配总工序数为7，关键路径上的工序数为4，分别标记为S1、S2、S3、S4. 关键路径上的每道工序均有3个装配班组,分别标记为M1～M12，及5个装配工位(满足装配工位数量大于等于装配班组数量的前提)可供分配，不同装配班组完成不同装配工序所需要的装配时间如表5所示。

基于该实例，分别用IDKH算法、DKH算法、EDA_W算法进行求解，每种算法均运行70次，每次最大运行时间仍为10 s.求解发现只有IDKH算法每次都能获得目前已知的最优解Cmax=195，其调度甘特图如图7所示。图7中，矩形框的数字表示产品序号j.

表5 某航天器的装配时间表Tab.5 Assembly time of a spacecraft h

5 结论

针对复杂产品装配调度问题，引入了一种高效的群体智能优化算法—KH算法，并采用局部搜索和重启策略增强了KH算法的局部寻优能力和全局搜索能力。针对已有参考文献[38,46]中3个不同实例和1个工程实例，IDKH算法求解结果明显比GSA、DKH算法、GA、EDA_W算法、DE算法等算法好，略优于EDA_I算法，且基本上每次都能获得目前已知的最优解，反映出IDKH算法在求解HFSP和复杂产品装配调度问题时是高效且稳定的。