考虑质量衰退的批量调度与预防性维护的联合优化

摘要

针对生产调度过程中由于设备退化引起产品质量劣化的问题，构建了考虑质量与设备状态之间的耦合关系以及生产批次可分的批量流调度与预防性维护的联合优化模型。假设设备存在“受控”和“失控”两种状态，利用设备“失控”状态下产生的不合格品作为媒介描述上游设备与下游设备之间的随机依赖性，以最小化期望完工时间以及总成本为优化目标建立了数学模型；设计改进的差分进化算法对模型进行求解；通过数值实验验证了算法的收敛性，并通过对比实验对模型性能进行分析，结果显示了该联合优化模型的有效性。

关键词

预防性维护; 随机依赖性; 批量流调度; 改进差分进化算法; 集成优化

传统的生产调度模型通常假设设备一直可用^［1-2］，但在实际生产中，设备会随着使用逐渐老化。对此，常采用预防性维护来提高设备可靠度，降低因设备故障带来的生产损失。然而，过多或过少的维护，都会对生产效益产生不良影响。因此，考虑生产和维护的相互影响，联合优化二者，才有更高的实用价值。

近年来，生产调度和维护计划的集成问题逐渐成为制造领域研究的热点^［3］。针对单机调度，Cassady等^［4］建立了单机调度与预防性维护的集成模型；蒋志高等^［5］首次提出虚拟维护的概念，解决了加工时间可变的单机调度与维护问题；Cui等^［6］研究了具有柔性预防性维护周期的单机调度问题。流水车间调度中，Allaoui等^［7］假设设备必须在周期T内预防性维护一次；Miyata等^［8］提出维护级别的概念，并设计了有效的启发式算法；陆志强等^［9］研究了流水线调度与维护集成模型的鲁棒性，并提出了三阶段启发式算法。

然而，上述文献都忽略了质量在生产调度与维护计划二者集成问题中的影响。为了提高经济效益，还要注意到因设备老化导致的产品质量劣化带来的损失。Hadidi等^［10］对生产计划、调度、维护和质量的相关性模型和集成模型进行了综述；Pandey等^［11］针对单机系统，开发了一个联合优化生产计划、维护策略和质量控制的模型；Jamshidi等^［12］提出了一个混合整数非线性模型来优化质量、维护、提前/延误和中断的成本。质量和设备的影响是相互的，工件质量的退化也会影响设备的退化状态，例如汽车装配线上由于来料质量缺陷引起的定位工具的故障^［13］。这两者的耦合关系，使得设备间产生随机依赖性。随机依赖性多见于生产计划、维护与质量控制的集成研究中^［14-16］，在生产调度、维护和质量三者之间的集成研究中则十分罕见。

综上所述，生产调度、维护与质量的集成研究较少，已有的集成文献也很少考虑随机依赖性。因此，本文针对串行流水线系统，考虑以质量衰退工件为媒介的设备间的随机依赖性，同时结合现代多品种和中小批量的生产特点，引入批量流概念，分批生产和传输^［17］，建立生产调度和预防性维护的联合优化模型。

1 问题描述及数学建模

1.1　问题描述与假设

考虑一个串行流水线系统，在其调度期内，共有$n$个加工批次，需按照相同的工艺路线依次通过$m$台设备进行加工。每个加工批次可等量分成若干个子批次，每个子批次在一台设备上加工完成即可作为独立的转移批量转移到下一台设备，使得同一加工批次的前后工序交叉进行，减少设备的空闲时间，加快工件在设备间的流动。每台设备加工完一个子批次后进行一次设备检测，同时所有产品在系统末端进行质量检测，合格品用于满足需求，不合格品则需要返工。随着生产运行，各设备会随着使用逐渐发生退化，相应的加工工件的质量也会随之下降，设备间存在随机依赖性。为保证系统的有效输出，需要制定合适的维护策略。为了更清晰地描述研究问题，给出如下假设：

（1） 在任一时刻，一台设备只可加工一个子批次，且一个子批次也仅能在一台设备上进行加工，不同加工批次的子批次不可交叉生产。

（2） 考虑划分子批次产生的能耗成本，各加工批次划分的最大子批次数为2。

（3） 所有工件均在零时刻到达，且设备在零时刻均可用。

（4） 设备分为受控和失控两个状态。当设备处于受控状态时，其生产的产品均合格，当设备的退化状态超过给定的失效阈值时，设备将失控，产生质量失效，开始以恒定的百分比产生不合格品。

（5） 只有通过检测才能发现设备的退化状态，检测时间忽略不计。若设备状态超过其给定的失效阈值，则进行大修，使得设备修复如新。否则，以预防性维护阈值为决策方式，超过它则进行预防性维护，设备修复非新，不考虑设备的随机失效。

（6） 质量检测仅在系统末端进行，对不合格品只考虑返工成本，不占用所研究系统的生产能力。

1.2　符号定义

模型中用到的参数及相应定义如表1所示。

表1 模型参数及其定义

Tab.1 The definition of model parameters

