基于能源效率的生产系统节能延迟窗机会维修模型

摘要

将系统有效产出与能源消耗结合起来，提出了能源效率指标（energy efficiency indicator， EEI），构建了联合维修模型。以多设备串联系统为研究对象，令各设备退化服从Tweedie指数散布过程，考虑了生产过程中的产出损失，用基于节能延迟窗口的机会维修策略，令传统时间窗中最早到达最优维修时刻决策点成为维修组合窗决策点，以最小化系统在规划期内的能源效率（energy efficiency， EE）为目标构建机会维修模型。最后通过蒙特卡洛仿真实验进行算例分析，结果表明在减少系统总能源消耗和提高系统有效产出方面，该机会维修模型切实可行。

关键词

机会维修; 能源效率; 节能延迟时间窗; Tweedie指数散布过程; 产出损失

复杂的制造生产活动需要消耗大量能源，然而随着能源成本上升，全球制造业竞争力增强，环境问题日益突出，能源在现代制造业扮演着愈加重要的角色^［1-2］。研究表明制造车间内有85%的能源都被消耗在生产无关的活动上，如机器堵塞、饥饿等，节能减耗势在必行^［3］。串联生产线中一旦有设备发生故障，将导致整个系统停机，因此预防性维修对提高生产系统效率作用日益重要^［4］。

早期维修方法中，大多数制造系统预防性维修研究主要基于可靠性等指标触发，目标集中于降低成本，减少质量损失，从而提高生产速率。Radhoui等^［5］针对考虑缓冲区的生产系统，开发了一个联合质量控制和预防性维护的模型。将不合格品数量与阈值比较，决定是否采取维护措施，并以单位时间内总成本最小为目标，利用生成的二阶响应面确定最优生产速率。能源消耗作为重要因素通常被忽略。在可持续发展理念下，能源消耗及能源利用效率在制造企业决策制定中占据越来越重要的地位^［6］。开展有效维护活动不仅可以减少因突发性停机造成的生产损失，提高设备可靠度，而且可以使设备恢复至较良好的状态，从而提高其能源利用效率。与此同时，能源效率（energy efficiency）作为一类新兴的关键绩效指标，开始在维修活动中受到广泛关注^［7］。Seow等^［8］提出能耗最小化设计（DFEM），在此基础上建立能源仿真模型，进而提高生产过程中能源利用效率，并为设计人员提供参考；Xia等^［9］提出了一种以能源为导向的维修方法，使生产线以节能的方式运行。在设备层对能源属性建模，扩展多属性模型。在系统层，则提出了一种新的节能窗口策略，目的是降低整个生产线的能耗；Zhou等^［10］将系统有效产出与能源消耗结合起来，引入节能机会窗维修机制，提出了一种新的能耗效率指标EEE（effective energy efficiency）。

目前，静态成组维修策略通过考虑制造系统在无限决策时间内稳定运行，根据长期运行积累的系统历史信息，制定的静态维修规则在决策时间内保持不变，这存在一定程度的不合理。苏春等^［11］采用滚动窗口方法完成维修活动的动态调度，以随机停机事件作为风力机成组维修的机会，并动态更新原有的维修计划；在俞梦琦等^［12］的机会维修决策模型中，以单设备达到最优维修周期作为其余设备提前维修的判断条件，将传统单窗策略扩展为双时间窗策略；在Lin等^［13］的研究中可以发现，在一定范围内将最先达到最佳维修周期设备的维修活动适当延迟，对当前周期内的该设备成本影响比较小，还可有效降低系统在规划期内整体的维修成本。

在目前的制造系统预防性维修文献中，大多在系统层和设备层进行独立研究，仅针对能耗或产出某一个因素展开研究，且大多数成果都是基于单一的随机过程建立预防性维护模型。然而，在实际生产中，每一个随机过程都有其特定的应用范围，随着设备退化机制的复杂化和产品功能的多样化，由单一随机过程构建的退化维护决策模型开始出现准确性与适用性不佳等问题。

目前少有论文结合Tweedie指数分布，并联合设备层和系统层对可持续发展目标开展研究。因此，本文针对现有研究的不足，面向多设备连续生产串联系统，令各设备的随机退化服从Tweedie指数散布过程，并提出能源效率指标EEI，通过引入设备能耗、产出和次品率等要素，求解系统层的目标函数，提出节约能耗型基于延迟时间窗口策略机会维修模型，在有效降低多设备随机退化系统的能耗同时并减少系统总生产损失。

1 问题描述

数学模型各参数见表1。如图1所示，是一个由$n$台设备组成的串联生产系统，其中第$l$台设备为额定加工速率最慢的设备$M_l$。为进一步研究问题域，假设：

表1  数学模型参数

Tab.1  Mathematical model parameter

参数	释义
$M_i$	串联生产系统中的第$i$个设备，$i=\mathrm{1,2},...,n$
$\alpha ,\beta$	Tweedie分布的偏移参数和散布参数
$D_i^{\mathrm{f}}$	设备$M_i$的失效阈值
$D_i^{\mathrm{p}}$	设备$M_i$的预防性维护阈值
$w_i^{\mathrm{p}}$	设备$M_i$的预防性维护时间
$t_{i,j}^{\mathrm{w}}$	设备$M_i$在第$j$个维护周期内的运行时间
$w_{i,j}^{\mathrm{s}}$	设备$M_i$在第$j$次更换时耗费时间
$T_{i,j}$	设备$M_i$的第$j$个维护周期
$t_{i,j}$	设备$i$进行第$j$次维护的时刻
$t_k^{\mathrm{*}}$	维护决策触发时刻
$T$	系统生产计划期
$k_i^{\mathrm{p}}$	维护周期个数
$W(X(t))$	产出损失
$W$	设备磨损量
$e_i^{\mathrm{w}}$	设备在正常运行时的单位能耗
$e_i^{\mathrm{p}}$	设备$M_i$在单位时间内预防性维护能耗
$E_{i,j}^{\mathrm{w}}$	设备$M_i$在第$j$个维护周期内的运行能耗
$E_{i,j}^{\mathrm{p}}$	设备$M_i$在第$j$个维护周期内的预防性维护能耗
$E_{i,j}^{\mathrm{s}}$	设备$M_i$在第$j$个维护周期内的完美更换能耗
$E_{\mathrm{s}\mathrm{y}\mathrm{s}}^{\mathrm{w}}$	系统运行能耗
$E_{\mathrm{s}\mathrm{y}\mathrm{s}}^{\mathrm{p}}$	预防性维护能耗
$E_{\mathrm{s}\mathrm{y}\mathrm{s}}^{\mathrm{s}}$	完美更换能耗
决策变量	释义
$\mathrm{\Delta }O$	设备维护组合的时间窗
$\mathrm{\Delta }o$	决策维护时机的时间窗

图1  串行生产线

Fig.1  Serial production line

⑴生产线连续生产，各设备生产过程相对稳定，当某台设备停机时，整个生产线停止。

⑵各台设备在加工的过程中逐渐发生退化，且退化到一定程度时设备发生故障。

⑶周期性检查各个设备状态，检查时无故障发生，检查时间不计，各个设备的退化状态信息完全可通过检查得到。

⑷设备在生产过程存在产出次品的情况，且设备退化程度的加剧会导致产出损失逐渐增大。

⑸在同一时间点有且仅有一台设备需要预防性维修或完美更换，即不会有两台及以上的设备同时进行维修。

⑹设备在运行和维修、完美更换时都要消耗一定的能源，且进行维修和完美更换活动时的单位能耗大于运行状态下的单位能耗。

⑺系统中各个设备的额定工作速率为$p_i,i=\mathrm{1,2},...,n$，其中$M_l$为瓶颈，$l=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}\left\{p_i,i=\mathrm{1,2},...,n\right\}$，即额定加工速率最慢的设备。

⑻生产系统的维修周期是指系统从此次维修开始时刻到下一次维修开始时刻之间的时间间隔。

