CN102541032B - 一种可重入制造系统瓶颈设备预测方法 - Google Patents

Abstract

本发明公开一种基于ANFIS(自适应模糊推理系统)的可重入制造系统瓶颈设备预测方法。首先确定影响系统瓶颈设备的参数及ANFIS输入与输出；其次采用固定工件移动量间隔的方法采集原型系统数据并采用数据比对方法对数据进行预处理得到ANFIS的输入输出数据集；最后用训练数据确定ANFIS参数，用测试数据进行预测得到下一时刻瓶颈设备编号，设备利用率及队列长度并对预测结果进行分析得出预测精度，并且每次进行瓶颈预测之前，采用模式匹配方法对ANFIS模型参数在线调整。利用本发明能够构建可重入制造系统瓶颈设备预测方案，准确地预测出下一时刻系统的瓶颈设备及相关参数，并且对后续基于瓶颈设备的调度方法有指导意义。

Description

一种可重入制造系统瓶颈设备预测方法

技术领域

本发明属于先进制造领域。具体涉及大规模可重入制造系统中，一种瓶颈设备预测方法。

技术背景

可重入制造系统是以半导体制造系统为典型背景提出的。半导体制造业作为其关键性基础产业，已经成为当前的热点产业。如何采用性能优良的控制策略，提升半导体制造系统整体效率及各方面性能已经成为半导体制造业关注的焦点。在半导体制造系统中，瓶颈设备是制约系统产量、生产周期和在制品水平的关键因素，因此如何快速有效的识别生产线中瓶颈设备，并且以生产线中关键性的瓶颈设备控制为核心，带动其他生产设备的调度与管理，成为提高整个生产线性能指标的关键任务。因此快速准确预测出瓶颈设备，对提高调度算法的性能及整个生产线性能具有重要意义。

制造系统瓶颈主要指加工强度最大，负荷最高，对系统产出影响最大的加工中心。在对瓶颈研究中主要存在如下问题：(1)影响瓶颈设备的参数有很多，选取哪些参数作为计算瓶颈及影响瓶颈的关键因素，直接影响瓶颈设备的识别的准确性；(2)不管采用何种瓶颈计算方法，瓶颈的确定根据瓶颈计算值排序后选择的结果，但是生产过程中各种不确定的事件可能导致计算值及排序发生变化，从而造成瓶颈漂移，很多基于瓶颈的调度方法并没有考虑这种更为复杂的情况；(3)基于瓶颈的调度方法，都是在已知瓶颈的情况下采取合理的调度方法，但是这种调度方法存在一定的滞后性，可能生产线已经出现了累加，导致性能的下降，所以上述基于瓶颈的方法很难对瓶颈设备负载进行有效控制以及维持生产线的合理运行。

发明内容

本发明针对可重入制造系统中的瓶颈设备预测问题，给出一种瓶颈设备预测方法。该方法综合考虑生产线上影响瓶颈设备的关键因素，采用ANFIS(自适应模糊推理系统)对生产线下一时刻瓶颈设备及相关参数进行快速有效的预测。

1.一种可重入制造系统瓶颈设备预测方法，其特征在于，包含以下步骤：

(1)确定影响系统中瓶颈设备的关键因素为工件类型、工件投料方式、设备加工时间、设备故障间隔时间、设备平均维修时间；以及设备利用率、缓存区队列长度、生产线在制品水平WIP；

(2)获取以下实际生产线数据：

(2.1.1)每台设备的设备加工时间、设备故障间隔时间、设备平均维修时间；

(2.1.2)产品工艺路径：总共O道工序；

(2.1.3)工件类型、工件投料方式；

根据上述数据信息构建模型，模型由M台设备及工件投料装置、工件回收装置构成，投料装置按照某种投料方式进行投料，投放W种类型工件进入生产线中，工件按照工艺路径在设备中根据加工时间进行加工，最后由回收装置收集成品；

对该模型进行仿真，设工件在任意一台设备上完成一道工序的加工为一步，每完成一步进行一次数据采集，得到R₀组数据样本：{QL_s(1)，QL_s(2)，L，QL_s(t)，L，QL_s(R₀)}，{η_s(1)，η_s(2)，Lη_s(t)L，η_s(R₀)}，{W(1)，W(2)，L，W(t)，L，W(R₀)}，{Q(1)，Q(2)，L，Q(t)，L，Q(R₀)}，{T(1)，T(2)，L，T(t)，L，T(R₀)}，t＝1，2，L，R₀，s为设备编号，s＝1，2，L，M；其中，QL_s(t)、η_s(t)分别为第t步数据采集所得到的设备s的缓冲区队列长度与设备利用率，W(t)为第t步数据采集得到的生产线WIP，Q(t)、T(t)分别为第t步数据采集得到的生产线的投料方式及产品类型；

