基于物料延期干扰管理的飞机装配线动态调度

摘要

为了提升飞机装配过程中应对物料供应延期干扰的能力，对物料供应延期情况下飞机装配线动态调度问题进行研究。通过对物料交付确定性和预测性两类信息的分析将作业分为四类进行调度决策，建立了基于时间和事件混合驱动的周期性滚动调度框架和对应的数学优化模型。针对该问题设计了一种改进离散粒子群算法，将Levy飞行和离散化解码机制嵌入算法，滚动求解调度模型。数值实验结果表明，本文提出的调度方法能充分利用物料供应信息，在确定本阶段决策的同时前瞻性地进行预测决策，相比于其他动态调度方法，该方法能更有效地应对各类物料供应干扰，满足飞机装配动态调度的需求。

关键词

飞机装配线; 物料供应延期; 干扰管理; 周期性滚动调度

飞机装配是飞机制造的重要环节，在现代飞机制造中，装配协调的工作量占到飞机制造总工作量的50%~60%，因此提高装配可靠性和效率具有重要意义^［1-2］。随着我国飞机总装技术的发展，装配方式由传统机库式转向脉动式、移动式，能有效缩短总装作业周期。但飞机总装厂还需考虑外部零件供应商的各种约束和影响。飞机装配零部件种类繁多^［3］，供应商的交付易受突发事件影响，导致物料供应不准时，增加了装配项目进度延误的风险。因此，研究一种能够有效应对物料供应延期干扰的调度方法具有重要意义。

物料供应延期干扰下的飞机装配线调度问题本质上是一类不确定环境下的项目调度问题，目前大多数关于飞机装配线调度相关问题的研究集中在模板装配计划的制定，强调装配工艺的先后顺序约束^［4］，或关注于装配作业执行中可能发生异常情况干扰。Zaman等^［5］针对项目工期不确定性，提出融合双搜索算子与双启发式算法的进化框架，启发式算法优化作业紧凑性；为增强鲁棒性，随机生成多场景进行仿真评价。Chand等^［6］根据针对项目资源中断问题，设计了基于遗传算法的超启发式算法GPHH，从历史数据学习启发式规则求解动态调度问题。Chakrabortty等^［7］研究了单中断和多中断的反应式重调度，通过资源分解重新安排作业，提出混合整数线性规划模型以最小化完工时间。但上述研究仅限于反应式调度，抗风险能力弱。由于大量不确定风险是可预知的，Ma等^［8］从前摄调度方法入手，针对资源和工期两种不确定环境，提出基于松弛的鲁棒性度量框架，分析资源与时间缓冲区作用。但前摄调度是一次性的，仍需通过动态调整才能适应环境变化，因此部分学者结合前摄与反应式调度，以增强初始计划鲁棒性并确保调度系统平稳运行。Davari等^［9］综合考虑前摄调度和反应调度所产生的成本，提出了结合最小化模板计划成本和反应调度成本加权和的目标函数，同时给出基于动态规划的模型来解决该问题。Brčić等^［10］假设已知作业不确定工期的分布，基于前摄-反应式调度框架，设计了基于部署和迭代策略搜索的元启发式方法进行优化求解。上述研究聚焦作业时间不确定性，忽略了物料供应不确定性和项目执行中的不确定性变化。针对物料供给延期干扰，陆志强等^［11］以最小化调度计划变动成本和资源额外投入成本加权和为目标，设计了一种将支持向量数据描述技术与传统的调度方法相结合的算法，实现了动态调度的快速响应。该研究虽考虑物料供应不确定性，但仅预估延期时间。若利用供应商历史交付信息预测延期，并按项目进程迭代调整，可进一步增强调度计划的鲁棒性。

近年来，物联网、云计算等技术使实时信息成为生产调度优化的关键。大型客机主制造商与供应商建立数据共享的合作伙伴关系^［12］，这使得主制造商能够构建一种综合供应商和物流信息的供应商交付预测系统，该系统能够提供交付预测概率，实时监测物料供应，使管理人员能提前灵活应对异常。随着交期临近，预测结果持续更新，生成的调度计划更具鲁棒性。本文针对供应延期干扰下的飞机装配线调度问题，借鉴文献［13-14］中的两阶段周期性滚动决策框架，建立供应延期干扰下的周期性动态调度框架，并设计改进离散粒子群算法求解。该算法既能利用确定性信息调整计划，又能预见性安排后续决策，同时保持计划调整的灵活性。