基于上述假设，首先对系统中设备的能耗、节能延迟时间窗、退化过程进行建模，接着根据最早达到预防性维修时间点的设备，计算出此次维修的时间窗组合后，以系统层目标函数最小化为目标，寻找最优的延迟窗决策点，由此构建随机退化维修模型。

2 随机退化延迟决策点模型构建

2.1　设备能耗

根据设备不同的工作状态，设备的单位时间能耗变化可分为几个阶段：启动、预热阶段、正常运行阶段、待机及停机阶段。设备在正常运行时的单位能耗为$e_i^{\mathrm{w}}$，那么预热时的单位能耗值为$g{e}_i^{\mathrm{w}}$，$g<\mathrm{1}$；当设备处于待机状态时的单位能耗值为$h{e}_i^{\mathrm{w}}$，此时$h<\mathrm{1}$或$h>\mathrm{1}$；停机状态下能耗为0。

2.2　设备连续退化过程及不完全维修建模

目前在描述生产系统退化轨迹的研究中，一般采用单随机过程建模，每种随机过程都有其特定的应用范围。在以往针对串行生产系统的维修建模中，大多研究不考虑能耗和产出损失，或直接假设能耗恒定，并未考虑生产系统中的能耗和次品率随着退化发生递增的情况，且都以单一指标作为维修决策依据，存在一定的局限性。

随着智能化高速发展，制造业中的产品结构也越来越多样化，本文在此背景下提出的串行生产系统考虑了能耗和产出损失都随设备退化而递增，并引入新型指标能耗效率（EEI），因此相较于传统串行生产系统，其退化机理更为复杂，采用单一随机过程建立的退化模型会存在准确性不佳、适用性不足等问题；利用模型选择准则从多个备选退化模型中确定合适模型的方法，容易出现模型误判且提升了工作的复杂度；而在现有模型上加入更多的影响因子，会显著增加计算与估计的难度，从而降低其易用性。针对上述不足，为了拓宽退化模型的应用范围，提高其可用性，基于指数散布型随机过程提出一类Tweedie指数散布退化模型（TEDM）^［14］。Tweedie指数散布过程（TEP）通过改变特定参数的值，可以将其转化为正态分布、伽马分布、逆高斯分布等常见的随机过程。由于TEDM模型似然函数结构的不确定性，基于鞍点近似法（SAM）极大似然估计方法给出了未知参数估计值。

2.2.1　Tweedie 指数散布过程

令$t$时刻的退化量为$X(t)$，若$X(t)~\mathrm{E}\mathrm{D}\mathrm{P}(\alpha t,\beta )$，其中$\alpha$称为偏移参数，$\beta$为散布参数，且两个参数之间互相独立，则随机过程$\left\{X(t)\left|t>\mathrm{0}\right.\right\}$称为指数散布型过程（exponential-dispersion process， EDP），且$X(t)$满足

（1）$X(\mathrm{0})=\mathrm{0}$；

（2）$\left\{X(t)\left|t>\mathrm{0}\right.\right\}$具有独立增量；

对任意的$\mathrm{\Delta }t\ge \mathrm{0}$，有$\mathrm{\Delta }X=X(t+\mathrm{\Delta }t)-X(t)~\mathrm{E}\mathrm{D}\mathrm{P}(\alpha (t+\mathrm{\Delta }t)-\alpha (t),\beta )$，且概率密度函数为

$$\begin{array}{l}f(\mathrm{\Delta }x\left|\alpha ,\beta \right.)=k(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\beta \right.)\\ e^{\left\{\beta \left[\mathrm{\Delta }x\cdot v(\alpha )-\mathrm{\Delta }t\cdot \varpi (v(\alpha ))\right]\right\}}\end{array}$$

（1）

式中：$k(\cdot )$是标准化函数；$\varpi (\cdot )$为适应度函数。

一般而言若一个指数散布型过程可用其方差函数定义，则称此指数散布型随机过程为Tweedie指数散布过程，其中方差函数可用式（2）表示

$$\mathrm{V}\mathrm{a}\mathrm{r}(\alpha )={\alpha }^{\chi },\chi \in \left(-\mathrm{\infty },\mathrm{0}\right]\bigcup \left[\mathrm{1},\mathrm{\infty }\right)$$

（2）

参数$\chi$不同时，其适应度函数$\varpi (\cdot )$的表达式可由其唯一确定。当$\chi$取不同值时，Tweedie指数散布过程可分别退化为常见的伽玛函数、逆高斯函数、正态分布以及复合泊松分布等。

在对Tweedie指数型散布随机过程（TEP）定义了解的基础上，考虑用TEP模型对本文中的生产系统进行退化过程建模，可以通过灵活变化参数$\chi$的值描述生产线上不同设备的退化过程。设初始时刻为全新状态，即初始退化量$X_i(\mathrm{0})=\mathrm{0}$，则在时刻$X_i(t)$退化量用TEP模型表示为

$$X_i(t)=X_i(\mathrm{0})+\mathrm{E}\mathrm{D}\mathrm{P}(\alpha t,\beta )=\mathrm{E}\mathrm{D}\mathrm{P}(\alpha t,\beta )$$

（3）

对于给定的失效阈值$D_i^{\mathrm{f}}$，相应的失效时间$L_i$的分布函数为

$$F_{{\alpha }_i(t),{\beta }_i}(t)=P(L_i<t)=P(X_i(t)>D_i^{\mathrm{f}})$$

（4）

TEP模型结构复杂，故设备首达失效时间分布难以确定，此处采用Birnbaum-Saunders分布近似设备首达时间分布，所以$L$关于模型参数$\mathit{\theta }=(\alpha ,\beta ,\chi )$的累积分布函数（CDF）为

$$F_L(t\left|\mathit{\theta }\right.)\cong \Phi \left[\sqrt[]{\frac{\beta }{{\alpha }^{\chi }}\left(\alpha \sqrt[]{t}-\frac{D^{\mathrm{f}}}{\sqrt[]{t}}\right)}\right]$$

（5）

式中：$\Phi (\cdot )$是标准正态分布的累积分布函数。设备的模型参数一般是未知的，可结合历史退化数据用极大似然估计方法得到。由于适应度函数$\varpi (\cdot )$会随着$\chi$波动而变化，所以模型中的显性表达式也无法确定，进一步参数的估计也无法继续。此处使用鞍点近似法先推导出概率密度函数（ PDF）的近似表达式，确定似然函数的表达式后，再用极大似然估计法得到模型参数$\mathit{\theta }=(\alpha ,\beta ,\chi )$的预估值。

2.2.2　TEDM模型未知参数的似然函数

模型未知参数为$\mathit{\theta }=(\alpha ,\beta ,\chi )$，假设共有$n$个样本，第$i$个样本在时间$t_{i,j}$时退化量为$X_{i,j}$，$\mathrm{1}\le i\le n$，$\mathrm{1}\le j\le G$，$G$为测量次数。设两次退化测量值之间的增量为$\mathrm{\Delta }{X}_{i,j}=X_{i,j}-X_{i,(j-\mathrm{1})},j\ge \mathrm{1}$；两次测量的时间增量可表示为$\mathrm{\Delta }{t}_{i,j}=t_{i,j}-t_{i,(j-\mathrm{1})},j\ge \mathrm{2}$，其中$\mathrm{\Delta }{t}_{i,\mathrm{1}}=t_{i,\mathrm{1}}$。对第$i$个样本，记$\mathrm{\Delta }{\mathit{X}}_i=(\mathrm{\Delta }{X}_{i,\mathrm{1}},\mathrm{\Delta }{X}_{i,\mathrm{2}},...,\mathrm{\Delta }{X}_{i,m}{)}^{\mathrm{T}}$，$f_{\mathrm{\Delta }{\mathit{X}}_i}(\mathrm{\Delta }{\mathit{x}}_i)$为退化量增量的联合概率密度函数。由于TEDM模型的增量具有独立平稳的性质，因此联合概率分布函数$f_{\mathrm{\Delta }{\mathit{X}}_i}(\mathrm{\Delta }{\mathit{x}}_i)$可表示为