参数	定义	参数	定义
$J_j$	加工批次编号，$j=\mathrm{1,2},...,n$	$z_{\left[i\right]}$	第$i$个位置的加工批次被划分成的子批次数
$M_k$	加工设备编号，$k=\mathrm{1,2},...,m$	$q_{\left[i\right]}$	第$i$个位置的加工批次的大小
$i$	加工批次顺序编号，$i=\mathrm{1,2},...,n$	$S_{0,\left[o\right],k}$	设备$M_k$第$o$个位置的子批次的计划开始时间
$o$	子批次顺序编号	$S_{\mathrm{w},\left[o\right],k}$	情景$\mathrm{w}$下，设备$M_k$第$o$个位置的子批次的实际开始时间
$t$	设备役龄变量	$C_{\mathrm{w},\left[o\right],k}$	情景$\mathrm{w}$下，设备$M_k$第$o$个位置的子批次的实际结束时间
$D_k\left(t\right)$	设备$M_k$在役龄$t$时的退化状态	$F_{\left[o\right],k}$	设备$M_k$在加工第$o$个子批次中的预测失效概率
$\mathrm{\Delta }{D}_k\left(t\right)$	时间$t$内$D_k\left(t\right)$的退化增量	$O_{\left[o\right]}$	加工第$o$个子批次时，系统的实际产出
$L_{k\mathrm{f}}$	设备$M_k$的失效阈值	${t_{L,}}_{\left[o\right],k}$	加工第$o$个子批次初，设备$M_k$的役龄
$L_{\mathrm{p}\mathrm{m}}$	预防性维护阈值	$t_{R,\left[o\right],k}$	加工第$o$个子批次末，设备$M_k$的役龄
${\alpha }_k\left(t\right),{\beta }_k$	Gamma过程的形状方程和尺度参数	$T_{\left[o\right],k}$	设备$M_k$加工第$o$个子批次的加工时间
$u_k,v_k$	Beta分布的参数	$T_{\mathrm{P}\mathrm{M}}$	设备单次预防性维护时间
$a_k$	形状方程${\alpha }_k\left(t\right)$的参数	$T_{\mathrm{O}\mathrm{H}}$	设备单次大修时间
${\rho }_{o,k}\left(t\right)$	加工第$o$个子批次时设备$M_k$的回归函数	$c_{\mathrm{e}\mathrm{c}}$	单次划分的能耗成本
$s_{o,k\mathrm{\text{'}},k}$	加工第$o$个子批次时设备$M_{k\mathrm{\text{'}}}$对$M_k$退化影响的加速因子	$c_{\mathrm{I}\mathrm{n}\mathrm{s}}$	单位设备状态检测成本
${\mu }_k$	设备$M_k$在不受控状态下产生不合格品的百分比	$c_{\mathrm{P},k}$	设备$M_k$的单位预防性维护成本
${\omega }_1,{\omega }_2$	优化目标的权重因子	$c_{\mathrm{O},k}$	设备$M_k$的单位大修成本
$x_{i,j}$	加工批次顺序的二元决策变量	$c_{\mathrm{s}\mathrm{c}\mathrm{r}}$	单位产品质量筛查成本
${\mathrm{\hslash }}_j$	表示加工批次$J_j$是否划分子批次的二元变量	$c_{\mathrm{p}\mathrm{e}\mathrm{n}}$	单位不合格品返工成本
$p_{j,k}$	设备$M_k$第$i$个位置的加工批次中产品的单位加工时间	$E_{\mathrm{C}}$	能耗总成本
$z_j$	加工批次$J_j$被划分成的子批次个数	$M_{\mathrm{C}}$	维护总成本
$q_j$	加工批次$J_j$的大小	$Q_{\mathrm{C}}$	质量总成本
$p_{\left[i\right],k}$	设备$M_k$第$i$个位置的加工批次中产品的单位加工时间

1.3　问题建模

（1）　维护策略模型

使用Gamma过程来对设备的退化状态进行建模。令$D_k\left(t\right)$表示设备$M_k$在役龄$t$时的退化状态，则设备退化状态$D_k\left(t\right)$在一段时间$\iota$内的非负增量$\mathrm{\Delta }{D}_k\left(t\right)$服从Gamma分布$Ga\left(\left.x\right|\left.{\alpha }_k\left(t\right),{\beta }_k\right)\right.$，其概率密度函数为

$$f_{{\alpha }_k\left(t\right),{\beta }_k}\left(x\right)=\frac{{{\beta }_k}^{{\alpha }_k\left(t\right)}{x}^{{\alpha }_k\left(t\right)-\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\beta }_{k}x\right)}{\Gamma \left({\alpha }_k\left(t\right)\right)}$$

设备每加工完一个子批次后进行状态检测并实施相应的维护。若设备处于失控状态，则刚加工的子批次需在系统末端进行质量筛查。若设备处于受控状态，则对应的子批次无需进行检测。但为了防止在下一个子批次加工期间出现质量失效，采用预测维护策略，将下一批次加工期间的预测失效概率作为预防性维护阈值。考虑到每台设备设立不同的维护阈值会使得求解难度随着系统内设备数的增加而上升，因而假设所有设备都采用相同的预防性维护阈值$L_{\mathrm{p}\mathrm{m}}$。则在加工第$o$个子批次中的预测失效概率可以表示为

$$\begin{array}{l}{F}_{\left[o\right],k}=P\left[\left.D_k\left(t_{R,\left[o\right],k}\right)\ge L_{k\mathrm{f}}\right|D_k\left(t_{R,\left[o-\mathrm{1}\right],k}\right)\right]=\\ P\left[\left.D_k\left(t_{R,\left[o-\mathrm{1}\right],k}\right)+\mathrm{\Delta }{D}_k\left(T_{\left[o\right],k}\right)\ge L_{k\mathrm{f}}\right|D_k\left(t_{R,\left[o-\mathrm{1}\right],k}\right)\right]=P\left[\mathrm{\Delta }{D}_k\left(T_{\left[o\right],k}\right)\ge L_{k\mathrm{f}}-D_k\left(t_{R,\left[o-\mathrm{1}\right],k}\right)\right]=\\ {\int }_{L_{k\mathrm{f}}-D_k\left(t_{R,\left[o-\mathrm{1}\right],k}\right)}^{\mathrm{\infty }}{f}_{{\alpha }_k\left(T_{\left[o\right],k}\right),{\beta }_k}\left(x\right)\mathrm{d}x=\\ \frac{\Gamma \left[{\alpha }_k\left(T_{\left[o\right],k}\right),{\beta }_k\left(L_{k\mathrm{f}}-D_k\left(t_{R,\left[o-\mathrm{1}\right],k}\right)\right)\right]}{\Gamma \left({\alpha }_k\left(T_{\left[o\right],k}\right)\right)}\end{array}$$