1 问题描述与数学模型

1.1　问题描述与基本假设

物料供应延期干扰下的飞机装配线调度问题，是为了在预防物料延期导致项目停滞，充分利用实时供应商交付信息，通过周期性滚动调度重新决策未开始作业的时间，并结合预测模型前瞻调度后续作业，增强装配计划的鲁棒性。基于飞机装配特点，假设如下：①调度问题视为资源受限项目调度，需满足时序和资源等基本约束；②工位共享多种可更新资源，如工人、设备、能源和存储空间；③模板计划预先确定；④物料动态到达且齐套配送；⑤时间轴离散化。

表1汇总了本文所用到的参数及决策变量，由于每个装配作业都有自己的紧前和紧后作业，所以在整个项目的首尾各增加一个虚作业，作业编号为0和$n+\mathrm{1}$，其执行工期与资源需求量为0。

表1  参数和决策变量

Tab.1  Parameters and decision variables

参数	定义
$J$	作业集合
$n$	作业数量
$j\in \left\{\mathrm{0,1},\dots ,n,n+1\right\}$	作业序号
$K$	可更新资源集合
$k\in K$	资源种类编号
$D$	时间集合
$t\in D$	离散时间节点编号
$d_j$	作业$j$的执行工期
$P_j$	作业$j$的紧前作业集合
$r_{jk}$	作业$j$对第$k$种资源的需求量
$R_k$	第$k$种资源总量
$s_j^0$	作业$j$的模板计划开始时间
$L_j^0$	作业$j$的物料计划到达时间
${L_{\mathrm{T}}}_j$	作业$j$的物料配送提前期
$A_j^{\theta }$	场景$\theta$下作业$j$的物料实际到达时间
$q\in Q$	滚动周期集合
$T_q$	第$q$个调度周期开始时刻
$s_j^R,j\in J_{\mathrm{F}}\bigcup J_{\mathrm{A}}$	$J_{\mathrm{F}}$和$J_{\mathrm{A}}$类作业$j$的实际开始时间
决策变量$x_{jt}^{\theta }$	0/1变量，场景$\theta$下作业$j$在时刻$t$开始取1，否则取0
辅助变量$s_j^{\theta }$	整数变量，场景$\theta$下作业$j$的开始时间

1.2　供应延期干扰下的滚动调度框架

基于前文所述，当某一作业的物料供应延期时间已知时，可确定该作业的最早开始时间并调度；其余可能发生延期的作业物料到达时间由供应商预测系统提供预测值，据此进行预测性调度。尽管这些干扰并未确定，但隐含在后续计划执行过程中，预先调度能降低重调度成本，进而降低项目总成本。为了降低问题复杂度以便于求解，本文采用滚动调度的方法，将动态调度问题按滚动机制分解成一系列多阶段的静态子问题；并选择混合机制，发生定义的突发事件时采用事件驱动，否则采用周期调度，以固定时间周期作为滚动窗口推动更迭。根据上述分析，对于任意重调度周期，给如下四个定义。

定义1：项目中已完成的所有作业构成的集合称为该时刻的已完成作业集，记为$J_{\mathrm{F}}$。

定义2：项目中已经开始但未完成的作业构成的集合称为该周期的已调度待执行作业集，记为$J_{\mathrm{A}}$。

定义3：物料送达时刻位于重调度周期内的作业集合称为周期内待调度作业集，记为$J_{\mathrm{D}}$。

定义4：物料送达时刻位于重调度周期后的作业集合称为周期外作业集，记为$J_{\mathrm{P}}$。

如图1所示，将第$q$个调度周期开始时刻记为$T_q$，$T_q$前开始执行的作业（即$J_{\mathrm{F}}$和$J_{\mathrm{A}}$类作业）无需安排；$T_q$到$T_{q+\mathrm{1}}$之间的作业（即$J_{\mathrm{D}}$类作业）为确定性决策；$T_{q+\mathrm{1}}$后的作业（即$J_{\mathrm{P}}$类作业）为预测性决策，决策模型为不确定模型，需通过构建场景的方法建模求解，决策结果用于评价确定性决策的优劣。

图1  ${\mathit{T}}_{\mathit{q}}$时刻的作业调度计划

Fig. 1  Job scheduling plan at ${\mathit{T}}_{\mathit{q}}$