$$f_{\mathrm{\Delta }{\mathit{X}}_i}(\mathrm{\Delta }{\mathit{x}}_i\left|\mathit{\theta }\right.)={\displaystyle \prod _{j=\mathrm{1}}^n}{f}_{\Delta X_{i,j}}\left(\Delta \right.x_{i,j}\left|\mathit{\theta }\right.)$$

（6）

考虑到样本之间相互独立，给定$n$个样本的退化数据$\mathrm{\Delta }\mathit{X}={\left\{\mathrm{\Delta }{\mathit{X}}_i\right\}}_i$，则关于TEDM未知参数$\mathit{\theta }$的似然函数为

$$\begin{array}{c}L(\mathit{\theta }\left|\mathrm{\Delta }\mathit{x}\right.)={\displaystyle \prod _{j=\mathrm{1}}^n}{f}_{\Delta {\mathit{X}}_{i,j}}\left(\Delta \right.{\mathit{x}}_{i,j}\left|\mathit{\theta }\right.)=\\ {\displaystyle \prod _{i=\mathrm{1}}^m}{\displaystyle \prod _{j=\mathrm{1}}^n}{f}_{\mathrm{\Delta }{X}_{i,j}}(\mathrm{\Delta }{x}_{i,j}\left|\mathit{\theta }\right.)\end{array}$$

（7）

当TEDM模型中的方差函数$\mathrm{V}\mathrm{a}\mathrm{r}(\alpha )={\alpha }^{\chi }$变化时，适应度函数$\varpi (\cdot )$的形式也随之变化，这极大地增加了似然函数$L(\mathit{\theta }\left|\mathrm{\Delta }\mathit{x}\right.)$的复杂度。一般的极值求解方法无法解决关于$L(\mathit{\theta }\left|\mathrm{\Delta }\mathit{x}\right.)$的极值优化问题，因此以下提出一种基于鞍点近似参数估计方法。

2.2.3　鞍点近似法（SAM）

鞍点近似法起源于复变函数，有效用于渐近分析。Barndorff-Nielson及Cox^［15］于1979年给出了基于极大似然估计法的密度函数的鞍点近似表达式，因此也被称为Barndorff-Nielson公式。鞍点近似法优势在于，即使样本量较小，也可以将原密度函数单调区间和趋势描述得很好，近似效果仍准确^［16］。所以当存在某个分布函数形式复杂并难以计算时，可以采用SAM得到密度函数的近似表达式。

该公式要求函数服从指数分布，且其分布函数的数学表达式满足一定的假设条件。所以，采用鞍点近似法来对TEDM模型中的退化量增量的概率密度函数进行简化，然后再利用极大似然近似法求出模型参数的估计值。

TEP模型的单位偏差可表示为

$$\begin{array}{c}d(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\alpha \right.)=\\ \mathrm{2}\left[\underset{\nu (\alpha )\in \mathrm{\Theta }}{\mathrm{s}\mathrm{u}\mathrm{p}}\left\{\nu (\alpha )\mathrm{\Delta }x/\mathrm{\Delta }t-\varpi (\nu (\alpha ))\right\}-\nu (\alpha )\mathrm{\Delta }x/\mathrm{\Delta }t+\varpi (\nu (\alpha ))\right]\end{array}$$

（8）

故真实退化增量$\mathrm{\Delta }X$的概率密度函数可写为

$$\begin{array}{l}{f}_{\mathrm{\Delta }X}(\mathrm{\Delta }x\left|\mathit{\theta }\right.)=k(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\beta \right.)\\ e^{\left\{\beta \mathrm{\Delta }t\underset{\nu (\alpha )\in \mathrm{\Theta }}{\mathrm{s}\mathrm{u}\mathrm{p}}\left[\mathrm{\Delta }x\nu (\alpha )/\mathrm{\Delta }t-\varpi (\nu (\alpha ))\right]\right\}}\\ e^{\left\{-\frac{\beta \mathrm{\Delta }t}{\mathrm{2}}d(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\alpha \right.)\right\}}\end{array}$$

（9）

其中，$\mathrm{\Theta }$为$\nu (\alpha )$的值域。

定理  　假设存在一个${\beta }_{\mathrm{0}}$，使得对于所有的$\forall \beta >{\beta }_{\mathrm{0}}$，概率密度函数式（9）都有界，则相应的鞍点近似表达式满足

$$\begin{array}{l}f(\mathrm{\Delta }x\left|\mathit{\theta }\right.)=\sqrt[]{\frac{\beta }{\mathrm{2}\mathrm{\pi }\mathrm{\Delta }{t}^{\mathrm{1}-\chi }\mathrm{\Delta }{x}^{\chi }}}{e}^{\left\{-\frac{\beta \mathrm{\Delta }t}{\mathrm{2}}d(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\alpha \right.)\right\}}\\ \left\{\mathrm{1}+O(\beta )\right\},\beta \mathrm{\Delta }t\to \mathrm{\infty }\end{array}$$

（10）

证明：式（9）的特征函数可表示为

$$\begin{array}{c}\varphi (z\left|\mathit{\theta }\right.)=E(e^{iz\mathrm{\Delta }x})=\\ \mathrm{e}\mathrm{x}\mathrm{p}\left\{\beta \mathrm{\Delta }t\left[\varpi (\nu (\alpha )+iz/\beta )-\varpi (\nu (\alpha ))\right]\right\}\end{array}$$

（11）

再由傅里叶反演定理，可得增量$\mathrm{\Delta }x$的概率密度函数为

$$f(\mathrm{\Delta }x\left|\mathit{\theta }\right.)=\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}{\int }_{-\infty }^{+\infty }\varphi (\left.z\right|\mathit{\theta })e^{-iz\mathrm{\Delta }x}\mathrm{d}z$$

（12）

因为式（12）中的被积函数是解析的，可将积分区域变为$i(\nu (\alpha )-\nu (\mathrm{\Delta }x/\mathrm{\Delta }t))+(-\mathrm{\infty },+\mathrm{\infty })$，则概率密度函数和单位偏差分别随之变为

$$\begin{array}{l}f(\mathrm{\Delta }x\left|\mathit{\theta }\right.)=\\ \frac{\beta }{\mathrm{2}\mathrm{\pi }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{e}^{\left[\beta \mathrm{\Delta }t\left\{\begin{array}{l}\varpi \left(\nu (\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t})+i\frac{z}{\beta }\right)-\\ \left(\nu (\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t})+i\frac{z}{\beta }\right)\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}+\nu (\alpha )\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}-\varpi (\nu (\alpha ))\end{array}\right\}\right]}\mathrm{d}\frac{z}{\beta }\end{array}$$

（13）

$$\begin{array}{c}d(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\alpha \right.)=\\ \mathrm{2}\left[\left\{\left.\left(\nu \left.\left(\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}\right.\right)\right.\right)\mathrm{\Delta }x/\mathrm{\Delta }t-\varpi \left(\nu \left(\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}\right)\right)\right\}-\nu (\alpha )\mathrm{\Delta }x/\mathrm{\Delta }t+\varpi (\nu (\alpha ))\right]\end{array}$$

（14）

将式（14）代入式（13）中的指数表达式中，并对$\nu (\mathrm{\Delta }x/\mathrm{\Delta }t)$用泰勒公式，得到式（15）