若$F_{\left[o\right],k}\ge L_{\mathrm{p}\mathrm{m}}$，则设备$M_k$需在生产第$o$个子批次前执行预防性维护。维护后，设备$M_k$的退化状态$D_k$有所改善，修复后设备的退化状态$D_k^{\mathrm{\text{'}}}$遵循区间$\left[\mathrm{0},D_k\right]$的Beta分布，其概率密度函数为

$$g_k\left(x\right)=\frac{\mathrm{1}}{D_k}\frac{\Gamma \left(u_k+v_k\right)}{\Gamma \left(u_k\right)\Gamma \left(v_k\right)}{\left(\frac{x}{D_k}\right)}^{u_k-\mathrm{1}}{\left(\mathrm{1}-\frac{x}{D_k}\right)}^{v_k-\mathrm{1}}\cdot {\mathrm{1}}_{\left\{\mathrm{0}\le x\le D_k\right\}}$$

式中：$u_k$和$v_k$可通过极大似然估计从历史数据中获取；${\mathrm{1}}_{\left\{\cdot \right\}}$是示性函数；$A$为真时${\mathrm{1}}_{\left\{A\right\}}=\mathrm{1}$，$A$为假时${\mathrm{1}}_{\left\{A\right\}}=\mathrm{0}$。

来料质量和设备退化的相互关系导致设备间的随机依赖性，使用加速因子$s_{o,k\mathrm{\text{'}},k}\left(s_{o,k\mathrm{\text{'}},k}\ge \mathrm{0}\right)$来描述第$o$个子批次从上游设备$M_{k\mathrm{\text{'}}}$输出对下游设备$M_k$退化的影响，若$s_{o,k\mathrm{\text{'}},k}=\mathrm{0}$，则表示设备$M_{k\mathrm{\text{'}}}$和$M_k$之间互不影响；否则设备$M_{k\mathrm{\text{'}}}$的失控将加速设备$M_k$退化。使用Gamma过程的形状函数来对该种随机依赖性建模，则有

$${\alpha }_k\left(t\right)=\left\{\begin{array}{c}{a}_{k}t,k=\mathrm{1}\\ a_{k}t\mathrm{e}\mathrm{x}\mathrm{p}\left[{\sum }_{i=\mathrm{1}}^{k-\mathrm{1}}\left({\rho }_{o,k\mathrm{\text{'}}}\left(t\right)s_{o,k\mathrm{\text{'}},k}\right)\right],k>\mathrm{1},k\mathrm{\text{'}}<k,\forall j\end{array}\right.$$

式中：$a_{k}t$是基准形状函数；回归函数${\rho }_{o,k\mathrm{\text{'}}}\left(t\right)$为上游设备是否失控的判断变量，其表达式可以写成

$${\rho }_{o,k\mathrm{\text{'}}}\left(t\right)=\left\{\begin{array}{c}\mathrm{0},D_k\left(t\right)<L_{k\mathrm{f}}\\ \mathrm{1},D_k\left(t\right)\ge L_{k\mathrm{f}}\end{array}\right.,k\mathrm{\text{'}}=\mathrm{1,2},...,k-\mathrm{1}$$

（2）　批量流调度与维护决策的集成模型

以子批次为对象进行调度，则总共有${\sum }_{j=\mathrm{1}}^n{z}_j$个子批次需要进行排序，且

$$z_j=\left\{\begin{array}{c}\mathrm{2},{\mathrm{\hslash }}_j=\mathrm{1}\\ \mathrm{1},{\mathrm{\hslash }}_j=\mathrm{0}\end{array}\right.$$

式中：二元变量${\mathrm{\hslash }}_j=\mathrm{1}$时表示划分子批次，否则为不划分。相同加工批次的子批次间不存在竞争排序关系。子批次顺序$o$与加工批次顺序$i$之间的关系满足$⌈\frac{o}{{\sum }_i{z}_{\left[i\right]}}⌉=\mathrm{1}$，$⌈x⌉$表示将$x$的值向上取整。设备$M_k$加工第$o$个子批次的加工时间$T_{\left[o\right],k}$为

$$T_{\left[o\right],k}=p_{\left[i\right],k}{q}_j/z_j,\forall o,k$$

用条件$A$和$B$分别表示加工第$o$个子批次前，设备$M_k$进行大修和预防性维护需要满足的条件，则有

$A$：$D_k\left(t_{R,\left[o-\mathrm{1}\right],k}\right)\ge L_{k\mathrm{f}}$；

$B$：$F_{\left[o\right],k}\ge L_{\mathrm{p}\mathrm{m}}$且$D_k\left(t_{R,\left[o-\mathrm{1}\right],k}\right)<L_{k\mathrm{f}}$。

综上，生产相关约束为

$${\sum }_{i=\mathrm{1}}^n{x}_{i,j}=\mathrm{1},j=\mathrm{1,2},...,n$$

（1）

$${\sum }_{j=\mathrm{1}}^n{x}_{i,j}=\mathrm{1},i=\mathrm{1,2},...,n$$

（2）

$$p_{\left[i\right],k}={\sum }_{j=\mathrm{1}}^n{x}_{i,j}{p}_{j,k},\forall i,k$$

（3）

$$z_{\left[i\right]}={\sum }_{j=\mathrm{1}}^n{x}_{i,j}{z}_j,\forall i$$

（4）

$$q_{\left[i\right]}={\sum }_{j=\mathrm{1}}^n{x}_{i,j}{q}_j,\forall i$$

（5）

$$S_{\mathrm{0},\left[o\right],k}\ge S_{\mathrm{0},\left[o-\mathrm{1}\right],k}+p_{\left[i\right],k}{q}_{\left[i\right]}/z_{\left[i\right]}+{\mathrm{1}}_{\left\{A\right\}}{T}_{\mathrm{O}\mathrm{H}}+{\mathrm{1}}_{\left\{B\right\}}\cdot T_{\mathrm{P}\mathrm{M}},\forall k,o>\mathrm{1}$$

（6）

$$S_{\mathrm{0},\left[o\right],k}\ge S_{\mathrm{0},\left[o\right],k-\mathrm{1}}+p_{\left[i\right],k-\mathrm{1}}{q}_{\left[i\right]}/z_{\left[i\right]},\forall o,k>\mathrm{1}$$

（7）

$$S_{\mathrm{0},\left[\mathrm{1}\right],\mathrm{1}}=\mathrm{0}$$

（8）

$$C_{\mathrm{w},\left[o\right],k}\ge S_{\mathrm{w},\left[o\right],k}+p_{\left[i\right],k}{q}_{\left[i\right]}/z_{\left[i\right]},\forall k,o$$

（9）

$$S_{\mathrm{w},\left[o\right],k}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{C_{\mathrm{w},\left[o-\mathrm{1}\right],k}+{\mathrm{1}}_{\left\{A\right\}}{T}_{\mathrm{O}\mathrm{H}}+{\mathrm{1}}_{\left\{B\right\}}\cdot T_{\mathrm{P}\mathrm{M}},C_{\mathrm{w},\left[o\right],k-\mathrm{1}}\right\},\forall k>\mathrm{1},o>\mathrm{1}$$

（10）

$$S_{\mathrm{w},\left[o\right],\mathrm{1}}=C_{\mathrm{w},\left[o-\mathrm{1}\right],\mathrm{1}}+{\mathrm{1}}_{\left\{A\right\}}{T}_{\mathrm{O}\mathrm{H}}+{\mathrm{1}}_{\left\{B\right\}}{T}_{\mathrm{P}\mathrm{M}},\forall o>\mathrm{1}$$

（11）

$$S_{\mathrm{w},\left[\mathrm{1}\right],k}=C_{\mathrm{w},\left[\mathrm{1}\right],k-\mathrm{1}},\forall k>\mathrm{1}$$

（12）

$$S_{\mathrm{w},\left[\mathrm{1}\right],\mathrm{1}}=\mathrm{0}$$

（13）

$$x_{i,j}\in \left\{\mathrm{0,1}\right\}$$， $$\forall i,j=\mathrm{1,2},...,n$$

（14）

$$p_{\left[i\right],k},z_{\left[i\right]},q_{\left[i\right]},S_{\mathrm{0},\left[o\right],k},S_{\mathrm{w},\left[o\right],k},C_{\mathrm{w},\left[o\right],k}\ge \mathrm{0},\forall i,k,o$$

（15）

约束（1）~（2）表示加工批次和加工位置的一一对应关系。约束（3）确定各设备按加工顺序的单位加工时间。约束（4）~（5）确定了按加工顺序的子批次数量及加工批次的大小。约束（6）~（13）为子批次加工时间的约束，当需要维护时通过右移策略进行重调度，求得实际加工时间$S_{\mathrm{w},\left[o\right],k}$。约束（14）~（15）为变量的取值约束。

（3）　目标函数计算

同时优化时间和成本两个目标。时间相关的优化目标为最小化完工时间的期望值，即$f_{\mathrm{1}}=E_{\mathrm{w}}\left(C_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)$。$C_{\mathrm{m}\mathrm{a}\mathrm{x}}$为最后一个产品的实际结束时间。

成本相关的优化目标为最小化能耗成本、维护成本及质量成本的期望值之和。划分子批次会产生设备准备能耗和传输能耗，能耗成本$E_{\mathrm{C}}$计算为

$$E_{\mathrm{C}}=c_{\mathrm{e}\mathrm{c}}{\sum }_{j=\mathrm{1}}^n{\mathrm{\hslash }}_j,j=\mathrm{1,2},...,n$$

维护相关成本包括系统检测成本、预防性维护成本和大修成本，则维护相关成本$M_{\mathrm{C}}$可以计算为

$$M_{\mathrm{C}}=mo{c}_{\mathrm{I}\mathrm{n}\mathrm{s}}+{\displaystyle \sum _{k=\mathrm{1}}^m}\sum _{o\ge \mathrm{2}}{c}_{\mathrm{P},k}{\mathrm{1}}_{\left\{B\right\}}+{\displaystyle \sum _{k=\mathrm{1}}^m}\sum _{o\ge \mathrm{2}}{c}_{\mathrm{O},k}{\mathrm{1}}_{\left\{A\right\}}$$

令$S$为串行系统中所有设备的集合，用$S_{\mathrm{F},\left[o\right],e}$表示在加工第$o$个子批次时，生产第$e\left(e=\mathrm{1,2},...n\right)$个产品时所有出现质量故障的设备的集合，$p_{\left[o\right],k}$为设备$M_k$加工第$o$个子批次的单位时间，则$S_{\mathrm{F},\left[o\right],e}=\left\{\left.k\right|D_k\left(t_{R,\left[o-\mathrm{1}\right],k}\right)+\mathrm{\Delta }{D}_k\left(e{p}_{\left[o\right],k}\right)>L_{k\mathrm{f}},k\in S\right\}$。那么，第$e$个产品为合格品的概率为${\prod }_{k\in S}\left(\mathrm{1}-{\mu }_k{\mathrm{1}}_{\left\{k\in S_{\mathrm{F},\left[o\right],e}\right\}}\right)$，则加工第$o$个子批次时，系统的实际产出为$O_{\left[o\right]}=\sum _o\prod _{k\in S}\left(\mathrm{1}-{\mu }_k{\mathrm{1}}_{\left\{k\in S_{\mathrm{F},\left[o\right],e}\right\}}\right)$。

质量相关的成本$Q_{\mathrm{C}}$包括产品质量筛查成本和不合格品返工成本，则有

$$Q_{\mathrm{C}}=\left(\left.{\sum }_o{c}_{\mathrm{s}\mathrm{c}\mathrm{r}}{q}_{\left[i\right]}/z_{\left[i\right]}\right|\exists k使得{\mathrm{1}}_{\left\{D_k\left(t_{R,\left[o\right],k}\right)\ge L_{k\mathrm{f}}\right\}}=\mathrm{1}\right)+{\sum }_o{c}_{\mathrm{p}\mathrm{e}\mathrm{n}}\left(q_{\left[i\right]}/z_{\left[i\right]}-O_{\left[o\right]}\right)$$

式中：$\left(\left.{\sum }_o{c}_{\mathrm{s}\mathrm{c}\mathrm{r}}{q}_{\left[i\right]}/z_{\left[i\right]}\right|\exists k使得{\mathrm{1}}_{\left\{D_k\left(t_{R,\left[o\right],k}\right)\ge L_{k\mathrm{f}}\right\}}=\mathrm{1}\right)$表示对任一子批次$o$，只要存在一台设备在加工完子批次$o$后，其状态超过给定的失效阈值，该子批次$o$就需要进行质量筛查。

基于以上计算，成本相关的优化目标可以表示为 $f_{\mathrm{2}}=E_{\mathrm{C}}+M_{\mathrm{C}}+Q_{\mathrm{C}}$。整合优化目标$f_{\mathrm{1}}$和$f_{\mathrm{2}}$，并统一量纲，引入权重因子${\omega }_{\mathrm{1}},{\omega }_{\mathrm{2}}$，则考虑时间和成本的多目标函数表达式为

$$\mathrm{m}\mathrm{i}\mathrm{n}f={\omega }_{\mathrm{1}}\frac{f_{\mathrm{1}}}{f_{\mathrm{1}}^{\mathrm{*}}}+{\omega }_{\mathrm{2}}\frac{f_{\mathrm{2}}}{f_{\mathrm{2}}^{\mathrm{*}}}$$

式中：${\omega }_{\mathrm{1}}+{\omega }_{\mathrm{2}}=\mathrm{1}$。$f_{\mathrm{1}}^{\mathrm{*}}$和$f_{\mathrm{2}}^{\mathrm{*}}$为单一目标最优值，即分别求解最小完工时间的期望值及最小相关成本值。

2 算法设计

本模型的求解是NP难问题，对此，设计了改进的差分进化算法。采用实值编码，$n$表示加工批次数，编码长度为$\left(n+\mathrm{1}\right)$，前$n$个编码用于确定加工批次的顺序及是否划分批量，取值范围为$\left(\mathrm{0,2}\right)$。整数部分用来确定是否划分批次；小数部分按升序排列加工批次的顺序。第$\left(n+\mathrm{1}\right)$个编码表示维护阈值的大小，取值范围为$\left(\mathrm{0,1}\right)$。给出一个实例来详细说明解码过程，如图1所示。

图1 解码过程示意图

Fig.1 Schematic diagram of the decoding process

接下来对算法求解步骤进行详细描述：

（1）种群$X\left(t\right)=\left\{\left.X_I\left(t\right)\right|I=\mathrm{1,2},...,I_{\mathrm{m}\mathrm{a}\mathrm{x}}\right\}$初始化，$I_{\mathrm{m}\mathrm{a}\mathrm{x}}$为种群中的个体数量，计算种群的适应度值，找到最优值对应的个体$X^{\mathrm{*}}\left(t\right)$，适应度值$f\left(X\left(t\right)\right)$通过蒙特卡洛仿真方法计算得到。

（2）变异操作

为了提高求解精度与加快收敛速度，利用鲸鱼算法^［18］的搜索机制来进行变异。同时，考虑到标准鲸鱼算法存在容易陷入局部最优、收敛精度低的问题，引入了柯西变异算子^［19］和非线性变化收敛因子^［20］。

对种群中的每个个体$X_I\left(t\right)$进行变异得到$X_{\mathrm{V}\mathrm{a},I}\left(t+\mathrm{1}\right)$，$Va\left(t+\mathrm{1}\right)=\left\{\left.X_{\mathrm{V}\mathrm{a},I}\left(t+\mathrm{1}\right)\right|I=\mathrm{1,2},...,I_{\mathrm{m}\mathrm{a}\mathrm{x}}\right\}$表示变异后得到的种群，$X_{\mathrm{V}\mathrm{a},I}\left(t+\mathrm{1}\right)$的计算与一个$\left[\mathrm{0,1}\right]$上的随机数$p$相关，即：

当$p<\mathrm{0.5}$时，$X_{\mathrm{V}\mathrm{a},I}\left(t+\mathrm{1}\right)=\left\{\begin{array}{c}{X}^{\mathrm{*}}\left(t\right)-A\left|C{X}^{\mathrm{*}}\left(t\right)-X_I\left(t\right)\right|,\left|A\right|<\mathrm{1}\\ X_I\left(t\right)+A\mathrm{t}\mathrm{a}\mathrm{n}\left(\left(r-\mathrm{1}/\mathrm{2}\right)\mathrm{\pi }\right),\left|A\right|\ge \mathrm{1}\end{array}\right.$

当$p\ge \mathrm{0.5}$时，$X_{\mathrm{V}\mathrm{a},I}\left(t+\mathrm{1}\right)=\left|X^{\mathrm{*}}\left(t\right)-X_I\left(t\right)\right|\cdot e^{bl}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\mathrm{\pi }l\right)+X^{\mathrm{*}}\left(t\right)$

式中：$r$为$\left[\mathrm{0,1}\right]$上的均匀分布；$b$为定义对数螺旋线形状的常数；$l$为区间$\left[-\mathrm{1,1}\right]$上的随机数。此外，$A=\mathrm{2}a{r}_{\mathrm{1}}-a$以及$C=\mathrm{2}{r}_{\mathrm{2}}$，$r_{\mathrm{1}}$和$r_{\mathrm{2}}$为区间$\left[\mathrm{0,1}\right]$上的随机数；$a$是收敛因子，影响算法的搜索精度，非线性变化收敛因子随迭代更新的公式如下^［20］：

$$a=N+\frac{\mathrm{1}-t/t_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{1}-\mu t}$$

式中：$N$为收敛因子基准常量；$t$为当前迭代次数；$t_{\mathrm{m}\mathrm{a}\mathrm{x}}$为最大迭代次数；$\mu \in \left(\mathrm{1}/t_{\mathrm{m}\mathrm{a}\mathrm{x}},\mathrm{1}\right)$为非线性调整系数。

（3）交叉操作

变异操作结束后，采用二项式交叉，表达式如下：

$$X_{\mathrm{M}\mathrm{u},I}\left(t+\mathrm{1}\right)=\left\{\begin{array}{c}{X}_{\mathrm{V}\mathrm{a},I}\left(t+\mathrm{1}\right),\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}<C_{R}or\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{i}\left(\mathrm{1},l_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)=I\\ X_I\left(t\right),\mathrm{否}\mathrm{则}\end{array}\right.$$

式中：$X_{\mathrm{M}\mathrm{u},I}\left(t+\mathrm{1}\right)$交叉后的个体；$C_{\mathrm{R}}$是交叉概率；$M_{\mathrm{u}}\left(t+\mathrm{1}\right)=\left\{\left.X_{\mathrm{M}\mathrm{u},I}\left(t+\mathrm{1}\right)\right|I=\mathrm{1,2},...,I_{\mathrm{m}\mathrm{a}\mathrm{x}}\right\}$为交叉后的种群。

（4）变异和交叉结束后，通过选择决定子代个体是否进入下一代。

（5）设置迭代上限作为算法终止条件，输出求解结果，否则返回步骤（2）。算法流程如图2所示。

图2 改进的差分算法流程

Fig.2 The improved differential evolution algorithm flow

3 算例分析

3.1　改进的差分进化算法与标准进化算法比较

考虑$m=\left\{\mathrm{3,5},\mathrm{8}\right\}$，$n=\left\{\mathrm{10,30,50}\right\}$共9种问题规模，表2给出了设备退化过程相关的参数，该生产系统的成本参数和时间参数如表3所示，假设加工时间$p_{j,k}$服从$\left[\mathrm{1,3}\right]$的均匀分布，加工批次大小$q_j$服从$\left[\mathrm{20,100}\right]$的均匀分布，权重因子${\omega }_{\mathrm{1}}={\omega }_{\mathrm{2}}=\mathrm{0.5}$。

表2 各设备的退化与维护参数

Tab.2 Parameters of different machines

参数	$L_{k\mathrm{f}}$	$D_k\left(0\right)$	${\mu }_k$	$c_{\mathrm{O},k}$($)	$c_{\mathrm{P},k}$($)	${\alpha }_k$	${\beta }_k$	$u_k$	$v_k$
取值	$U\left(\mathrm{5,20}\right)$	0	$U\left(\mathrm{0.1,0.6}\right)$	$U\left(\mathrm{150,300}\right)$	$U\left(\mathrm{50,100}\right)$	$U\left(\mathrm{0,0.1}\right)$	$U\left(\mathrm{0.5,1}\right)$	$U\left(\mathrm{0.2,2}\right)$	$U\left(\mathrm{1,7}\right)$

表3 其他成本和时间参数

Tab.3 Cost and time parameters in production

$c_{\mathrm{e}\mathrm{c}}$($)	$c_{\mathrm{I}\mathrm{n}\mathrm{s}}$($)	$c_{\mathrm{s}\mathrm{c}\mathrm{r}}$ ($)	$c_{\mathrm{p}\mathrm{e}\mathrm{n}}$($)	$T_{\mathrm{P}\mathrm{M}}$	$T_{\mathrm{O}\mathrm{H}}$
20	5	0.5	30	40	100

由于目标函数由系统多次仿真的期望值表示，因此需要确定仿真的抽样样本数。通过选取4个具有不同决策变量的测试实例通过大量仿真来确定合理的样本数量。根据最终的实验结果，仿真样本量确定为50。

在改进的差分进化算法（IDE）中，将种群大小设置为30，将最大迭代次数设置为50，并将交叉因子$C_{\mathrm{R}}$设置为0.4，标准差分进化算法（DE）变异因子设置为0.5，其余参数与IDE相同。确定所有参数后，对算例进行求解，绘制各规模算例下的IDE迭代曲线，如图3所示，IDE与DE的求解结果对比如图4所示。从图3可以看出，IDE算法具有收敛性，对中小规模的算例求解效果较好，而且由于仿真计算时间的限制，该算法也不适用于更大规模问题的求解。从图4可以看出，IDE相较于DE的求解结果较优，在设备数较少和规模小时求解效果更好。

图3 改进的差分进化算法的迭代曲线

Fig.3 Iterative curve of improved differential evolution algorithm

图4 改进的差分进化算法与标准差分进化算法的求解结果

Fig.4 The results of the improved differential evolution algorithm and the differential evolution algorithm

3.2　与其他维护模型比较

记上文已建立的模型为P0；不考虑预防性维护仅考虑大修的维护模型记为P1；每次检查后都进行一次预防性维护或者大修的维护模型记为P2以及忽略设备间的随机依赖性，单独求解每台设备的最优维护阈值的维护模型记为P3，将这几种维护模型进行比较，其对比结果如表4所示。

表4 不同维护模型下的结果

Tab.4 The results obtained under different maintenance models

$m\times n$		$f$	$f_1$	$f_2$($)	$E_{\mathrm{C}}$($)	$M_{\mathrm{C}}$($)			$Q_{\mathrm{C}}$($)
						$C_{\mathrm{I}\mathrm{n}\mathrm{s}}$	$C_{\mathrm{P}\mathrm{M}}$	$C_{\mathrm{O}\mathrm{H}}$	$C_{\mathrm{I}\mathrm{n}\mathrm{s}\mathrm{Q}}$	$C_{\mathrm{r}\mathrm{e}\mathrm{w}}$
3×10	P0	1.001	2 293.10	2 873.71	140	255	1 105.40	515.48	110.30	747.53
	P1	1.433	2 313.18	5 327.01	120	240	0	1 356.48	228.09	3 382.44
	P2	1.103	2 483.70	3 219.70	120	240	1 737.66	468.68	103.37	549.99
	P3	1.018	2 314.90	2 946.52	140	255	1 115.60	507.20	115.10	813.62
3×30	P0	1.046	6 271.40	6 449.46	340	705	1 047.58	1 539.24	401.35	2 416.29
	P1	1.393	6 062.08	11 008.29	360	720	0	2 810.00	554.74	6 563.55
	P2	1.213	6 238.40	8 575.82	100	525	2 817.80	1 214.16	478.49	3 440.37
	P3	1.076	6 407.90	6 683.56	420	765	1 211.80	1 448.20	367.70	2 470.86
3×50	P0	1.003	8 706.20	13 661.76	520	1 140	4 390.88	3 218.52	454.21	3 938.15
	P1	1.292	8 689.18	21 526.07	540	1 155	0	7 229.32	831.72	11 770.03
	P2	1.178	9 418.12	17 287.22	80	810	8 340.64	3 091.88	529.42	4 435.28
	P3	1.020	9 295.80	13 183.95	320	990	4 422.40	3 220.60	479.45	3 751.50
5×10	P0	1.009	2 934.24	8 106.06	160	450	2 853.28	2 229.00	189.60	2 224.18
	P1	1.237	2 982.26	11 715.97	200	500	0	4 901.28	265.59	5 849.10
	P2	1.064	2 916.26	9 092.99	80	350	3 814.72	1 950.48	221.82	2 675.97
	P3	1.046	2 873.40	8 932.41	80	350	2 408.10	2 778.00	232.50	3 083.81
5×30	P0	1.015	6 476.70	21 260.53	340	1 175	7 778.54	6 032.32	504.86	5 429.81
	P1	1.199	6 437.20	29 220.48	360	1 200	0	12 880.24	691.28	14 088.96
	P2	1.121	6 761.00	24 803.10	80	850	9 272.20	5 786.98	594.03	8 219.89
	P3	1.031	6 840.00	20 723.37	420	1 275	7 661.90	6 161.20	470.35	4 734.92
5×50	P0	1.011	10 926.44	36 961.82	560	1 950	12 710.60	8 432.80	893.41	12 415.01
	P1	1.201	10 608.62	52 142.11	680	2 100	0	18 588.60	1 236.02	29 537.49
	P2	1.102	11 291.00	42 477.62	500	1 875	20 056.90	7 430.84	870.73	11 744.15
	P3	1.026	11 154.14	37 278.18	460	1 825	13 301.26	8 321.60	922.83	12 447.49
10×10	P0	1.048	5 131.56	24 872.89	80	700	2 260.16	8 073.48	325.86	13 433.39
	P1	1.076	5 006.80	26 909.29	60	650	0	9 849.00	331.84	16 018.45
	P2	1.056	5 053.88	25 682.52	80	700	3 375.58	8 508.04	333.32	12 685.58
	P3	1.067	5 033.20	26 367.84	60	650	1 816.60	9 471.80	348.00	14 021.44
10×30	P0	1.044	9 271.36	55 241.46	380	2 450	7 132.46	16 543.10	829.34	27 906.56
	P1	1.181	9 577.68	67 689.47	100	1 750	0	22 160.90	854.20	42 824.37
	P2	1.077	9 819.30	55 480.96	380	2 450	11 868.36	14 484.42	810.09	25 488.09
	P3	1.059	9 817.20	53 597.79	380	2 450	8 346.30	15 432.80	827.20	26 161.49
10×50	P0	1.041	13 744.84	82 974.58	580	3 950	20 366.50	28 271.60	1 267.88	28 538.60
	P1	1.118	13 133.32	99 687.39	280	3 200	0	42 416.92	1 435.17	52 355.30
	P2	1.091	13 066.90	95 661.66	220	3 050	26 945.40	29 298.32	1 355.76	34 792.18
	P3	1.042	13 328.76	85 849.10	520	3 800	20 930.50	29 025.52	1 300.03	30 273.05

注：  $C_{\mathrm{I}\mathrm{n}\mathrm{s}}$-系统检测成本；$C_{\mathrm{P}\mathrm{M}}$-预防性维护成本；$C_{\mathrm{O}\mathrm{H}}$-大修成本；$C_{\mathrm{I}\mathrm{n}\mathrm{s}\mathrm{Q}}$-产品筛查成本；$C_{\mathrm{r}\mathrm{e}\mathrm{w}}$-不合格品返工成本

由表4可以看出模型P0的求解结果是最优的。通过对比发现，随着设备的增多，P0的大修成本占总成本的比例增加，这是因为设备越多，设备间的随机依赖关系可能使得下游设备的退化加速，造成设备失效的同时产生不合格品，返工成本的占比增加也侧面印证了这一点。P1的总维护成本高于P0的总维护成本，说明安排预防性维护能够降低维护成本；此外，P1的返工成本远大于P0的返工成本，这是由于设备失效次数增多导致不合格品数量增加。P2的总维护成本总是高于P0的总维护成本表明过度维护也会增加维护成本。P3的求解结果仅次于P0，说明忽视设备间的随机依赖性也会使得求解效果不佳。

3.3　与其他生产模型比较

将不进行划分批次的模型记为P4；将单独决策模型，即先确定生产顺序及批次划分，再制定维护策略的模型记为P5，把二者的求解结果与P0的求解结果进行对比，如表5所示。

表5 不同生产模型下的结果

Tab.5 The results obtained under different production models

$m\times n$		$f$	$f_1$	$f_2$($)	$E_{\mathrm{C}}$($)	$M_{\mathrm{C}}$($)			$Q_{\mathrm{C}}$($)
						$C_{\mathrm{I}\mathrm{n}\mathrm{s}}$	$C_{\mathrm{P}\mathrm{M}}$	$C_{\mathrm{O}\mathrm{H}}$	$C_{\mathrm{I}\mathrm{n}\mathrm{s}\mathrm{Q}}$	$C_{\mathrm{r}\mathrm{e}\mathrm{w}}$
3×10	P0	1.001	2 293.10	2 873.71	140	255	1 105.40	515.48	110.30	747.53
	P4	1.167	2 457.32	3 620.19	0	150	825.74	730.60	210.44	1 703.41
	P5	1.111	2 576.00	3 154.75	120	240	957.60	658.60	147.00	1 031.55
3×30	P0	1.046	6 271.40	6 449.46	340	705	1 047.58	1 539.24	401.35	2 416.29
	P4	1.132	5 860.00	7 957.96	0	450	1 252.12	1 232.40	564.90	4 458.54
	P5	1.057	6 058.70	6 818.11	240	630	1 332.50	1 072.40	409.00	3 134.21
3×50	P0	1.003	8 706.20	13 661.76	520	1 140	4 390.88	3 218.52	454.21	3 938.15
	P4	1.043	8 953.64	14 347.60	0	750	4 105.20	3 442.98	625.54	5 423.88
	P5	1.015	8 972.40	13 555.38	320	990	3 825.20	3 688.60	535.50	4 196.08
5×10	P0	1.009	2 934.24	8 106.06	160	450	2 853.28	2 229.00	189.60	2 224.18
	P4	1.179	2 955.80	10 848.86	0	250	1 478.88	2 561.36	271.06	6 287.56
	P5	1.014	2 922.60	8 266.64	80	350	2 433.60	2 443.20	216.50	2 743.34
5×30	P0	1.015	6 476.70	21 260.53	340	1 175	7 778.54	6 032.32	504.86	5 429.81
	P4	1.063	6 119.48	24 517.25	0	750	4 701.06	7 392.20	686.68	10 987.31
	P5	1.045	6 761.60	21 581.26	300	1 125	6 125.90	6 779.60	540.10	6 710.66
5×50	P0	1.011	10 926.44	36 961.82	560	1 950	12 710.60	8 432.80	893.41	12 415.01
	P4	1.133	10 769.04	46 530.29	0	1 250	8 675.46	10 764.88	1 206.56	24 633.39
	P5	1.086	11 063.10	42 051.65	240	1 550	9 514.20	9 823.60	1 063.35	19 860.50
10×10	P0	1.048	5 131.56	24 872.89	80	700	2 260.16	8 073.48	325.86	13 433.39
	P4	1.058	4 942.36	26 382.70	0	500	1 326.64	8 491.84	346.90	15 717.32
	P5	1.049	4 769.40	26 872.20	120	800	2 331.00	10 767.20	339.85	12 514.15
10×30	P0	1.044	9 271.36	55 241.46	380	2 450	7 132.46	16 543.10	829.34	27 906.56
	P4	1.154	9 292.24	66 575.68	0	500	2 222.86	19 846.74	856.48	42 149.60
	P5	1.094	9 252.90	60 498.28	320	2 300	6 163.20	19 648.00	830.35	31 236.73
10×50	P0	1.041	13 744.84	82 974.58	580	3 950	20 366.50	28 271.60	1 267.88	28 538.60
	P4	1.064	12 808.48	92 869.69	0	2 500	13 762.94	32 565.28	1 430.94	42 610.53
	P5	1.045	13 470.64	85 387.29	440	3 600	16 358.70	30 201.92	1 333.31	33 453.36

由表5可以发现， P5的求解结果次优于P0的求解结果，表明了集成模型较单独决策模型的优越性。P4的质量成本远大于P0的质量成本，原因是不划分批次的情况下，由于批次加工无法中断，且批次较大，设备在加工中失效的次数增多，相应的不合格品的数量也随之增加。选取一个5×10的算例分析批次划分次数对模型求解结果的影响，如图5所示，可以看出，随着批次划分次数的增多，期望完工时间不规律变化，此外，维护成本呈增加的趋势、质量成本呈下降趋势，但总成本是趋于减小的，说明适量选择划分批次次数才能同时优化总成本和期望完工时间。

图5 批次划分次数对求解结果的影响

Fig.5 Influence of splitting times on the result

4 结语

本文考虑产品质量与设备退化状态之间的相互影响，以最小化完工时间和成本为目标，建立了考虑设备间随机依赖性的生产调度与维护的集成模型。此外，还引入了批量流的概念，通过算例分析得知适量划分批次可以优化总成本和期望完工时间。设计改进的差分进化算法对模型进行求解，数值实验的结果表明，该模型可有效规划生产排程、进行维护决策；并通过对比试验显示了该集成模型较单独决策模型的优越性。

本文以“批次在设备上加工完成”这一事件来触发设备的检测，存在检测后不进行维护的情况，即设备检测较为频繁，未来研究可以考虑对检测计划进行研究，以节省检测成本。此外，还可以考虑每台设备的加工顺序不同的状况，不等量分批的情形，以及引入机会维护的可能。

作者贡献声明

陆志强：对论文进行指导工作。

牟小涵：论文的撰写、模型建立、实验分析、数据处理的工作。