滚动周期由企业根据实际情况调整，本文假定已设为合适值。若交付预测未更新，则当前周期调度沿用上一周期计划。在周期内物料到达信息在大多数情况下为确定性信息，但在少数情况下会出现突发事件。将由于工厂内部物流异常等情况导致的$J_D$类物料供应延期情况定义为突发事件。

由于突发事件的延期时间短、概率低，无需为此专门设计调度算法，若不影响计划则忽略，否则按规则调整，干扰事件发生的时间点为调度决策点。图2展示了整个滚动调度框架的结构。

图2  滚动调度基本框架

Fig. 2  Basic framework for rolling scheduling

1.3　数学模型

根据上述调度逻辑，模板装配计划的制定耦合了装配资源配置计划，供应延期时模板计划失效，重调度将产生额外成本并延长工期。因此目标函数需综合考虑实际执行计划与模板计划的偏差和装配项目工期两类指标，采用加权求和的方法转化为单目标函数。记一次重调度决策为$x$，一次重调度的目标函数公式如式（1）所示，其中${\omega }_{\mathrm{1}}$和${\omega }_{\mathrm{2}}$表示指标的权重，满足${\omega }_{\mathrm{1}}+{\omega }_{\mathrm{2}}=\mathrm{1}$。

$$f(x)={\omega }_{\mathrm{1}}{\displaystyle \sum _{j=\mathrm{0}}^{n+\mathrm{1}}}\left|s_j-s_j^{\mathrm{0}}\right|+{\omega }_{\mathrm{2}}{s}_{n+\mathrm{1}}$$

（1）

根据上述前提条件，参考文献［14］的建模方法，对于作业所需物料到达时间不确定的情况采用蒙特卡洛抽样的方法，从候选值中随机抽样，从而获得仿真样本。

含有$n$个作业物料到达时间的序列表示为$\theta =\left(A_{\mathrm{1}},A_{\mathrm{2}},\cdots ,A_n\right)$，将每次抽样后获得的序列$\theta$定义为一个场景，即一种可能的物料到达情况。重复抽样$M$次后，物料到达场景集合表示为$\Theta =\left({\theta }_{\mathrm{1}},{\theta }_{\mathrm{2}},\cdots ,{\theta }_M\right)$。通过物料到达场景集合将不确定性问题转化为若干个场景下的优化问题，其中每个场景下的优化问题均包含对$J_{\mathrm{D}}$类作业的固定决策（记为$S_{\mathrm{1}}$）和$J_P$类作业的预测决策（记为$S_{\mathrm{2}}^{\theta }$，$\theta \in \Theta$），分别对应固定成本$f\left(S_{\mathrm{1}}\right)$和期望成本$E[f(S_{\mathrm{2}}^{\theta })]$，其中$f\left(S_{\mathrm{1}}\right)={\displaystyle \sum _{j=\mathrm{0}}^{n+\mathrm{1}}}\left|s_j-s_j^{\mathrm{0}}\right|$，$f\left(S_{\mathrm{2}}^{\theta }\right)=s_{n+\mathrm{1}}^{\theta }$。目标函数可表示为两者加权和，通过不同场景下的预测决策，评估确定性决策的优劣以及应对不同场景的反应能力，达到兼顾决策反应性和预测性的目的。模型的决策变量为项目中各作业的开始时间，决策变量统一用场景进行表述。数学模型表述如下：

$$\mathrm{m}\mathrm{i}\mathrm{n}Z={\displaystyle \sum _{\theta }^M}\left({\omega }_{\mathrm{1}}{\displaystyle \sum _{j=\mathrm{0}}^{n+\mathrm{1}}}\left|s_j^{\theta }-s_j^{\mathrm{0}}\right|+{\omega }_{\mathrm{2}}{s}_{n+\mathrm{1}}^{\theta }\right)/M$$

（2）

$$\mathrm{s}.\mathrm{t}.$$

$$s_j^{\theta }=s_j^R,\forall j\in J_F\bigcup J_A,\theta \in \Theta$$

（3）

$$s_j^{{\theta }_{\mathrm{1}}}=s_j^{{\theta }_{\mathrm{2}}},\forall j\in J_D,{\theta }_{\mathrm{1}},{\theta }_{\mathrm{2}}\in \Theta$$

（4）

$$s_j^{\theta }={\displaystyle \sum _{t=\mathrm{0}}^D}\left(x_{jt}^{\theta }t\right),\forall j\in J,\theta \in \Theta$$

（5）