$$\begin{array}{l}\varpi \left(\nu \left(\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}\right)+i\frac{z}{\beta }\right)-\left(\nu \left(\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}\right)+i\frac{z}{\beta }\right)\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}+\\ \nu (\alpha )\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}-\varpi \left(\nu (\alpha )\right)=\\ \varpi \left(\nu \left(\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}\right)+i\frac{z}{\beta }\right)-\left(\nu \left(\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}\right)+i\frac{z}{\beta }\right)\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}+\\ \nu (\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t})\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}-\varpi \left(\nu (\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t})\right)-\\ \frac{\mathrm{1}}{\mathrm{2}}d(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\alpha ,\beta \right.)\approx \\ \frac{\mathrm{1}}{\mathrm{2}}\left(i{\left(\frac{z}{\beta }\right)}^{\mathrm{2}}\right){\varpi }^{"}\left(\nu \left(\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}\right)\right)-\frac{\mathrm{1}}{\mathrm{2}}d(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\beta \right.)\end{array}$$

（15）

再将式（15）代入式（13），概率密度函数变为

$$\sqrt[]{\frac{\beta }{\mathrm{2}\mathrm{\pi }\mathrm{\Delta }tVar\left(\mathrm{\Delta }x/\mathrm{\Delta }t\right)}}{e}^{\left[-\frac{\beta \mathrm{\Delta }t}{\mathrm{2}}d(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\alpha \right.)\right]}$$

（16）

其中，$Var\left(\mathrm{\Delta }x/\mathrm{\Delta }t\right)={(\mathrm{\Delta }x/\mathrm{\Delta }t)}^{\chi }$。在$\nu \left(\mathrm{\Delta }x/\mathrm{\Delta }t\right)$处，被积函数的指数有一个鞍点。将式（2）代入到式（16）中，得到

$$\begin{array}{l}f(\mathrm{\Delta }x\left|\mathit{\theta }\right.)=\sqrt[]{\frac{\beta }{\mathrm{2}\mathrm{\pi }\mathrm{\Delta }{t}^{\mathrm{1}-\chi }\mathrm{\Delta }{x}^{\chi }}}\\ e^{\left\{-\frac{\beta \mathrm{\Delta }t}{\mathrm{2}}d(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\alpha \right.)\right\}}\end{array}$$

（17）

其中，

$$\begin{array}{l}d(\mathrm{\Delta }x\left|\mathrm{\Delta }t,\alpha \right.)\\ \left\{\begin{array}{l}{\left(\mathrm{\Delta }x/\mathrm{\Delta }t-\alpha \right)}^{\mathrm{2}},\chi =\mathrm{0}(\mathrm{正态}\mathrm{分布})\\ \mathrm{2}\left\{\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}\mathrm{l}\mathrm{n}\left(\frac{\mathrm{\Delta }x}{\alpha \mathrm{\Delta }t}\right)-\left(\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}-\alpha \right)\right\},\chi =\mathrm{1}(\mathrm{泊松}\mathrm{分布})\\ \mathrm{2}\left\{\mathrm{l}\mathrm{n}\left(\frac{\alpha \mathrm{\Delta }t}{\mathrm{\Delta }x}\right)+\frac{\mathrm{\Delta }x}{\alpha \mathrm{\Delta }t}-\mathrm{1}\right\},\chi =\mathrm{2}(\mathrm{伽玛}\mathrm{分布})\\ \mathrm{2}\left[\frac{{\left\{\mathrm{m}\mathrm{a}\mathrm{x}\left(\mathrm{\Delta }x/\mathrm{\Delta }t,\mathrm{0}\right)\right\}}^{\mathrm{2}-\chi }}{(\mathrm{1}-\chi )(\mathrm{2}-\chi )}-\frac{{\alpha }^{\mathrm{1}-\chi }\mathrm{\Delta }x}{\left(\mathrm{1}-\chi \right)\mathrm{\Delta }t}+\frac{{\alpha }^{\mathrm{2}-\chi }}{\mathrm{2}-\chi }\right],\\ \chi \ne \mathrm{0,1},\mathrm{2}\end{array}\right.\end{array}$$

（18）

定理得证。

用鞍点近似法得到概率密度函数的近似表达式的基础上，可再用极大似然近似法求得模型未知参数$\alpha ,\beta ,\chi$的值。

对于给定$n$个样本，$G$次测量，第$i$个样本在时间$t_{ij}$时退化量为$X_{i,j}(\mathrm{1}\le i\le n,\mathrm{1}\le j\le G)$，$n$个样本的退化量测量数据$\mathrm{\Delta }\mathit{X}={\left\{\mathrm{\Delta }{\mathit{X}}_i\right\}}_i$，关于模型未知参数$\mathit{\theta }$的对数似然函数的近似表达式为

$$\begin{array}{l}\mathrm{𝓁}(\mathit{\theta }\left|\mathrm{\Delta }\mathit{x}\right.)\approx \mathrm{l}\mathrm{n}{\displaystyle \prod _{i=\mathrm{1}}^n}{\displaystyle \prod _{j=\mathrm{1}}^G}\sqrt[]{\frac{\beta }{\mathrm{2}\mathrm{\pi }\mathrm{\Delta }{t_{i,j}}^{\mathrm{1}-\chi }\mathrm{\Delta }{x_{i,j}}^{\chi }}}\\ e^{\left\{-\frac{\beta \mathrm{\Delta }{t}_{i,j}}{\mathrm{2}}d(\mathrm{\Delta }{x}_{i,j}\left|\mathrm{\Delta }{t}_{i,j},\alpha \right.)\right\}}\end{array}$$

（19）

对式（19）求偏导，并令偏导数为0，则得到$\widehat{\alpha }$的极大似然估计值为

$$\widehat{\alpha }={\sum }_{i=\mathrm{1}}^n{\sum }_{j=\mathrm{1}}^G\mathrm{\Delta }{x}_{i,j}/{\sum }_{i=\mathrm{1}}^n{\sum }_{j=\mathrm{1}}^G\mathrm{\Delta }{t}_{i,j}$$

（20）

将$\widehat{\alpha }$代入式（19）中就可以得到$\beta$和$\chi$的估计值，即$\widehat{\beta }$和$\widehat{\chi }$。

接下来考虑对生产线上各设备运用的预防性维修（PM）策略及维修效果。普遍来说，在真正的生产活动中，即使对达到维修阈值的设备进行了PM活动，其修复后效果往往并不能使其恢复到起始全新状态，而只是令其恢复到任一中间状态。

定义　  ^［17］    设$X$，$Y$为两个随机变量，假设对$\forall \alpha$，均有$P(X>\alpha )\ge P(Y>\alpha )$，则称变量$X$随机的大于变量$Y$，记作$X{\ge }_{\mathrm{s}\mathrm{t}}Y$。

定理：假设随着维修次数的增大，第$j+\mathrm{1}$次PM活动可使设备恢复到退化状态$X(M_{j+\mathrm{1}})$都随机的大于前一次的维修效果$X(M_j)$，即$X(M_{j+\mathrm{1}}){\ge }_{\mathrm{s}\mathrm{t}}X(M_j)$，那么第$j+\mathrm{1}$个维修周期的期望长度（维修间隔）比第$j$个维修周期的期望长度短。

其中，关于设备的预防性维修时间建模如下：假设设备$M_i$的预防性维修时间$w_i^{\mathrm{p}}$服从参数为${\lambda }_i$的指数分布，一般来说，设备的退化程度越大，则PM所需时间越长，这里假设${\lambda }_i\propto \mathrm{1}/X_i^{\mathrm{p}}$。于是，$w_i^{\mathrm{p}}$的概率密度函数为

$$g(w_i^{\mathrm{p}})={\lambda }_i{e}^{-{\lambda }_i{w}_i^{\mathrm{p}}}$$

（21）

$${\lambda }_i=d/X_i^{\mathrm{p}}$$

（22）

式（22）中：$d$为比例系数，$d>\mathrm{0}$。

本节利用Tweedie指数型散布随机过程对各台设备运行状态下的退化过程进行了建模，考虑PM修复效果恢复非新，且随着维修次数的增加，设备维修周期的期望间隔越来越短，每次PM时间和维修时间点的退化阈值呈正相关。