除({Q(1)，Q(2)，L，Q(t)，L，Q(R₀)}，{T(1)，T(2)，L，T(t)，L，T(R₀)}外，对每一个数据样本进行如下数据修正：

设采集得到的数据样本为{g(t)}，t＝1，2，L，R₀，且实际生产线中对应数据的样本均值为ξ，方差为σ，定义数据偏差为D(t)＝g(t)-ξ，若|D(t)|≤2σ，t＝1，2，L，R₀，则认为g(t)能够反映实际生产线状况，否则剔除第t步所采集的所有数据，得到修正后的数据；

(3)对上述修正过的数据按下式计算每个设备的瓶颈度，并进行数据再修正：

${\mathrm{BV}}_s\left(n\right)=\sqrt{{\mathrm{\θ}}_1{\mathrm{\η}}_s\left(n\right)+{\mathrm{\θ}}_2\frac{{\mathrm{QL}}_s\left(n\right)}{T}}---\left(1\right)$

式中QL_s(n)，η_s(n)，分别为第n步采集得到的设备s的缓冲区队列长度及设备利用率，BV_s(n)第n步计算得到的设备s瓶颈度，s为设备编号，n＝1，2，L，R₀，但其中不包括步骤(2)剔除的数据；T为缓冲数预估上限，θ₁与θ₂为经验参数，θ₁，θ₂∈[0，1]，θ₁+θ₂＝1；

设瓶颈度阈值B∈[0.8，0.95]，若第n步每台设备的瓶颈度BV_s(n)均小于B，则剔除第n步所采集数据，得到再修正的数据样本；若存在BV_s(n)≥B，取瓶颈度最大值所对应的设备作为系统中瓶颈设备，得到第n步瓶颈设备编号及相应的设备利用率、缓冲区队列长度数据；

(4)确定每台设备的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间，以及生产线中WIP，投料方式，工件类型为三个ANFIS的输入，下一时刻瓶颈设备编号、缓冲区队列长度、设备利用率分别为三个ANFIS的输出；

(5)采样数据经过数据预处理后，得到H+1组数据，按如下方式构建三个数据集 $[x_1^j,x_2^j,x_3^{j}L,x_{5\×M+1}^j,{\mathrm{pattern}}^j,{\mathrm{type}}^j,y_1^{j+1}],$ $[x_1^j,x_2^j,x_3^{j}L,x_{5\×M+1}^j,{\mathrm{pattern}}^j,{\mathrm{type}}^j,y_2^{j+1}],$ $[x_1^j,x_2^j,x_3^{j}L,x_{5\×M+1}^j,{\mathrm{pattern}}^j,{\mathrm{type}}^j,y_3^{j+1}],$ 及确定需要三个ANFIS；

式中j代表第j组数据，j＝1，2，L，H；

分别为第j组数据中设备1的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间；

分别为第j组数据中设备2的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间；以此类推；

第j组数据中生产线WIP，M为设备数量；

第j+1组数据中瓶颈设备编号；

第j+1组数据中瓶颈设备缓冲区队列长度；

第j+1组数据中瓶颈设备的设备利用率；

pattern^j：第j组数据中的投料方式；

type^j：第j组数据中的工件类型；

以上3个数据集均由H组数据组成，取每个数据集V组为训练数据，V＝[0.6H，0.7H]，剩余H-V组为测试数据，按如下步骤对数值型数据进行聚类，确定3个ANFIS结构，并对类别型数据进行编码，ANFIS结构确定包括每个输入变量对应隶属度函数及规则的个数；

(5.1)对每个数据集中的类别型数据按如下方式进行编码，形成一个行向量：

(5.1.1)设投料方式有A种，对第j组数据，令

若第j组数据中生产线采用某种投料方式，则对应的行向量中的元素为1，其余均为0；

(5.1.2)设工件类型有C种，对第j组数据，令

第j组数据中生产线存在几种类型工件，则对应的行向量中的元素为1，其余均为0；

(5.1.3)将pattern^j，type^j行向量构成向量s^j ，即 $s^j=[{\mathrm{pattern}}^j,{\mathrm{type}}^j]=[p_1^j,p_2^j,L,p_A^j,q_1^j,q_2^j,L,q_C^j],$ 则完成对类别数据的编码；

(5.2)对每个数据集的训练数据中数值型数据采用二叉树及模糊C均值进行聚类，并根据聚类后子空间形成初始模糊推理系统，步骤如下：

(5.2.1)令二叉树的根节点为LN(1)，LN(1)为训练数据中的数值型数据，LN(1)＝[X，Y₁]，X＝[x¹L x^jL x^V]^T，X为由x^j构成的矩阵，

x^j为由

构成的行向量，

为数据集中的第j组数据，分别代表每台设备的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间以及WIP，

Y₁为由构成的列向量，j＝1，2，L，V，

为第j+1组数据中的瓶颈设备编号；初始化叶节点集合为Ψ＝{LN(M)}，整棵树的节点M＝1；

(5.2.2)设定误差限R，R∈[0.001，0.05]，对叶节点集合Ψ中所有节点采用二叉树法进行判定，集合中每个节点是否需要被继续划分，若节点M的线性拟合误差e(M)大于R，则按照步骤(5.2.3)继续划分，若所有节点的线性拟合误差均小于R，则转到步骤(5.2.4)，；

(5.2.3)采用模糊C均值对e(M)大于误差限R的节点LN(M)进行聚类，详细聚类步骤如下：

(5.2.3.1)给定聚类中心数为c＝2，设定容许误差E_max及加权指数m＝2，E_max∈[0.001，0.05]，随机初始化聚类中心；

(5.2.3.2)若模糊C均值的聚类准则函数值J_m≤E_max，则结束聚类；

(5.2.3.3)节点LN(M)被划分成两个叶节点LN(M+1)，LN(M+2)；

(5.2.3.4)更新节点集合：将Ψ与{LN(M+1)，LN(M+2)}合并，更新Ψ；

(5.2.3.5)返回第(5.2.2)步；

(5.2.4)得到叶节点数据集Ψ＝{LN(1)，LN(2)，L，LN(K₁)}，LN(1)，LN(2)，L，LN(K₁)为聚类后形成的数据子空间，K₁为通过上述聚类过程得到的聚类中心的个数；

(5.2.5)根据聚类中心个数确定模糊系统的模糊规则数为K₁；

将(5.2.1)中的 $Y_1={\left[y_1^{2}L{y}_1^{j+1}L{y}_1^{V+1}\right]}^T,$ 分别替换成 $Y_2={\left[y_2^{2}L{y}_2^{j+1}L{y}_2^{V+1}\right]}^T,$ $Y_3={\left[y_3^{2}L{y}_3^{j+1}L{y}_3^{V+1}\right]}^T,$ Y₁为由

构成的列向量，Y₃为由构成的列向量，j＝1，2，L，V，

为第j+1组数据中的瓶颈设备利用率，为第j+1组数据中的瓶颈设备缓冲区队列长度，重复上述聚类过程，得到每个ANFIS的中模糊规则数K₁，K₂，K₃；

(6)按照如下步骤确定三个ANFIS中的参数：

(6.1)初始化参数集

其中

称为前件参数，i＝1，2，L，5×M+1，k＝1，2，LK₁，K₁为模糊规则数量，v_i，k，σ_i，k均为(0，1)中的随机数；

为(A+C)×K₁维矩阵，矩阵中每个元素均为(0，1)中的随机数；

称为后件参数，k＝1，2，L K₁，均为(0，1)中的随机数；

(6.2)第一个数据集中训练数据输入到ANFIS中，确定第一个ANFIS中参数：

(6.2.1)对数值型数据从ANFIS第一层输入，第一层有5×M+1个输入量，输入第j组数据为

j＝1，2，L，V，

为数据集中的第j组数据，分别代表每台设备的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间以及WIP，第一层对输入变量模糊化，输出为对应模糊集的隶属度，这里隶属度函数采用高斯钟型函数：

$O_{1,h}^j={\mathrm{\μ}}_{A,i}^k\left(x_i^j\right)=\exp (-\frac{1}{2}{\left(\frac{x_i^j-v_{i,k}}{{\mathrm{\σ}}_{i,k}}\right)}^2)---\left(2\right)$

其中i＝1，2，L，5×M+1；k＝1，2，L K₁；h＝(5×M+1)×K₁，K₁为模糊规则数量，

为隶属度函数，

代表第一层输出，其中

为前件参数；

(6.2.2)第二层，采用乘法规则计算每条模糊规则的激活强度w_k，第二层各节点的输出为

k＝1，2L，K₁；

$O_{2,k}^j=w_k=\underset{i=1}{\overset{26}{\mathrm{\Π}}}{\mathrm{\μ}}_{A,i}^k\left(x_i^j\right)---\left(3\right)$

(6.2.3)第三层，归一化得到各条模糊规则的激活强度

输出为

k＝1，2L，K₁；

$O_{3,k}^{\′j}=\stackrel{\‾}{w_k}=\frac{w_k}{\underset{k=1}{\overset{K_1}{\mathrm{\Σ}}}{w}_k}---\left(4\right)$

对编码后的类别输入经过一个激活强度转换矩阵T₁，且

输入到第三层F^j＝s^jT₁

F^j＝s^jT₁                           (5)

其中j＝1，2，L，V， $s^j=[{\mathrm{pattern}}^j,{\mathrm{type}}^j]=[p_1^j,p_2^j,L,p_A^j,q_1^j,q_2^j,L,q_C^j]$ 为1×(A+C)行向量，

为K₁维行向量，T₁为(A+C)×K₁维矩阵；

则第三层的总输出为；

k＝1，2L，K₁，j＝1，2，L，V；

(6.2.4)第四层，每个节点的传递函数为线性函数，计算得到每条规则的输出O_4，k

$O_{4,k}^j=O_{3,k}^j(d_0^k+d_1^k{x}_1^j+d_2^k{x}_2^j+L+d_{5\×M+1}^k{x}_{5\×M+1}^j)=O_{3,k}^j\·g_k\left(x^j\right)---\left(6\right)$

式中 $g_k\left(x^j\right)=(d_0^k+d_1^k{x}_1^j+d_2^k{x}_2^j+L+d_{5\×M+1}^k{x}_{5\×M+1}^j),$ $P_1^c=\{d_0^k,d_1^k,d_2^k,L,d_{5\×M+1}^k\}$ 为后件参数， $x^j=[x_1^j,L,x_{5\×M+1}^j],$ ，k＝1，2，L，K₁，j＝1，2，L，V；

(6.2.5)计算所有规则输出之和

$O_5^j=\underset{k=1}{\overset{K_1}{\mathrm{\Σ}}}{O}_{4,k}^j---\left(7\right)$

(6.3)设定最大训练次数epoc及误差限error，epoc∈[100，500]，error∈[0，0.005]，训练数据的总输出误差为：

$E_0=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{e}_t=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\frac{1}{2}{[O_5^j-y_1^{j+1}]}^2---\left(8\right)$

式中

为Y₁中第j+1组数据中的瓶颈设备编号，J＝V，Y₁为第一个ANFIS中训练数据中的输出数据，

为第j组训练数据输入到ANFIS中得到的输出；

若E₀≤error，或者达到最大训练次数epoc，结束训练并转到(6.5)，否则进行(6.4)

(6.4)训练参数集

(6.4.1)固定参数集

采用误差反向传播算法对进行修正；

(6.4.2)固定参数集

采用最小二乘法对

进行修正；

(6.4.3)固定参数集

采用最小二乘法对进行修正；

(6.4.4)得到修正后的参数集

返回步骤(6.2)；

(6.5)得到确定参数集

用第二个数据集的训练数据对第二个ANFIS进行训练，重复上述训练过程，保持ANFIS的输入数据不变，输出数据将替换Y₁为Y₂，

J＝V，

为Y₂中第j+1组数据中的瓶颈设备利用率，K₁将替换为K₂，确定第二个ANFIS的模糊规则数为K₂，训练ANFIS，确定第二个ANFIS的参数集

用第三个数据集的训练数据对第三个ANFIS进行训练，重复上述训练过程，保持ANFIS的输入数据不变，输出数据将替换Y₁为Y₃，

J＝V，

为Y₃中第j+1组数据中的瓶颈设备缓冲区队列长度，K₁将替换为K₃，确定第二个ANFIS的模糊规则数为K₃，训练ANFIS，确定第三个ANFIS的参数集

按照以上步骤用训练数据完成对三个ANFIS的训练，训练完成后；将后H-V组测试数据输入到ANFIS中，得到预测输出。

还可以在步骤(6)后包括以下步骤：

(7)将生产线所获得的实时数据，与原始数据进行分析对比，若其值低于原始数据值区间上限的120％或高于区间下限的80％，则将数据直接输入到ANFIS，得到下一时刻预测输出值；若数据值超出原始数据值区间上限的120％或高于区间下限的80％，则认为该组数据是在新的工况下获得，将该数据与原始训练数据组合重新作为ANFIS的训练数据集，返回到步骤(5)，对ANFIS中的参数进行在线调整，得到下一时刻瓶颈设备相关参数的预测值。

本发明一种可重入制造系统瓶颈设备的预测方法，与现有技术相比，具有以下明显的优势和有益效果：设计目标清晰，过程明确，现实可行。本发明能够构建可重入制造系统瓶颈设备预测方案，准确地预测出下一时刻系统的瓶颈设备及相关参数，并且对后续基于瓶颈设备的调度方法有指导意义。

附图说明

图1是本发明中瓶颈设备预测框图；

图2是本发明中瓶颈设备预测实施方法流程图；

图3是本发明中ANFIS示意图；

图4是本发明中瓶颈设备编号的训练输出与实际输出对比曲线；

图5是本发明中瓶颈设备编号的预测输出与实际输出对比曲线；

图6是本发明中瓶颈设备利用率的训练输出与实际输出对比曲线；

图7是本发明中瓶颈设备利用率的预测输出与实际输出对比曲线；

图8是本发明中缓冲区队列长度的训练输出与实际输出对比曲线；

图9是本发明中缓冲区队列长度的预测输出与实际输出对比曲线。

具体实施方式

为了更好地理解本发明的技术方案，在此以一个半导体生产线Minifab为例进行分析。Minifab生产3种类型工件，由5台设备M_a，M_b，M_c，M_d，M_e构成(M_a，M_b)，(M_c，M_d)，(M_e)3个设备群，其中M_a，M_b为批加工设备，一次可加工3卡硅片，M_c，M_d，M_e为单件加工设备，一次只能加工一卡硅片，工件的工艺流程共包含六道工序。

(1)确定以下因素为影响瓶颈设备的关键因素：

(1.1)工件类型：A₁，A₂，A₃共3种；

(1.2)工件投料方式：固定时间投料，固定在制品投料，共2种；

(1.3)设备加工时间：M_a加工时间为230分钟，M_b加工时间为240分钟，M_c加工时间为50分钟，M_d加工时间为40分钟，M_e加工时间为40分钟；

(1.4)设备故障间隔时间：M_a间隔时间为1400分钟，M_b间隔时间为1500分钟，M_c间隔时间为700分钟，M_d间隔时间为700分钟，M_e间隔时间为800分钟；

(1.5)设备平均维修时间：M_a维修时间为40分钟，M_b维修时间为70分钟，M_c M_c维修时间为120分钟，M_d维修时间为100分钟，M_e维修时间为30分钟；

(1.6)设备利用率：M_a、M_b利用率∈[0.85，095]，M_c、M_d利用率∈[0.40，0.55]，M_e利用率∈[0.9，1]；

(1.7)缓存区队列长度：取值范围[0，15]，即缓冲区最多容纳15个工件；

(1.8)WIP(生产线在制品水平)：取值范围[25，35]；

(2)获取实际生产线数据：

(2.1.1)设备加工信息：

设备	平均加工时间/分钟	故障间隔时间/分钟	维修时间/分钟
				Ma	230	1400	40
Mb	240	1500	70
				Mc	50	700	120
Md	40	700	100
				Me	40	800	30

其中M_a，M_b为批加工设备，最大加工批量3卡；

(2.1.2)产品工艺路径：ENTER-M_a，M_b-M_c，M_d-M_e-M_c，M_d-M_a，M_b-M_e-EXIT；总共六道工序；

(2.1.3)工件信息：A₁，A₂，A₃共3种类型工件，工件需求量不限；

(2.1.4)投料方式：选择固定时间投料或固定在制品投料。

根据上述数据信息构建MiniFab模型，模型由5台设备及工件投料装置、工件回收装置构成，投料装置按照固定在制品或固定时间进行投料，投放3种类型工件进入生产线中，工件按照工艺路径在设备中根据加工时间进行加工，最后由回收装置收集成品。

对该模型进行仿真，仿真时间以每天24小时不间断仿真，共仿真180天，并去掉系统的warm-up时间(30天)。设工件在任意一台设备上完成一道工序的加工为一步，例如工件在设备M_a上完成第一道工序的加工称为一步，每完成一步进行一次数据采集，得到数据样本：{QL_s(1)，QL_s(2)，L，QL_s(t)，L，QL_s(2000)}，{η_s(1)，η_s(2)，L，η_s(t)，L，η_s(2000)}，{W(1)，W(2)，L，W(t)，L，W(2000)}，{Q(1)，Q(2)，L，Q(t)，L，Q(2000)}，{T(1)，T(2)，L，T(t)，L，T(2000)}，t＝1，2，L，2000，s＝M_a，M_b，M_c，M_d，M_e；其中，QL_s(t)，η_s(t)分别为第t步数据采集所得到的设备s的缓冲区队列长度与设备利用率，W(t)为第t步数据采集得到的生产线WIP，Q(t)，T(t)分别为第t步数据采集得到的生产线的投料方式及产品类型

对每一个数据样本({Q(1)，Q(2)，L，Q(t)，L，Q(2000)}，{T(1)，T(2)，L，T(t)，L，T(2000)}除外)进行如下数据修正：

设采集得到的数据样本为{g(t)}，t＝1，2，L，2000，且实际生产线中对应数据的样本均值为ξ，方差为σ，定义数据偏差为D(t)＝g(t)-ξ，若|D(t)|≤2σ，t＝1，2，L，2000，则认为g(t)能够反映实际生产线状况，否则剔除第t步所采集的所有数据，得到修正后的数据。

(3)对上述修正过的数据按下式计算每个设备的瓶颈度，并进行数据再修正：

${\mathrm{BV}}_s\left(n\right)=\sqrt{{\mathrm{\θ}}_1{\mathrm{\η}}_s\left(n\right)+{\mathrm{\θ}}_2\frac{{\mathrm{QL}}_s\left(n\right)}{T}}---\left(1\right)$

式中QL_s(n)，η_s(n)，分别为第n步采集得到的设备s的缓冲区队列长度及设备利用率，BV_s(n)第n步计算得到的设备s瓶颈度，s＝M_a，M_b，M_c，M_d，M_e，n＝1，2，L，2000，T为缓冲数预估上限，取T＝15，θ₁与θ₂为经验参数，这里取θ₁＝0.4，θ₂＝0.6。

对设备M_a，M_b，M_c，M_d，M_e分别编号为1，2，3，4，5；设瓶颈度阈值B＝0.85，若第n步每台设备的瓶颈度BV_s(n)均小于B，则剔除第n步所采集数据，得到再修正的数据样本；若存在BV_s(n)≥B，取瓶颈度最大值所对应的设备作为系统中瓶颈设备，得到第n步瓶颈设备编号及相应的设备利用率、缓冲区队列长度数据。

(4)确定每台设备的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间，以及生产线中WIP，投料方式，工件类型为三个ANFIS的输入，下一时刻瓶颈设备编号、缓冲区队列长度、设备利用率分别为三个ANFIS的输出。

(5)采样数据经过数据预处理后，得到1301组数据，按如下方式构建三个数据集 $[x_1^j,x_2^j,x_3^{j}L,x_{5a+b}^j,L{x}_{26}^j,{\mathrm{pattern}}^j,{\mathrm{type}}^j,y_1^{j+1}],$ $[x_1^j,x_2^j,x_3^{j}L,x_{5a+b}^j,L{x}_{26}^j,{\mathrm{pattern}}^j,{\mathrm{type}}^j,y_2^{j+1}],$ $[x_1^j,x_2^j,x_3^{j}L,x_{5a+b}^j,L{x}_{26}^j,{\mathrm{pattern}}^j,{\mathrm{type}}^j,y_3^{j+1}],$ 及确定需要三个ANFIS。

式中j代表第j组数据，j＝1，2，L，1300；

分别为第j组数据中设备M_a的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间；

分别为第j组数据中设备M_b的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间；

分别为第j组数据中设备M_c的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间；

分别为第j组数据中设备M_d的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间；

分别为第j组数据中设备M_e的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间；

第j组数据中生产线WIP；

第j+1组数据中瓶颈设备编号；

第j+1组数据中瓶颈设备缓冲区队列长度；

第j+1组数据中瓶颈设备的设备利用率；

pattern^j：第j组数据中的投料方式；

type^j：第j组数据中的工件类型；

以上3个数据集均由1300组数据组成，取每个数据集前800组为训练数据，后500组为测试数据，按如下步骤对数值型数据进行聚类，确定3个ANFIS结构，ANFIS结构确定包括每个输入变量对应隶属度函数及规则的个数，并对类别型数据进行编码：

(5.1)对每个数据集中的类别型数据按如下方式进行编码，形成一个行向量：

(5.1.1)设投料方式有固定时间投料，固定在制品投料，共2种，对第j组数据，令

若第j组数据中生产线投料方式为固定时间投料，则

若第j组数据中生产线投料方式为固定在制品投料，则

(5.1.2)设工件类型有A₁，A₂，A₃，共3种，对第j组数据，令

第j组数据中生产线存在几种类型工件则对应的行向量中的元素为1，其余均为0；

(5.1.3)将pattern^j，type^j行向量构成向量s^j，即则完成对类别数据的编码；

(5.2)对每个数据集的训练数据中数值型数据采用二叉树及模糊C均值进行聚类，并根据聚类后子空间形成初始模糊推理系统，步骤如下：

(5.2.1)令二叉树的根节点为LN(1)，LN(1)为800组训练数据中的数值型数据，LN(1)＝[X，Y₁]，X＝[x¹L x^jL x^J]^T，X为由x^j构成的矩阵，

x^j为由

构成的行向量，

为数据集中的第j组数据，分别代表每台设备的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间以及WIP，

Y₁为由

构成的列向量，j＝1，2，L，J，J＝800，

为第j+1组数据中的瓶颈设备编号；初始化叶节点集合为Ψ＝{LN(M)}，整棵树的节点M＝1；

(5.2.2)设定误差限R＝0.01，对叶节点集合Ψ中所有节点采用二叉树法进行判定，集合中每个节点是否需要被继续划分，若节点M的线性拟合误差e(M)大于R，则按照步骤(5.2.3)继续划分，若所有节点的线性拟合误差均小于R，则转到步骤(5.2.4)，；

(5.2.3)采用模糊C均值对e(M)大于误差限R的节点LN(M)进行聚类，详细聚类步骤如下：

(5.2.3.1)给定聚类中心数为c＝2，设定容许误差E_max＝0.01及加权指数m＝2，随机初始化聚类中心；

(5.2.3.2)若模糊C均值的聚类准则函数值J_m≤E_max，则结束聚类；

(5.2.3.3)节点LN(M)被划分成两个叶节点LN(M+1)，LN(M+2)；

(5.2.3.4)更新节点集合：将Ψ与{LN(M+1)，LN(M+2)}合并，更新Ψ；

(5.2.3.5)返回第(5.2.2)步；

(5.2.4)得到叶节点数据集Ψ＝{LN(1)，LN(2)，L，LN(K₁)}，LN(1)，LN(2)，L，LN(K₁)为聚类后形成的数据子空间，K₁为通过上述聚类过程得到的聚类中心的个数；

(5.2.5)根据聚类中心个数确定模糊系统的模糊规则数为K₁；

将(5.2.1)中的 $Y_1={\left[y_1^{2}L{y}_1^{j+1}L{y}_1^{J+1}\right]}^T,$ 分别替换成 $Y_2={\left[y_2^{2}L{y}_2^{j+1}L{y}_2^{J+1}\right]}^T,$ $Y_3={\left[y_3^{2}L{y}_3^{j+1}L{y}_3^{J+1}\right]}^T,$ Y₁为由构成的列向量，Y₃为由构成的列向量，j＝1，2，L，J，J＝800，

为第j+1组数据中的瓶颈设备利用率，

为第j+1组数据中的瓶颈设备缓冲区队列长度，重复上述聚类过程，得到每个ANFIS的中模糊规则数K₁，K₂，K₃分别为30，16，12。

(6)按照如下步骤确定三个ANFIS中的参数：

(6.1)初始化参数集

其中

称为前件参数，i＝1，2，L，26；k＝1，2，L K₁，K₁为模糊规则数量，v_i，k，σ_i，k均为(0，1)中的随机数；

为5×K₁维矩阵，矩阵中每个元素均为(0，1)中的随机数；

称为后件参数，k＝1，2，L K₁，

均为(0，1)中的随机数。

(6.2)第一个数据集中训练数据输入到ANFIS中，确定第一个ANFIS中参数：

(6.2.1)对数值型数据从ANFIS第一层输入，第一层有26个输入量，输入第j组数据为

j＝1，2，L，800，

为数据集中的第j组数据，分别代表每台设备的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间以及WIP，第一层对输入变量模糊化，输出为对应模糊集的隶属度，这里隶属度函数采用高斯型函数：

$O_{1,h}^j={\mathrm{\μ}}_{A,i}^k\left(x_i^j\right)=\exp (-\frac{1}{2}{\left(\frac{x_i^j-v_{i,k}}{{\mathrm{\σ}}_{i,k}}\right)}^2)---\left(2\right)$

其中i＝1，2，L，26；k＝1，2，L K₁；h＝26×K₁，K₁为模糊规则数量，

为隶属度函数，

代表第一层输出，其中

为前件参数；

(6.2.2)第二层，采用乘法规则计算每条模糊规则的激活强度w_k，第二层各节点的输出为

k＝1，2L，K₁；

$O_{2,k}^j=w_k=\underset{i=1}{\overset{26}{\mathrm{\Π}}}{\mathrm{\μ}}_{A,i}^k\left(x_i^j\right)---\left(3\right)$

(6.2.3)第三层，归一化得到各条模糊规则的激活强度输出为

k＝1，2L，K₁；

$O_{3,k}^{\′j}=\stackrel{\‾}{w_k}=\frac{w_k}{\underset{k=1}{\overset{K_1}{\mathrm{\Σ}}}{w}_k}---\left(4\right)$

对编码后的类别输入经过一个激活强度转换矩阵T₁，且

输入到第三层F^j＝s^jT₁

F^j＝s^jT₁                    (5)

其中j＝1，2，L，800， $s^j=[{\mathrm{pattern}}^j,{\mathrm{type}}^j]=[p_1^j,p_2^j,q_1^j,q_2^j,q_3^j]$ 为1×5行向量， $F^j=[f_1^j,f_2^j,L,f_{K_1}^j]$ 为K₁维行向量，T₁为5×K₁维矩阵；

则第三层的总输出为；

k＝1，2L，K₁，j＝1，2，L，800；

(6.2.4)第四层，每个节点的传递函数为线性函数，计算得到每条规则的输出O_4，k

$O_{4,k}^j=O_{3,k}^j(d_0^k+d_1^k{x}_1^j+d_2^k{x}_2^j+L+d_{26}^k{x}_{26}^j)=O_{3,k}^j\·g_k\left(x^j\right)---\left(6\right)$

式中 $g_k\left(x^j\right)=(d_0^k+d_1^k{x}_1^j+d_2^k{x}_2^j+L+d_{26}^k{x}_{26}^j),$ $P_1^c=\{d_0^k,d_1^k,d_2^k,L,d_{26}^k\}$ 为后件参数， $x^j=[x_1^j,L,x_{26}^j],$ ，k＝1，2，L，K₁，j＝1，2，L，800；

(6.2.5)计算所有规则输出之和

$O_5^j=\underset{k=1}{\overset{K_1}{\mathrm{\Σ}}}{O}_{4,k}^j---\left(7\right)$

(6.3)设定最大训练次数epoc＝100及误差限error＝0.00001，训练数据的总输出误差为：

$E_0=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{e}_t=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\frac{1}{2}{[O_5^j-y_1^{j+1}]}^2---\left(8\right)$

式中

为Y₁中第j+1组数据中的瓶颈设备编号，J＝800，Y₁为第一个ANFIS中训练数据中的输出数据，

为第j组训练数据输入到ANFIS中得到的输出；

若E₀≤error，或者达到最大训练次数epoc，结束训练并转到(6.5)，否则进行(6.4)

(6.4)训练参数集

(6.4.1)固定参数集

采用误差反向传播(BP)算法对

进行修正；

(6.4.2)固定参数集

采用最小二乘法对

进行修正；

(6.4.3)固定参数集

采用最小二乘法对

进行修正；

(6.4.4)得到修正后的参数集返回步骤(6.2)；

(6.5)得到确定参数集

用第二个数据集的训练数据对第二个ANFIS进行训练，重复上述训练过程，保持ANFIS的输入数据不变，输出数据将替换Y₁为Y₂，

J＝800，为Y₂中第j+1组数据中的瓶颈设备利用率，K₁将替换为K₂，确定第二个ANFIS的模糊规则数为K₂，训练ANFIS，确定第二个ANFIS的参数集

用第三个数据集的训练数据对第三个ANFIS进行训练，重复上述训练过程，保持ANFIS的输入数据不变，输出数据将替换Y₁为Y₃，

J＝800，

为Y₃中第j+1组数据中的瓶颈设备缓冲区队列长度，K₁将替换为K₃，确定第二个ANFIS的模糊规则数为K₃，训练ANFIS，确定第三个ANFIS的参数集

按照以上步骤用训练数据完成对三个ANFIS的训练，训练完成后。将后500组测试数据输入到ANFIS中，得到预测输出。瓶颈设备编号的训练输出与实际输出对比曲线如图4所示，预测输出与实际输出对比曲线如图5所示；瓶颈设备利用率的训练输出与实际输出对比曲线如图6所示，预测输出与实际输出对比曲线如图7所示；瓶颈设备缓冲区队列长度的训练输出与实际输出对比曲线如图8所示，预测输出与实际输出对比曲线如图9所示；

由图5瓶颈设备编号预测输出与实际输出不重合的点数为13，预测正确率为(500-13)/500＝97.4％；图7瓶颈设备利用率的预测输出与实际输出的平均相对误差为0.09618％；由图9瓶颈设备队列长度预测输出与实际输出不重合的点数为22，预测正确率为(500-22)/500＝95.6％。该预测方法可以准确的预测出瓶颈设备的编号，设备利用率与缓冲区队列长度，对基于瓶颈调度方法有着重要的指导意义。

(7)将实际生产线运行获得的实时数据进行工况分析，将其与原始数据进行比较，数值高于原始数据值区间上限的120％或低于区间下限的80％，则认为该组数据是在新的工况下获得，将该数据与原始训练数据组合重新作为ANFIS的训练数据集，重新对ANFIS进行训练，进行在线调整，使其能在不同的生产状况下保证对生产线瓶颈设备编号及相关参数的准确预测。

Claims (2)

1.一种可重入制造系统瓶颈设备预测方法，其特征在于，包含以下步骤： 

(1)确定影响系统中瓶颈设备的关键因素为工件类型、工件投料方式、设备加工时间、设备故障间隔时间、设备平均维修时间；以及设备利用率、缓冲区队列长度、生产线在制品水平WIP； 

(2)获取以下实际生产线数据： 

(2.1.1)每台设备的设备加工时间、设备故障间隔时间、设备平均维修时间； 

(2.1.2)产品工艺路径：总共O道工序； 

(2.1.3)工件类型、工件投料方式； 

根据上述数据信息构建模型，模型由M台设备及工件投料装置、工件回收装置构成，工件投料装置按照某种投料方式进行投料，投放W种类型工件进入生产线中，工件按照工艺路径在设备中根据加工时间进行加工，最后由回收装置收集成品； 

对该模型进行仿真，设工件在任意一台设备上完成一道工序的加工为一步，每完成一步进行一次数据采集，得到R₀组数据样本：{QL_s(1),QL_s(2),…,QL_s(t),…,QL_s(R₀)}，{η_s(1),η_s(2),…η_s(t)…,η_s(R₀)}，{W(1),W(2),…,W(t),…,W(R₀)}，{Q(1),Q(2),…,Q(t),…,Q(R₀)}，{T(1),T(2),…,T(t),…,T(R₀)}，t=1,2,…,R₀，s为设备编号，s=1,2,…,M；其中，QL_s(t)、η_s(t)分别为第t步数据采集所得到的设备s的缓冲区队列长度与设备利用率，W(t)为第t步数据采集得到的生产线WIP，Q(t)、T(t)分别为第t步数据采集得到的生产线的投料方式及产品类型； 

除({Q(1),Q(2),…,Q(t),…,Q(R₀)},{T(1),T(2),…,T(t),…,T(R₀)}外，对每一个数据样本进行如下数据修正： 

设采集得到的数据样本为{g(t)},t=1,2,…,R₀，且实际生产线中对应数据的样本均值为ξ，方差为σ，定义数据偏差为D(t)=g(t)-ξ，若|D(t)|≤2σ，t=1,2,…,R₀，则认为g(t)能够反映实际生产线状况，否则剔除第t步所采集的所有数据，得到修正后的数据； 

(3)对上述修正过的数据按下式计算每个设备的瓶颈度，并进行数据再修正： 

式中QL_s(n)，η_s(n)，分别为第n步采集得到的设备s的缓冲区队列长度及设备利用率，BV_s(n)第n步计算得到的设备s瓶颈度，s为设备编号，n=1,2,…,R₀，但其中不包括步骤(2)剔除的数据；T为缓冲数预估上限，θ₁与θ₂为经验参数，θ₁,θ₂∈[0,1]，θ₁+θ₂=1； 

设瓶颈度阈值B∈[0.8,0.95]，若第n步每台设备的瓶颈度BV_s(n)均小于B，则剔除第n步所采集数据，得到再修正的数据样本；若存在BV_s(n)≥B，取瓶颈度最大值所对应的设备作为系统中瓶颈设备，得到第n步瓶颈设备编号及相应的设备利用率、缓冲区队列长度数据； 

(4)确定每台设备的缓冲区队列长度、设备利用率、设备加工时间、设备故障时间间隔、设备平均维修时间，以及生产线在制品水平WIP，工件投料方式，工件类型为三个ANFIS的输入，下一时刻瓶颈设备编号、缓冲区队列长度、设备利用率分别为三个ANFIS的输出； 

(5)采样数据经过数据预处理后，得到H+1组数据，按如下方式构建三个数据集

及确定需要三个ANFIS； 

式中j代表第j组数据，j=1,2,…,H； 

第j组数据中，

为设备1的缓冲区队列长度、

为设备1的设备利用率、

为设备1的设备加工时间、

为设备1的设备故障时间间隔、

为设备1的设备平均维修时间； 

第j组数据中，

为设备2的缓冲区队列长度、

为设备2的设备利用率、

为设备2的设备加工时间、

为设备2的设备故障时间间隔、

为设备2的设备平均维修时间；以此类推； 

第j组数据中生产线在制品水平(WIP)，M为设备数量； 

第j+1组数据中瓶颈设备编号； 

第j+1组数据中瓶颈设备缓冲区队列长度； 

第j+1组数据中瓶颈设备的设备利用率； 

pattern^j：第j组数据中的投料方式； 

type^j：第j组数据中的工件类型； 

以上3个数据集均由H组数据组成，取每个数据集V组为训练数据，V=[0.6H,0.7H]，剩余H-V组为测试数据，按如下步骤对数值型数据进行聚类，确定3个ANFIS结构，并对类别型数据进行编码，ANFIS结构确定包括每个输入变量对应隶属度函数及规则的个数； 

(5.1)对每个数据集中的类别型数据按如下方式进行编码，形成一个行向量： 

(5.1.1)设投料方式有A种，对第j组数据，令

若第j组数据中生产线采用某种投料方式，则对应的行向量中的元素为1，其余均为0； 

(5.1.2)设工件类型有C种，对第j组数据，令

第j组数据中生产线存在几种类型工件，则对应的行向量中的元素为1，其余均为0； 

(5.1.3)将pattern^j,type^j行向量构成向量s^j，即 

则完成对类别数据的编码； 

(5.2)对每个数据集的训练数据中数值型数据采用二叉树及模糊C均值进行聚类，并根据聚类后子空间形成初始模糊推理系统，步骤如下： 

(5.2.1)令二叉树的根节点为LN(1)，LN(1)为训练数据中的数值型数据，LN(1)=[X,Y₁]，X=[x¹…x^j…x^V]^T，X为由x^j构成的矩阵，

x^j为由

构成的行向量，

为数据集中的第j组数据，分别代表每台设备的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间以及WIP，

Y₁为由

构成的列向量，j=1,2,…,V，为第j+1组数据中的瓶颈设备编号；初始化叶节点集合为ψ={LN(M)}，整棵树的节点M=1； 

(5.2.2) 设定误差限R，R∈[0.001,0.05]，对叶节点集合ψ中所有节点采用二叉树法进行判定，集合中每个节点是否需要被继续划分，若节点M的线性拟合误差e(M)大于R，则按照步骤(5.2.3)继续划分，若所有节点的线性拟合误差均小于R，则转到步骤(5.2.4)； 

(5.2.3)采用模糊C均值对e(M)大于误差限R的节点LN(M)进行聚类，详细聚类步骤如下： 

(5.2.3.1)给定聚类中心数为c=2，设定容许误差E_max及加权指数m=2，E_max∈[0.001,0.05]，随机初始化聚类中心； 

(5.2.3.2)若模糊C均值的聚类准则函数值J_m≤E_max，则结束聚类； 

(5.2.3.3)节点LN(M)被划分成两个叶节点LN(M+1),LN(M+2)； 

(5.2.3.4)更新节点集合：将ψ与{LN(M+1),LN(M+2)}合并，更新ψ； 

(5.2.3.5)返回第(5.2.2)步； 

(5.2.4)得到叶节点数据集ψ={LN(1),LN(2),…,LN(K₁)}，LN(1),LN(2),…,LN(K₁)为聚类后形成的数据子空间，K₁为通过上述聚类过程得到的聚类中心的个数； 

(5.2.5)根据聚类中心个数确定模糊系统的模糊规则数为K₁； 

将(5.2.1)中的分别替换成 Y₁为由

构成的列向量，Y₃为由构成的列向量， 

j=1,2,…,V，

为第j+1组数据中的瓶颈设备利用率，为第j+1组数据中的瓶颈设备缓冲区队列长度，重复上述聚类过程，得到每个ANFIS的中模糊规则数K₁，K₂，K₃； 

(6)按照如下步骤确定三个ANFIS中的参数： 

(6.1)初始化参数集

其中称为前件参数， 

i=1,2,…,5×M+1，k=1,2,…K₁，K₁为模糊规则数量，v_i,k,σ_i,k均为(0,1)中的随机数； 

为(A+C)×K₁维矩阵，矩阵中每个元素均为(0,1)中的随机数； 

称为后件参数，k=1,2,…K₁，

均为(0,1)中的随机数； 

(6.2)第一个数据集中训练数据输入到ANFIS中，确定第一个ANFIS中参数： 

(6.2.1)对数值型数据从ANFIS第一层输入，第一层有5×M+1个输入量，输入第j组数据为

j=1,2，…,V，

为数据集中的第j组数据，分别代表每台设备的队列长度、设备利用率、加工时间、故障时间间隔、平均维修时间以及WIP，第一层对输入变量模糊化，输出为对应模糊集的隶属度，这里隶属度函数采用高斯型函数： 

其中i=1,2,…,5×M+1;k=1,2,…K₁;h=(5×M+1)×K₁，K₁为模糊规则数量，

为隶属度函数，代表第一层输出，其中

为前件参数； 

(6.2.2)第二层，采用乘法规则计算每条模糊规则的激活强度w_k，第二层各节点的输出为

(6.2.3)第三层，归一化得到各条模糊规则的激活强度

输出为 

对编码后的类别输入经过一个激活强度转换矩阵T₁，且

输入到第三层F^j=s^jT₁

F^j=s^jT₁                                       (5) 

其中j=1,2,…,V，

为1×(A+C)行向量，

为K₁维行向量，T₁为(A+C)×K₁维矩阵； 

则第三层的总输出为；

j=1,2,...,V; 

(6.2.4)第四层，每个节点的传递函数为线性函数，计算得到每条规则的输出O4,k 

式中

为后件参数，

(6.2.5)计算所有规则输出之和

； 

(6.3)设定最大训练次数epoc及误差限error，epoc∈[100,500]，error∈[0,0.005]，训练数据的总输出误差为： 

式中

为Y₁中第j+1组数据中的瓶颈设备编号，J=V，Y₁为第一个ANFIS 中训练数据中的输出数据，

为第j组训练数据输入到ANFIS中得到的输出； 

若E_o≤error，或者达到最大训练次数epoc，结束训练并转到(6.5)，否则进行(6.4) 

(6.4)训练参数集

(6.4.1)固定参数集

采用误差反向传播算法对

进行修正； 

(6.4.2)固定参数集

采用最小二乘法对

进行修正； 

(6.4.3)固定参数集

采用最小二乘法对

进行修正； 

(6.4.4)得到修正后的参数集

返回步骤(6.2); 

(6.5)得到确定参数集

用第二个数据集的训练数据对第二个ANFIS进行训练，重复上述训练过程，保持ANFIS的输入数据不变，输出数据将替换Y₁为

J=V， 

为Y₂中第j+1组数据中的瓶颈设备利用率，K₁将替换为K₂，确定第二个ANFIS的模糊规则数为K₂，训练ANFIS，确定第二个ANFIS的参数集

用第三个数据集的训练数据对第三个ANFIS进行训练，重复上述训练过程，保持ANFIS的输入数据不变，输出数据将替换Y₁为

为Y₃中第j+1组数据中的瓶颈设备缓冲区队列长度，K₁将替换为K₃，确定第二个ANFIS的模糊规则数为K₃，训练ANFIS，确定第三个ANFIS的参数集 

按照以上步骤用训练数据完成对三个ANFIS的训练，训练完成后；将后H-V组测试数据输入到ANFIS中，得到预测输出。 

2.根据权利要求1所述的一种可重入制造系统瓶颈设备预测方法，其特征在于，还可以在步骤(6)后包括以下步骤： 

(7)将生产线所获得的实时数据，与原始数据进行分析对比，若其值低于原始数据值区间上限的120%或高于区间下限的80%，则将数据直接输入到ANFIS，得到下一时刻预测输出值；若数据值超出原始数据值区间上限的120%或高于区间下限的80%，则认为该组数据是在新的工况下获得，将该数据与原始训练数据组合重新作为ANFIS的训练数据集，返回到步骤(5)，对ANFIS中的参数进行在线调整，得到下一时刻瓶颈设备相关参数的预测值。 

Images