$$s_j^{\theta }\ge s_p^{\theta }+d_p,\forall j\in J,p\in P_j,\theta \in \Theta$$

（6）

$${\displaystyle \sum _{j=\mathrm{1}}^n}\left(r_{jk}{\displaystyle \sum _{\tau =t-p_j+\mathrm{1}}^t}{x}_{j\tau }^{\theta }\right)\le R_k,\forall k\in K,t\in D,\theta \in \Theta$$

（7）

$${\displaystyle \sum _{t=\mathrm{0}}^D}{x}_{jt}^{\theta }=\mathrm{1},\forall j\in J,\theta \in \Theta$$

（8）

$${\displaystyle \sum _{t=\mathrm{0}}^D}\left(x_{jt}^{\theta }t\right)-L{T}_j\ge A_j^{\theta },\forall j\in J,\theta \in \Theta$$

（9）

$${\displaystyle \sum _{t=\mathrm{0}}^D}\left(x_{jt}^{\theta }t\right)-L{T}_j\ge T_q,\forall j\in J_D\bigcup J_P,\theta \in \Theta ,q\in Q$$

（10）

$$x_{jt}^{\theta }=\left\{\mathrm{0,1}\right\},\forall j\in J,t\in D,\theta \in \Theta$$

（11）

其中，目标函数式（2），由式（1）改写而来，受式（3）~（11）的约束。式（3）表示任意场景下任意$J_A$和$J_F$类作业的决策保持一致；式（4）表示对$J_D$类作业的决策在任意场景下保持一致；式（5）表示决策变量$x_{jt}^{\theta }$与$s_j^{\theta }$之间的关系；式（6）表示任意场景下需要满足作业间的时序约束；式（7）表示任意场景下需要满足作业执行时的资源约束；式（8）表示任意场景下任意作业一旦开始不能中断；式（9）表示任意场景下任意作业的物料开始配送时刻不应早于物料到达时刻；式（10）表示任意场景下任意$J_{\mathrm{D}}$和$J_{\mathrm{P}}$类作业的物料开始配送时刻不应早于干扰发生时刻；式（11）表示决策变量$x_{jt}^{\theta }$的可行域。

2 算法设计

基于上述模型，本文设计了以一种改进离散粒子群算法（improved discrete particle swarm optimization，IDPSO）。该算法调整并优化了基本粒子群算法的解码部分，以适应问题约束，并引入Levy飞行作为粒子更新机制，通过模拟鸟类和鱼类的运动模式避免局部最优，提高搜索效率。算法设计思路如下：基于粒子群算法的一般步骤，IDPSO算法搜索$J_{\mathrm{D}}$类作业的确定性决策和不同场景下$J_{\mathrm{P}}$类作业的预测性决策，通过计算目标函数值评估确定性决策结果的优劣。IDPSO算法框架及流程如图3所示，具体步骤如下：

图3  IDPSO算法流程图

Fig. 3  The flow chart of IDPSO

步骤1：初始化各项参数，设置迭代次数$i=\mathrm{0}$，初始化粒子位置；

步骤2：根据2.1节所述的解码方案对每个粒子进行解码；