2.3　设备产出损失

由假设3可知，随着设备的退化，产品的质量也相应损失。本文用关于设备退化水平的次品率函数来衡量生产系统中由质量损失引起的相应损失为

$$p(X(t))=p_{\mathrm{0}}+\tau \left[\mathrm{1}-e^{(-\delta \cdot X{(t)}^{\mu }}\right]$$

（23）

式中：$p_{\mathrm{0}}$为设备全新时初始状态下的次品率，其中$\tau$为生产线上质量退化的边界值；$\delta ,\mu$均为正常数。

设备质量损失为

$$Q(X(t))=p(X(t))t_{i,j}^{\mathrm{w}}$$

（24）

式中：$t_{i,j}^{\mathrm{w}}$表示设备$i$在第$j$个维修周期的运作时间。

参考文献^［18］中关于质量损失与产出损失两者间函数关系，可得

$$W(X(t))=b_{\mathrm{1}}\left[\mathrm{1}-e^{-b_{\mathrm{2}}Q(X(t))}\right]$$

（25）

式中：$W(X(t))$代表由质量损失引起的产出损失；$b_{\mathrm{1}}$，$b_{\mathrm{2}}$为相关系数。

3 基于节能延迟时间窗口的预防性维修策略

3.1　设备能耗效率目标函数建模

基于节能延迟时间窗口这一概念及前文构建的状态退化模型，提出一种基于延迟时间窗口的机会维修策略，以便在降低生产系统能源消耗的同时尽可能提高系统的有效产出，设备的维修活动通常分为预防性维修、完美更换，这意味着维修周期$T_{i,j}$中能耗包括设备运行能耗$E_{i,j}^{\mathrm{w}}$、预防性维修能耗$E_{i,j}^{\mathrm{p}}$、完美更换能耗$E_{i,j}^{\mathrm{s}}$。一般说来，完美更换一台严重老化的机器比一台新机器要耗费更多的时间，完美更换的时间与设备已被维修过次数、退化量有关，所以，假设设备$i$在第$j$次更新时完美更换耗费时间$w_{i,j}^{\mathrm{s}}$为

$$w_{i,j}^{\mathrm{s}}=f\sqrt[]{X(t_{i,j})}/v^{j-\mathrm{1}}$$

（26）

式中：$f$，$v$为更换时间参数，且满足$f>\mathrm{1}$，$\mathrm{0}<v<\mathrm{1}$。$f$，$v$可通过已有数据进行函数拟合求出。

设备磨损是耗能的过程，而在其剧烈磨损时（需更换），可认为此时能耗等同于设备完美更换时的能耗，设备完美更换能耗与磨损关系可写为$E_{i,J}^{\mathrm{s}}=E_{i,j}^{\mathrm{w}}+c{W}^{\gamma }$。其中，$W$为磨损量。引入磨损量函数$W=o{w}_{i,j}^{\mathrm{s}}+\varphi e^{\phi t}$，其中$o\mathrm{、}\varphi \mathrm{、}\phi$都是常数。故有

$$E_{i,j}^{\mathrm{s}}=E_{i,j}^{\mathrm{w}}+c(o{w}_{i,j}^{\mathrm{s}}+\varphi e^{\phi t}{)}^{\gamma }$$

（27）

在此基础上，单台设备$i$在第$j$个维修周期中的能源消耗效率$\mathrm{E}\mathrm{E}{\mathrm{I}}_{i,j}$可表示为

$$\begin{array}{l}\mathrm{E}\mathrm{E}{\mathrm{I}}_{i,j}=\frac{\mathrm{设备}\mathrm{运行}\mathrm{能耗}+\mathrm{设备}\mathrm{预防}\mathrm{维护}\mathrm{能耗}+\mathrm{设备}\mathrm{完美}\mathrm{更换}\mathrm{能耗}}{\mathrm{设备}\mathrm{产出}-\mathrm{质量}\mathrm{损失}\mathrm{引起}\mathrm{的设}\mathrm{备产}\mathrm{出损}\mathrm{失}}=\\ \end{array}$$$$\frac{E_{i,j}^{\mathrm{w}}+E_{i,j}^{\mathrm{p}}+E_{i,j}^{\mathrm{s}}}{p_i{t}_{i,j}^{\mathrm{w}}-W(X(t_{i,j}))}$$

（28）

将式（25）~式（27）代入式（28）得

$$\mathrm{E}\mathrm{E}{\mathrm{I}}_{i,j}=\frac{\mathrm{2}{e}_i^{\mathrm{w}}{t}_{i,j}^{\mathrm{w}}+e_i^{\mathrm{p}}{w}_{i,j}^{\mathrm{p}}+c(o{w}_{i,j}^{\mathrm{s}}+\varphi e^{\phi t}{)}^{\gamma }}{p_i{t}_{i,j}^{\mathrm{w}}-b_{\mathrm{1}}\left[\mathrm{1}-e^{-b_{\mathrm{2}}H(X(t_{i,j}))}\right]}$$

（29）

3.2　系统最优维修策略决策建模

串行系统中的维修周期即系统从此次维修开始时刻到下次维修开始时刻这其中的时间间隔。串行系统规划期内的总能耗同样由运行能耗$E_{\mathrm{s}\mathrm{y}\mathrm{s}}^{\mathrm{w}}$、预防性维修能耗$E_{\mathrm{s}\mathrm{y}\mathrm{s}}^{\mathrm{p}}$和完美更换能耗$E_{\mathrm{s}\mathrm{y}\mathrm{s}}^{\mathrm{s}}$组成，其中系统中各个设备在规划期$T$内所经历的维修周期总次数为$k_i^{\mathrm{p}}$，本文的目标函数为最小化系统的能源效率指数$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}$，在上节构建的Tweedie指数型散布随机过程退化模型前提下，提出一种基于延迟时间窗口的机会维修策略，以达到在降低串行生产系统能源消耗的同时尽可能增大系统的有效产出。系统最优维修策略决策建模可表示为

$$\begin{array}{l}\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}=\frac{\mathrm{系统}\mathrm{运行}\mathrm{能耗}+\mathrm{系统}\mathrm{预防}\mathrm{维护}\mathrm{能耗}+\mathrm{系统}\mathrm{完美}\mathrm{更换}\mathrm{能耗}}{\mathrm{系统}\mathrm{产出}-\mathrm{质量}\mathrm{损失}\mathrm{引起}\mathrm{的系}\mathrm{统产}\mathrm{出损}\mathrm{失}}=\\ \frac{E_{\mathrm{s}\mathrm{y}\mathrm{s}}^{\mathrm{w}}+E_{\mathrm{s}\mathrm{y}\mathrm{s}}^{\mathrm{p}}+E_{\mathrm{s}\mathrm{y}\mathrm{s}}^{\mathrm{s}}}{p_{l}T-{\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^{k_i^{\mathrm{p}}}}{b}_{\mathrm{1}}\left[\mathrm{1}-e^{-b_{\mathrm{2}}Q(X_{i,j})}\right]}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^{k_i^{\mathrm{p}}}}\left.\left[(\mathrm{2}{e}_i^{\mathrm{w}}{t}_{i,j}^{\mathrm{w}}+e_i^{\mathrm{p}}{w}_{i,j}^{\mathrm{p}}+c(o{w}_{i,j}^{\mathrm{s}}+\varphi e^{\phi t}{)}^{\gamma }\right.\right]}{p_{l}T-{\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^{k_i^{\mathrm{p}}}}{b}_{\mathrm{1}}\left[\mathrm{1}-e^{-b_{\mathrm{2}}Q(X_{i,j})}\right]}\end{array}$$

（30）

以总能源消耗效率$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}$最小为最终目标函数进行求解，即可得到$\mathrm{\Delta }O$和$\mathrm{\Delta }o$的最优值，即系统的最优维修策略。

3.3　基于节能延迟时间窗口的维修决策

时间窗口维修策略的原理是以设备的随机停机事件作为多设备串联系统成组维修的机会，对设备维修活动进行有机集成，并动态调整原有的维修计划。由于整个系统会随着设备的维修活动停机，为减少生产系统的停机时间和停机损失，本文采用节能延迟时间窗（energy saving delay time window）方法完成维修活动的动态调度。

对于一个由$n$台设备组成的串联系统，一旦对系统中任一设备进行维修，整个系统都会停机。此时维修时间窗应起到整合的作用，串联系统的维修作业根据维修时间窗进行作业合并。传统的单窗策略中，停机时刻由系统中最早达到最优维修时刻的设备决定，其他需要机会维修的设备也在该停机时刻同时进行维修活动。而在节能延迟时间窗维修策略中，将这一停机时刻改为了机会维修决策触发时刻。

如图2所示，设$t_{i,j}$为设备$i$进行第$j$次维修的时刻，对应系统整体进行第$k$次维修所在时刻$t_k$，即$t_{i,j}=t_k$。分别用$\mathrm{\Delta }O$表示设备维修组合的时间窗，$\mathrm{\Delta }o$表示决策维修时机的时间窗。由节能延迟时间窗概念，该设备维修周期$j$对应系统的第$k$个维修阶段，因此，$t_{i,j}^{\mathrm{w}}$亦为设备$i$在系统第$k$个维修阶段内的最优维修时刻。即$t_{i,j}^{\mathrm{w}}=t_k$。此时可用时间窗$\mathrm{\Delta }O$确定机会维修设备组合。系统在第$k$个维修阶段内的机会维修范围所在区间为$\left[t_k^{\mathrm{*}},t_k^{\mathrm{*}}+\mathrm{\Delta }O\right]$。此时如果设备$i$的最优维修时刻$t_{i,j}^{\mathrm{w}}$满足$t_{i,j}^{\mathrm{w}}\le t_k^{\mathrm{*}}+\mathrm{\Delta }O$，则对设备$i$进行机会维修，否则不进行维修。

图2  节能延迟时间窗口

Fig2  Energy saving delay time window

下面考虑节能延迟时间窗的计算方法。决策过程如下：

（1） 系统初始化，从零时开始运行；

（2） 设设备$i$的预防性维修阈值$D_i^{\mathrm{p}}=\epsilon D_i^{\mathrm{f}}$，由式（9）可得到设备$i$达到预防性维修阈值时刻，即单设备$i$在第$j$个维修周期内的最优维修时刻$t_{i,j}^{\mathrm{w}}$；

（3） 以系统里当前最早到达最优维修时刻的设备$i$为准，其最优维修时刻作为现阶段的系统维修决策触发点，即$\mathrm{\Delta }O$的起始时刻，记为$t_k^{\mathrm{*}}$，令其他最优维修时刻落在区间$\left[t_k^{\mathrm{*}},t_k^{\mathrm{*}}+\mathrm{\Delta }O\right]$内的设备参与本次机会维护；

（4） 在确定系统第$k$个维修阶段内的维修设备组合后，进行该阶段维修时刻的决策。延迟时间窗策略下的维修时刻的决策触发时间点是$t_k^{\mathrm{*}}$，令其作为$\mathrm{\Delta }o$的起点时刻，使设备实际的维修时刻向后推迟$\mathrm{\Delta }o$，即系统在第$k$个维修阶段的最后维修决策时刻为$t_k=t_k^{\mathrm{*}}+\mathrm{\Delta }o$。通过对比不同$\mathrm{\Delta }O$和$\mathrm{\Delta }o$值下的系统能源效率指标值（$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}$），其中令$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}$最小的$\mathrm{\Delta }O$和$\mathrm{\Delta }o$的值为最优解$(\mathrm{\Delta }{O}^{\mathrm{*}},\mathrm{\Delta }{o}^{\mathrm{*}})$，进一步得到多设备生产系统机会维修最优维修时刻$t_k$。

接下来的目标就是要寻找最优解$(\mathrm{\Delta }{O}^{\mathrm{*}},\mathrm{\Delta }{o}^{\mathrm{*}})$，建立以单位时间内的系统能源效率指标值（energy efficiency indicator， EEI）最小化为目标的数学优化约束模型。

$$\begin{array}{l}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}(\mathrm{\Delta }O,\mathrm{\Delta }o)\\ \left\{\begin{array}{l}\mathrm{0}\le \mathrm{\Delta }O\le \mathrm{50},\mathrm{\Delta }O\in N^{+}\\ \mathrm{0}\le \mathrm{\Delta }o\le \mathrm{10},\mathrm{\Delta }o\in N^{+}\\ d>\mathrm{0}\\ f>\mathrm{1}\\ \mathrm{0}<v<\mathrm{1}\end{array}\right.\end{array}$$

（31）

4 算例分析

4.1　蒙特卡洛仿真算法

由于各台设备的预防性维修时间$t_{i,j}^{\mathrm{w}}$、维修时长$w_{i,j}^{\mathrm{p}}$及维修周期个数$k_i^{\mathrm{p}}$都是随机变量，很难用解析方法求解。随着系统预防性维修周期的变化，系统的实际产出和预防性维修能耗均会随之变化，本文目标是使单位能耗转化的系统实际产出最大。因此采用蒙特卡洛仿真方法对该串联生产系统进行建模，在假设各台设备退化服从Tweedie指数型散布随机过程的基础上，随机产生各设备达到维修阈值所经历的工作时间及相应的预防性维修时间，并通过最早达到最优维修时刻的设备来确定系统的机会维修组合窗起始点。根据蒙特卡洛抽样原理可知，当仿真时长足够长时，在特定时间窗阈值下的系统有效能源效率值将趋于真实值。因此，采用仿真算法进行迭代计算以寻求最优的$\mathrm{\Delta }O$和$\mathrm{\Delta }o$阈值，从而最小化整个系统的能源效率，如图3，具体步骤如下。

图3  基于节能延迟时间窗口的机会维修策略仿真流程

Fig.3  Simulation process of opportunistic maintenance based on energy saving delay time window

步骤1  　确定各个设备所服从TEP的参数${\alpha }_i$、${\beta }_i$和${\chi }_i$，更换时间参数$f$和$v$，磨损参数$o\mathrm{、}\varphi \mathrm{、}\phi$。参数可通过各设备历史故障数据拟合得到。并确定各设备正常工作时单位时间的能耗$e_i^{\mathrm{w}}$，单位时间预防性维修能耗$e_i^{\mathrm{p}}$。初始化各设备的运行时间、运行能耗、维修时间、维修能耗均为0。

步骤2  　设设备$i$的预防性维修阈值$D_i^{\mathrm{p}}=\epsilon D_i^{\mathrm{f}}$，由式（9）可得到各个设备达到预防性维修阈值时刻$t_{i,j}^{\mathrm{w}}$。

步骤3  　记$M_{\mathrm{r}}$为系统中最早到达最优维修时刻的设备，将$t_{\mathrm{r},j}^{\mathrm{w}}$作为当前阶段的系统机会维修组合窗决策点，即$\mathrm{\Delta }O$的起始时刻，记作$t_k^{\mathrm{*}}$，所有最优维修时刻满足$t_{i,j}^{\mathrm{w}}\le t_k^{\mathrm{*}}+\mathrm{\Delta }O$的设备即为参与此次机会维修的设备组合。

步骤4  　对系统内满足这次机会维修所有设备的组合窗$\mathrm{\Delta }O$确定后，更新参加此次维修活动的设备维修周期$j=j+\mathrm{1}$。

步骤5  　计算系统实际实施维修的时刻$t_k=t_k^{\mathrm{*}}+\mathrm{\Delta }o$。根据式（21）、式（22）和式（26）产生设备$M_i$进行预防性维修所需时间$w_i^{\mathrm{p}}$和完美更换时间$w_{i,j}^{\mathrm{s}}$。

步骤6  　分别记录$M_i$正常运行和预防性维修所消耗的能源：$E_{i,j}^{\mathrm{w}}=E_{i,j}^{\mathrm{w}}+t_{i,j}^{\mathrm{w}}{e}_i^{\mathrm{w}}$，$E_{i,j}^{\mathrm{p}}=E_{i,j}^{\mathrm{p}}+w_{i,j}^{\mathrm{p}}{e}_i^{\mathrm{p}}$。

步骤7  　更新系统维修阶段$k=k+\mathrm{1}$。

步骤8  　若$t_{{}_k}^{\mathrm{*}}>T$，根据式（30）计算出整个串联生产系统的能效指标并记录$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^{k_i^{\mathrm{p}}}}\left.\left[(\mathrm{2}{e}_i^{\mathrm{w}}{t}_{i,j}^{\mathrm{w}}+e_i^{\mathrm{p}}{w}_{i,j}^{\mathrm{p}}+c(o{w}_{i,j}^{\mathrm{s}}+\varphi e^{\phi t}{)}^{\gamma }\right.\right]}{p_{l}T-{\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^{k_i^{\mathrm{p}}}}{b}_{\mathrm{1}}\left[\mathrm{1}-e^{-b_{\mathrm{2}}Q(X_{i,j})}\right]}$，结束整个仿真过程，否则，返回步骤2。