步骤3：计算每个粒子的适应度值，确定各粒子的最优位置$x_{id}^{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$和全局最优位置$x_d^{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$；

步骤4：判断是否满足算法终止条件$i<\mathrm{I}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{M}\mathrm{a}\mathrm{x}$，若不满足则转步骤5；否则转步骤6；

步骤5：根据2.2节所述粒子更新机制更新粒子的速度和位置， $i=i+\mathrm{1}$，转步骤2；

步骤6：算法结束，输出结果。

2.1　编码与解码设计

本文采用基于优先权列表的编码方案，其中粒子位置值代表作业优先权，位置序号即作业编号。针对传统粒子群算法不适用于离散优化及存在解码冗余的问题，本方案通过优先权排序消除冗余，并通过共享最优解特征引导种群向全局最优进化。本文的解码方案步骤如下，需要注意的是，当算法第一次迭代时，由于此时还没有决策出当前最优解，在第一次解码时需要跳过步骤3。

步骤1：将向量$X=\left\{x_{\mathrm{0}},x_{\mathrm{1}},\cdots ,x_{n+\mathrm{1}}\right\}$中的连续值从小到大进行排序，得到排序结果向量$L=\left\{l_{\mathrm{0}},l_{\mathrm{1}},\cdots ,l_{n+\mathrm{1}}\right\}$。

步骤2：将向量中每个连续值替换为其排名的序号$X=\left\{l_{\mathrm{0}},l_{\mathrm{1}},\cdots ,l_{n+\mathrm{1}}\right\}$。

步骤3：以概率$Q$随机选择一些位置$p_j$，根据当前最优解相应位置$X^{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$上的值进行替换，得到新的作业优先值编码$X^{\text{'}}=\left\{l_{\mathrm{0}},l_{\mathrm{1}},\cdots ,l_j^{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}},\cdots ,l_{n+\mathrm{1}}\right\}$。

步骤4：根据串行进度生成机制对优先值列表进行解码，比较被选作业优先值和它的紧前任务优先值的大小，若小于紧前任务，则无须调整； 若大于紧前任务，则交换两者的优先值。

步骤5：若出现重复的解码结果，则对重复的结果进行一次2-opt邻域交换操作。

步骤6：当所有粒子解码完毕且没有重复解时，算法结束，输出结果。

2.2　基于Levy飞行的粒子更新机制

粒子群算法（PSO）是一种基于群体智能的优化方法，在PSO中根据式（12）更新速度，根据式（13）更新位置，当粒子运动速度超过规定值$V_{\mathrm{m}\mathrm{a}\mathrm{x}}$时，使用式（14）重置速度。其中$x_{id}^k$，$v_{id}^k$表示在$k$时刻粒子$i$在$d$维的位置和速度，$x_{id}^{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$表示粒子$i$经过的历史最佳位置，$x_d^{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$表示群体所经过的历史最佳位置，$\omega$表示粒子的惯性权重，$c_{\mathrm{1}}$和$c_{\mathrm{2}}$称为学习因子，$\xi$和$\eta$是0―1之间的随机数。为克服PSO在离散优化中易产生冗余解、多样性不足的缺陷，采用Levy飞行策略替代均匀分布，防止粒子陷入局部最优。

$$v_{id}^{k+\mathrm{1}}=\omega v_{id}^k+c_{\mathrm{1}}\xi \left(x_{id}^{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}-x_{id}^k\right)+c_{\mathrm{2}}\eta \left(x_d^{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}-x_{id}^k\right)$$

（12）

$$x_{id}^{k+\mathrm{1}}=x_{id}^k+v_{id}^{k+\mathrm{1}}$$

（13）

$$v_{id}^{k+\mathrm{1}}=V_{\mathrm{m}\mathrm{a}\mathrm{x}}$$

（14）

在Levy飞行的改进方面参考文献［15-16］，根据式（15）得到飞行的步长，其中参数$u$，$v$由式（16）~（18）决定，根据文献［16］设置${\sigma }_v=\mathrm{1}$，$\lambda =\mathrm{1.5}$，将式（12）改写为式（19）。

$$L\left(\lambda \right)=\frac{u}{{\left|v\right|}^{\frac{\mathrm{1}}{\lambda }}}$$

（15）

$$u\sim N\left(\mathrm{0},{\sigma }_u^{\mathrm{2}}\right)$$

（16）

$$v\sim N\left(\mathrm{0},{\sigma }_v^{\mathrm{2}}\right)$$

（17）

$${\sigma }_u={\left[\frac{\Gamma \left(\mathrm{1}+\lambda \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\pi \lambda }{\mathrm{2}}\right)}{\Gamma \left(\frac{\mathrm{1}+\lambda }{\mathrm{2}}\right)\lambda {\mathrm{2}}^{\frac{\lambda -\mathrm{1}}{\mathrm{2}}}}\right]}^{\frac{\mathrm{1}}{\lambda }}$$

（18）

$$v_{id}^{k+\mathrm{1}}=\omega v_{id}^k+c_{\mathrm{1}}L\left(\lambda \right)\left(x_{id}^{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}-x_{id}^k\right)+c_{\mathrm{2}}L\left(\lambda \right)\cdot \left(x_d^{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}-x_{id}^k\right)$$

（19）

2.3　参数设置

为了确定算法所用参数设置，以120规模的5个算例作为对象进行预实验。由于算例规模越大，计算复杂度越高，设置种群规模$P_{\mathrm{S}}=\mathrm{0.5}n$；根据文献［17］的研究结论设置$\omega =\mathrm{1}$，$c_{\mathrm{1}}=\mathrm{2}$，$c_{\mathrm{1}}=\mathrm{1}$，$V_{\mathrm{m}\mathrm{a}\mathrm{x}}$为历史最大目标值$Z_{\mathrm{m}\mathrm{a}\mathrm{x}}$，$Q=\mathrm{0.1}$。根据IDPSO算法预实验的收敛曲线设置迭代次数$\mathrm{I}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{M}\mathrm{a}\mathrm{x}=\mathrm{50}$。

3 数值实验

以文献［11］中模板计划实例为参考，结合本问题的特点划分作业数量为20、30、60（小规模）和90、120（大规模）的5种类型，分别生成10组算例，初始目标权重设置为${\omega }_{\mathrm{1}}=\mathrm{0.5}$，${\omega }_{\mathrm{2}}=\mathrm{0.5}$，采用最迟完工时间（LFT）优先规则作为初始模板计划的生成方法。假设受供应延期干扰的作业占项目中作业总数的10%，其余作业的物料到达时间均已知，提前期在［0，20］内随机选取；受干扰作业的物料实际延期时间在［5，15］内随机选取。突发事件的发生概率为0.05，延期时间在［1，5］内随机选取，当突发事件发生时采用右移规则（right shift，RS）进行调整。作业$j$预估的物料到达时间与实际的偏差$L{T}_j^{\text{'}}$根据时间的推进不断更新，如式（20）所示，其中$\gamma$表示衡量预测算法准确性的系数，本文中设为5。

$$L{T}_j^{\text{'}}=N\left(\mathrm{0},\left(L_j^{\mathrm{0}}-A_j^{\theta }\right)/\gamma \right)$$

（20）

滚动周期在小规模算例中设置为7，大规模算例中设置为14。在每个决策点设置场景池，大小为2 000，实验时从中抽取100个场景作为物料到达场景集合。计算后验精确解使用Gurobi软件。

3.1　算法调度求解有效性验证

为了验证IDPSO算法的有效性，本节设计如下对比实验，利用后验信息（即假设所有物料延期时间已知）进行反应式调度，单次求解数学模型。由于Gurobi软件不适合求解大规模算例，仅对小规模中各5组不同算例进行求解，其中$Z$为目标值，$\mathrm{G}\mathrm{A}\mathrm{P}$为其他算法与Gurobi求解结果的差值百分比，结果如表2所示。

表2  算法求解性能实验结果

Tab.2  Experiment result of algorithm performance

算例		实验1		实验2			实验3					实验4
		$Z_{\mathrm{G}}$		$Z_{\mathrm{I}\mathrm{D}\mathrm{P}\mathrm{S}\mathrm{O}}$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$		$Z_{\mathrm{D}\mathrm{E}}$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$	$Z_{\mathrm{T}\mathrm{S}}$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$		$Z_{\mathrm{R}\mathrm{S}}$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$
J20	1	63.15		63.68	0.84		63.67	0.82	63.88	1.16		65.88	4.33
	2	65.37		65.37	0		65.37	0	65.37	0		66.72	2.07
	3	69.52		69.57	0.07		69.60	0.12	69.90	0.55		71.72	3.16
	4	63.32		63.32	0		63.32	0	63.33	0.03		64.18	1.37
	5	66.83		66.83	0		66.95	0.17	67.90	1.60		69.72	4.31
	均值	65.64		65.75	0.18		65.78	0.22	66.08	0.67		67.64	3.05
J30	1	147.00		147.62	0.42		147.97	0.66	152.92	4.02		154.60	5.17
	2	153.75		154.91	0.75		155.03	0.83	156.68	1.91		163.57	6.38
	3	166.13		166.33	0.12		167.00	0.52	167.43	0.78		174.47	5.02
	4	159.28		160.13	0.53		161.82	1.59	162.13	1.79		167.67	5.26
	5	153.37		156.42	1.99		157.43	2.65	159.33	3.89		163.85	6.84
	均值	155.91		157.08	0.76		157.85	1.25	159.70	2.48		166.43	5.73
J60	1	361.23		365.67	1.23		366.18	1.37	368.33	1.97		373.05	3.27
	2	384.18		388.33	1.08		395.35	2.91	396.92	3.31		412.05	7.25
	3	371.57		376.13	1.23		375.97	1.18	375.33	1.01		406.85	9.50
	4	378.40		382.50	1.08		383.32	1.30	382.50	1.08		393.67	4.03
	5	375.10		379.88	1.27		385.05	2.65	385.68	2.82		414.45	10.49
	均值	374.10		378.50	1.18		381.17	1.88	381.75	2.04		400.01	6.91

实验1：使用Gurobi软件单次求解数学模型。

实验2：使用IDPSO算法单次求解数学模型。

实验3：与其他元启发式算法的对比：使用文献［18］提出的差分进化算法（DE）和文献［19］提出的禁忌搜索算法（TS）作为对比算法，单次求解数学模型。

实验4：采用传统的右移策略单次求解数学模型。

与实验1相对比，IDPSO的GAP值低（0.18%~1.18%），部分算例得精确解，求解质量好。与两种元启发式算法DE和TS相比，DE的求解结果较好，在各规模下的GAP值较低，但是计算速度较慢；TS计算速度较快，但因算法易陷入局部最优导致求解质量一般。而IDPSO算法嵌入Levy飞行机制，能有效扩大搜索范围，跳出局部最优。与传统的RS算法相比，尽管RS的求解速度最快，但求解质量较差。在计算时间方面，IDPSO算法在各规模算例中均表现优异且计算时间较短。Gurobi求解器虽在J20规模高效，但随规模增大，计算时间显著增加。J20规模下，四种算法用时约0.05min；J60规模时，IDPSO、DE、TS分别用时0.25min、0.40min、0.12min，Gurobi增至6.78min，是J20规模时的170倍。

3.2　物料供应干扰下的对比实验

为了验证本文滚动调度框架的有效性，本节设计了在物料供应延期干扰环境下的如下对比实验。将基于后验信息使用Gurobi（小规模）或IDPSO算法（大规模）进行求解的结果作为对比基准。

实验5：采用本文提出的滚动调度框架和IDPSO算法滚动求解数学模型。

实验6：根据本文提出的滚动调度框架采用右移策略滚动求解数学模型。

实验7：根据完全反应式调度策略采用IDPSO算法滚动求解数学模型。

实验8：根据预测―反应式调度策略（即在制定模板计划时考虑物料配送的不确定性，根据预测结果进行预测性调度，在项目执行过程中采用反应式调度）采用IDPSO算法滚动求解数学模型。

表3和表4结果显示，本文方法在各作业规模下的

表3  小规模算例实验结果

Tab.3  Experiment result of small-scale example

算例		后验结果		实验5				实验6		实验7			实验8
		$Z$		$Z_5$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$		$Z_6$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$		$Z_7$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$		$Z_8$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$
J20	1	54.12		55.48	2.51		63.19	16.760		64.14	18.51		59.58	10.09
	2	73.16		75.73	3.51		82.16	12.30		85.41	16.74		77.94	6.53
	3	68.64		70.61	2.87		74.48	8.51		80.44	17.19		73.58	7.20
	4	63.36		64.73	2.16		69.60	9.85		72.28	14.08		67.65	6.77
	5	64.97		64.98	0.02		85.85	32.14		80.40	23.75		74.25	14.28
	均值	64.85		66.31	2.21		75.06	15.91		76.53	18.06		70.60	8.97
J30	1	124.71		124.74	0.02		137.14	9.96		143.93	15.41		124.74	0.02
	2	116.91		120.05	2.68		125.15	7.05		134.23	14.81		126.16	7.91
	3	175.73		181.14	3.08		185.33	5.46		185.7	5.68		188.46	7.25
	4	203.84		206.69	1.40		223.01	9.41		225.70	10.73		218.66	7.27
	5	181.05		187.90	3.78		201.19	11.12		221.70	22.45		195.04	7.73
	均值	160.45		164.10	2.19		174.36	8.60		182.25	13.81		170.61	6.04
J60	1	462.22		467.70	1.19		525.90	13.78		475.60	2.89		479.98	3.84
	2	344.38		360.68	4.73		391.73	13.75		409.97	19.04		369.05	7.16
	3	345.00		353.82	2.56		385.12	11.63		379.58	10.02		371.48	7.68
	4	406.52		417.50	2.70		445.83	9.67		437.60	7.65		433.47	6.63
	5	386.70		387.43	0.19		421.33	8.96		433.70	12.15		395.15	2.19
	均值	388.96		397.43	2.27		433.98	11.56		427.29	10.35		409.83	5.50

表4  大规模算例实验结果

Tab.4  Experiment result of large-scale example

算例		后验结果		实验5			实验6			实验7			实验8
		$Z$		$Z_5$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$		$Z_6$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$		$Z_7$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$		$Z_8$	$\mathrm{G}\mathrm{A}\mathrm{P}/\%$
J90	1	613.23		625.73	2.04		654.60	6.75		656.43	7.04		644.98	5.18
	2	393.08		411.03	4.57		459.88	16.99		471.40	19.93		425.85	8.34
	3	410.45		416.30	1.43		457.90	11.56		485.50	18.28		428.25	4.34
	4	474.73		495.50	4.37		567.63	19.57		554.67	16.84		494.03	4.07
	5	488.47		504.03	3.19		563.50	15.36		537.17	9.97		530.60	8.63
	均值	475.99		490.52	3.12		540.70	14.05		541.03	14.41		504.74	6.11
J120	1	754.55		779.80	3.35		815.10	8.02		837.25	10.96		805.55	6.76
	2	523.10		528.60	1.05		577.30	10.36		598.60	14.43		554.45	5.99
	3	640.55		666.90	4.11		707.10	10.39		723.10	12.89		687.45	7.32
	4	807.40		837.75	3.76		905.95	12.21		909.30	12.62		872.35	8.04
	5	662.65		690.35	4.18		783.85	18.29		797.65	20.37		705.45	6.46
	均值	677.65		700.68	3.29		757.86	11.85		773.18	14.25		725.05	6.92

GAP值为2.21%~3.29%，部分算例结果与后验结果接近，证明算法有效利用预测信息进行前瞻性调度，求解稳定。实验7的GAP值为10.35%~18.06%，求解效果一般且稳定性差，因其为完全反应式调度，不考虑后续延期，导致额外的偏差成本。实验8的GAP值为5.50%~8.97%，求解结果和稳定性优于实验7但逊于实验5，因为前摄―反应式调度框架虽在制定计划时考虑延期信息，但在项目执行中仍需根据实际情况调整。

图4展示了J30规模下某算例在各重调度阶段时的偏差项值，本文算法框架能及早根据预测值吸收后续的不确定性，随着项目推进偏差项值显著减少。

图4  各实验重调度阶段的目标函数偏差项值

Fig. 4  Deviation value of objective function in each experimental rescheduling stage

3.3　不同干扰强度下的性能对比

供应延期干扰强度对算法性能和结果有很大影响：干扰极小时，问题可近似为一般RCPSP；干扰极大时，需更好地平衡质量鲁棒性和解鲁棒性。

本节以J30和J90规模的5组算例，通过改变干扰作业数量（占5%、10%、15%）验证方法有效性，采用实验5~8的框架和算法，结果如图5所示。干扰作业少且规模小时，改进不明显；随着干扰作业增多或规模扩大，本文方法优势显著。

图5  不同干扰强度下的结果对比

Fig. 5  Comparison of results under different interference intensities

3.4　不同目标函数权值对比

通过调整目标函数中权重系数${\omega }_{\mathrm{1}}$和${\omega }_{\mathrm{2}}$的取值，模拟不同企业需求下算法性能的优劣。以J30和J90规模的5组算例为测试对象，采用实验5~8的框架和算法，记录不同权重下的$Z$值。

实验结果如图6所示，随着偏差项权重增加，各方法$Z$值上升，可见本问题中偏差项影响更大；观察偏差项和工期项的差异关系发现，不同权重下各方法求解的工期值差异较小；而对于偏差值，不同方法在不同权重下结果差异明显，采用前摄方法的结果更优说明了利用预测信息的调度计划鲁棒性更强。本文方法在各权重下均优于其他三种，能兼顾质量鲁棒性和解鲁棒性，适用于各类需求的飞机装配线动态调度。

图6  不同目标函数权值结果对比

Fig. 6  Comparison of different weights for objective functions

4 结语

本文研究了基于物料供应延期干扰的飞机装配线作业调度问题，建立了基于时间和事件混合驱动的滚动调度框架，利用确定性和预测性信息进行重调度。数值实验表明，利用最新供应延期信息能提升预调度针对性，减少偏差成本。设计了嵌入基于Levy飞行粒子更新机制的改进粒子群算法，求解质量和效率高。后续将研究更多干扰模式下的飞机装配线动态调度问题以及各类干扰管理方法的比较和选择。

作者贡献声明

陆志强：提出研究问题，指导研究方案，论文审阅；

陆晨瑶：设计研究方案，建立数学模型并进行仿真实验，论文撰写与修订。