4.2　基于仿真的遗传算法优化

遗传算法（genetic algorithm，GA）是根据大自然中生物体进化规律而设计提出的一类自适应随机迭代搜索算法，具有自组织、自学习、自适应的智能特征，寻优过程具有并行性，可有效避免陷入局部最优，目前已被广泛应用于机器学习、信号处理、自适应控制和人工生命、组合优化等领域优化问题的求解。遗传算法通过模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程来搜索最优解，利用计算机仿真运算，将问题的求解过程转换成类似生物进化中的染色体基因的自然选择、交叉、变异等过程，实现个体适应性的提高，使问题的解不断进化，以求得最优或近似最优解。在求解较为复杂的组合优化问题时，相对一些常规的优化算法，通常能够较快地获得较好的优化结果。

由于本文模型中含有较多带有不确定因素的随机变量，以及整数变量和连续变量，是一个随机混合整数优化问题，难以用解析方法求解，所以用设计启发式并嵌入蒙特卡洛仿真方法求解。具体的算法流程描述如下：

遗传算法中，把一个问题的可行解从其解空间转换到遗传算法所能处理的搜索空间的转换方法称为编码。目前的编码方法主要有二进制编码、十进制编码和实数编码。实数编码是将问题的解用实数表示。实数编码解决了二进制和十进制编码对算法精度和存储量的影响，同时便于优化中引入问题的相关信息。实数编码可以直观表达解的含义，极大减小搜索空间，因此运算速度提高，节约了运算时间。本文中优化问题的决策变量组合$(\mathrm{\Delta }O,\mathrm{\Delta }o)$，通过实数编码方案将其编码为染色体$\mathit{X}$，每个决策变量构成$\mathit{X}$的每个基因。维护的解用一串整数表示，如$m_{\mathrm{1}}$，$m_{\mathrm{2}}$，…$m_i$…，$m_n$，$n$代表需要维护的设备数量，这一整数串代表了设备维护的先后顺序，由于串联系统中每个设备生产速率和能耗速率不同，不同的维护顺序会导致系统能耗效率不同。

为保证种群的多样性，采用随机化的方式生成初始种群，$g_p$表示算法第$p$次迭代后的种群。

将系统能耗效率解$\mathit{X}=(x_{\mathrm{1}},x_{\mathrm{2}})$的适应度定义为$f(\mathit{X})=E(x_{\mathrm{1}},x_{\mathrm{2}}{)}^{-\mathrm{1}}$，其中$E(x_{\mathrm{1}},x_{\mathrm{2}})$由蒙特卡洛仿真算法计算得出。

选用轮盘赌法进行选择，步骤如下：

⑴计算每条染色体$\mathit{X}$的适应度值$f(\mathit{X})=E(x_{\mathrm{1}},x_{\mathrm{2}}{)}^{-\mathrm{1}}$；

⑵计算群体的适应度总和，设种群数为$M$：

$$F={\displaystyle \sum _{i=\mathrm{1}}^M}f({\mathit{X}}_i)$$

⑶计算每个染色体$\mathit{X}$的选择概率$P_i$：

$$P_i=f({\mathit{X}}_i)/F$$

⑷计算每个染色体$\mathit{X}$的累积概率$Q_i$：

$$Q_i={\displaystyle \sum _{k=\mathrm{1}}^i}{P}_k\mathrm{k}=\mathrm{1,2},...$$

⑸生成［0，1］之间的随机数$r$，如果$r>Q_{\mathrm{1}}$，就选择染色体${\mathit{X}}_{\mathrm{1}}$，否则选择第$i$个染色体${\mathit{X}}_i(\mathrm{2}\le i\le M)$，使得$Q_i\le r\le Q_{\mathrm{1}}$。

在本文优化中，采取次序交叉。首先随机产生两个交换点，定义两交换点之间的区域为匹配区，再对两个父代的匹配区进行交换。如图4，次序交叉的具体操作步骤如下：

图4  交叉示意图

Fig.4  Illustration of crossover

⑴随机产生两个交换点，交换中间的基因片段，并将此基因片段放入子代里，基因在染色体中的位置不变。

⑵在父代个体中删除从另一个父代个体交换过来的基因片段中的基因，将剩余的基因按顺序保存。

⑶将剩余基因按逆序方式，从子代第2交叉位置后开始填入。

在生物自然进化的过程中，染色体可能产生某些复制差错，变异产生新的染色体，表现出新的性状，产生新的个体。变异操作主要有互换变异、逆序变异、插入变异等，本文中的变异操作采用逆序变异，首先随机确定两个位置，然后将这两个位置之间的基因顺序颠倒。如图5，父代个体为［1 2 5 10 9 6 4 12 8 11 3］，逆序位置随机确定为2和6，则变异后的子代个体为［1 2 4 6 9 10 5 12 8 11 3］。

图5  逆序变异示意图

Fig.5  Illustration of inverse mutation

如图6，当算法总的迭代次数达到设定迭代次数时，遗传算法终止。

图6  基于仿真的优化算法流程

Fig.6  Optimization algorithm flow based on simulation

遗传算法运行参数中包含种群规模大小、最大迭代次数（算法终止条件）、交叉概率和变异概率等。合理选择参数能提高算法有效性，提升遗传算法性能，并增大运行效率，有效优化蒙特卡洛得到的近似结果。针对上文遗传算法中提到的各项相关参数，进行以下设置：

⑴交叉：过程中仅包含交叉概率一个参数。交叉概率决定着种群中进行交叉的频率，取值范围通常在0.5~1。交叉概率设置较大时，对种群的扰动效果也较大，可以使新产生的个体更好地融入群体中，但此时算法的收敛速度会变慢。另一方面，较小的交叉概率扰动不足，可能会使算法早熟，收敛到局部最优解。针对以上两个原因，将交叉概率设置为$p_c=\mathrm{0.7}$。

⑵变异：变异过程中包含变异概率$p_m$（决定变异的染色体条数）和变异率$m_u$（决定每条染色体变异的基因数）。变异概率取值范围一般为0.001~0.1。取值太大时易使算法退化为随机搜索算法，难以收敛，而过小则对增加种群多样性的效果不大。综合考虑本文研究模型的规模，将变异概率设为0.1，变异率设为0.15。

⑶种群大小和终止条件：种群大小应根据研究问题的规模决定。若种群较大，为保证种群多样性，但运算量也会增大，搜索过程变慢；而种群大小过小时不利于全局收敛。种群大小的取值范围一般为20~200。针对能耗效率生产系统模型，种群规模的大小初步设定为20。终止条件也与研究问题的规模有关。考虑到研究问题规模受限于时间和空间，遗传算法不可能无限长时间运行，因此，最大迭代次数设定为50。

4.3　数值仿真实验

利用蒙特卡洛仿真方法对生产系统进行建模，实例基于文献 ^［19］中的第四章节算法实例部分开发得来，以5台不同类型的生产设备（车床、钻机、塔式铣床、立式铣床、磨床）组成的批量生产系统为例，5台设备组成的串行连续制造生产线，$T=\mathrm{500}$，$o=\mathrm{0.2}$，$\varphi =\mathrm{0.5}$，$\phi =\mathrm{0.8}$，$\gamma =\mathrm{0.3}$，$\epsilon =\mathrm{20}\mathrm{\%}$，$f=\mathrm{2}$，$v=\mathrm{0.6}$。

TEP模型可以通过参数的变化而涵盖多种常见的随机过程，因此其在工程中的适用范围很广，TEP模型中系数$\chi$与几类常见随机过程模型的转换关系见表2。

表2  Tweedie指数型散布随机过程

Tab.2  Tweedie exponential-dispersion stochastic processes

$\chi$	退化量$X$服从的分布
$\chi <0$	极稳定分布
$\chi =0$	正态分布
$\chi =1$	泊松分布
$1<\chi <2$	复合泊松分布
$\chi =2$	伽马分布
$2<\chi <3$	逆高斯分布

由于各个设备退化过程是累积增长的渐变过程，具有非负独立增量性质，而伽马过程是一个独立增量服从相同尺度参数的连续随机过程。为此，选取伽马过程对设备的连续退化过程进行建模。令${\chi }_i=\mathrm{2}$，即设备的退化过程服从参数不同的伽马分布，各设备的参数见表3。

表3  串联系统中各设备参数

Tab.3  Production parameters of each machine in serial system

参数	$M_1$	$M_2$	$M_3$	$M_4$	$M_5$
$e_i^{\mathrm{w}}$	6	5	8	5	7
$e_i^{\mathrm{p}}$	12	15	16	13	14
${\alpha }_i$	5	2	2.7	1.05	1.25
${\beta }_i$	0.5	1.5	3.2	2.8	1.65
$p_i$	1.2	0.8	1.5	0.9	1.3
$D_i^{\mathrm{f}}$	150	130	145	135	155

最终，经遗传算法迭代后，得到最优解为$\mathrm{\Delta }{O}^{\mathrm{*}}=\mathrm{18}d,\mathrm{\Delta }{o}^{\mathrm{*}}=\mathrm{5}d$，系统能源效率指数收敛于$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}(\mathrm{18,5})=\mathrm{1.452}$，见图7。

图7  遗传算法迭代进化过程

Fig.7  The Evolution diagram of genetic algorithm

图8给出了延迟窗策略下串行生产系统中各个设备在维修周期计次下运行时间变化情况。

图8  设备运行时间

Fig.8  Operating time in devices

设备运行时间为0时代表设备在此维修周期中不参与机会维修。根据图8，我们可得系统在1~17个维修周期中停机时刻，见表4。

将预防性维修阈值从20%提高到60%，记录基于不同时间窗下的$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}$和维修集成窗最优值$\mathrm{\Delta }{O}^{\mathrm{*}}$绘制成图9。可以看出，无论哪种维修策略下的$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}$是都是随着PM阈值先增后减的，PM阈值较低时进入机会维修范围内设备较多，导致维修能耗的增加。PM阈值增大，设备到达PM阈值时刻所需时间越长，进行维修活动的次数减少，维修能耗因此降低。最优维修组合时间窗$\mathrm{\Delta }{O}^{\mathrm{*}}$在一定范围内，由于预防性维修阈值变大，设备运行时间、维修时间以及更新时间都会增加，需要通过增大维修集成时间窗$\mathrm{\Delta }{O}^{\mathrm{*}}$的值使各设备尽可能地同时维修，以减少系统的整体停机次数，从而减少由生产质量损失导致的产出损失。当维修阈值增加到一定数值时，此时所有设备都在集成窗维修范围内，$\mathrm{\Delta }{O}^{\mathrm{*}}$不受维修阈值变化影响，取值趋于稳定。

图9  不同时间窗下的$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}$和维修集成窗最优值$\mathrm{\Delta }{O}^{\mathrm{*}}$随PM 阈值的变化情况

Fig.9  $\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}$of two different time window and the optimal value of maintenance integration window$\mathrm{\Delta }{O}^{\mathrm{*}}$under different PM thresholds

一般来说，基于传统时间窗的预防性维修中，设备突发停机经常会使正常的生产安排紊乱，降低串行生产系统的有效产出。基于文献［19］中的第五章节算法实例部分开发得来，以7台不同类型的生产设备（车床、钻机、塔式铣床、钻机、塔式铣床、立式铣床、磨床）组成的批量生产系统为例，将设备数增加至8台，具体参数见表5并以传统单窗维修策略为对照，为验证延迟时间窗的有效性，控制其他参数及变量不变，传统时间窗口策略中的$t_k^{\mathrm{*}}$即为维修触发点。模拟结果显示，传统时间窗策略下，$\mathrm{\Delta }{O}^{\mathrm{*}}=\mathrm{27}d$时，$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}=\mathrm{2.751}$；延迟时间窗下，$\mathrm{\Delta }{O}^{\mathrm{*}}=\mathrm{21}d,\mathrm{\Delta }{o}^{\mathrm{*}}=\mathrm{4}d$时有$\mathrm{E}\mathrm{E}{\mathrm{I}}_{\mathrm{s}\mathrm{y}\mathrm{s}}=\mathrm{2.036}$。可知传统时间窗策略下的预防性维修活动能源效率更高，即延迟窗策略在维修时的能耗节约效率方面更优于单窗策略。

表5  增加设备后的参数设置

Tab.5  Parameters setting after increasing the number of machines

参数	$e_i^{\mathrm{w}}$	$e_i^{\mathrm{p}}$	${\alpha }_i$	${\beta }_i$	$p_i$	$D_i^{\mathrm{f}}$
$M_1$	6	12	5	0.5	1.2	150
$M_2$	5	15	2	1.5	0.8	130
$M_3$	8	16	2.7	3.2	1.5	145
$M_4$	5	13	1.05	2.8	0.9	135
$M_5$	7	14	1.25	1.65	1.3	155
$M_6$	5	14	3	2	1.2	152
$M_7$	4	11	2.8	2.5	1.4	148
$M_8$	6	13	3.6	1.45	1.1	149

5 结论

（1）以多设备串联生产系统为研究对象，在令系统中各设备退化服从Tweedie指数散布过程及不完全维修的基础上，考虑了生产过程中的产出损失，引入了“延迟时间窗”的机会维修策略，令传统时间窗中最早到达最优维修时刻决策点成为维修组合窗决策点，提出了一种基于延迟时间窗的延迟式预防性维修策略。

（2）通过最小化系统在规划期内的能源消耗效率为目标构建机会维修模型，利用蒙特卡洛仿真算法迭代求解出最佳的维修窗组合和预防性维修策略，并探究了不同的维修阈值对不同机会维修窗口下能耗效率指标和维修集成窗最优值的影响，结果表明在减少系统总能源消耗和提高系统有效产出方面，该机会维修模型切实可行。

（3）本文主要研究了串联生产系统中设备退化、能耗因素、质量损失与维护策略的优化问题。但在实际生产过程中，串并混联作为现代制造系统最普遍使用的配置型式，较串联系统而言，串并联系统加工过程中能耗情况、质量损失与维护决策要复杂得多，未来研究中可以考虑具体的质量和能耗特性指标在多阶段串并联生产系统中的传递性。

作者贡献声明

张文娟：提出研究选题，设计研究思路和论文框架；

朱雨彤：设计研究思路，完成实验设计并实施，分析数据，撰写论文。