ARIMA算法在工业控制器故障预测的应用

冯 剑,姚罕琦,黄啸虎,胡钦炫

(浙江中控技术股份有限公司，浙江 杭州 310053)

0 引言

工业控制器在现代工业控制系统中，有着举足轻重的地位，是控制系统的大脑。但在实际的工业场景中，高温、湿度、化学腐蚀等恶劣环境，都会让控制器的硬件元件(如电容、电阻、三极管等)加速老化损坏，从而引起控制器功能失效，进而造成整个工业过程控制失效。这可能导致巨大的生产故障和生命财产损失。所以在控制器故障发生的早期，即还未对生产过程造成损害之前，就及时检测出控制器故障，并提前对其进行故障排除和有针对性的维修，对工业安全生产有着十分重要的意义。

工业控制器是复杂的、集成度很高的电子系统，由多种不同功用的电子元件组成。这些电子元件的特点是很难建立寿命模型来描述当前状态与寿命之间的关系，而且元器件之间会互相影响，当1个元器件失效或者半失效会影响其他元器件的状态。工业控制器的上述特性，增加了人们对其故障预测的建模难度。

现代工业系统可以采集海量的过程运行数据，为基于数据的故障预测方法奠定基础。其原理是通过各种数据分析方法挖掘出其中的隐含信息并进行预测操作，从而避免“难以获得动态模型、不适合预测间歇性故障”等瓶颈出现，是1种实用的故障预测方法[1]。其中，时间序列分析是1种重要的现代统计分析方法，被广泛应用于自然、社会、科研等各个领域。它描述的是某一变量自身的统计规律性，根据序列自身的统计相关关系，揭示相应系统的动态结构特性及其发展变化规律[2]。如：有的研究者提取车辆设备振动特征的时序数据，建立了振动故障序列自回归滑移平均模型(auto-regressive moving average model,ARMA)，为进一步提取故障征兆信息及故障发展趋势预测提供了条件[3-5]；有的研究者将隶属监督学习的机器学习算法和时间序列分析算法相结合进行预测，即利用时序算法计算得到学习样本，在此基础上建立神经网络模型或者包含专家经验的模糊数学隶属度函数，从而为预测工业故障开辟了新的道路[6-9]；还有的研究者使用非监督学习算法结合时间序列分析进行预测实践，即结合无监督学习的无需先验知识及人工打标签等优点到时序算法中，取得了一定的效果[10-13]。

在上述研究成果的基础上，本文结合聚类算法，提出了1种用自回归积分滑动平均模型(auto-regressive integrated moving average model,ARIMA)预测工业控制器硬件故障的方法，并通过试验验证了该方法的可行性。

1 预测研究原理及方案

1.1 预测原理

工业控制器的硬件故障，可以从其各个电子原件的运行参数异常来判断。如输入/输出(input/output,I/O)管脚原件老化异常，其电压会偏离常规的范围。本文选取了能表现控制器硬件故障的特征(CPU电压、I/O管脚电压、电源电压)，对每个特征建立相应的ARIMA，从而预测未来值。当未来值不在正常范围内，则认为控制器硬件出现故障。

ARIMA的数据样本要求具有平稳性，否则预测结果会不准确。考虑到各个电子元件在运行期间会相互影响，导致它们的运行参数会有相关性，如CPU电压在某个范围波动，则I/O电压也一定会在某个范围内波动。所以本文将待预测的所有特征的实时值合成1个特征向量作为计算单元，再使用k-medoids聚类算法将每个单位时间内的多个特征向量进行聚类，获得该单位时间的聚类中心，即可获得每个特征在该单位时间的中心值。对连续时间单位进行聚类操作，即可获得每个特征的连续时间序列作为ARIMA的学习样本。

1.2 ARIMA

ARIMA的基本思想是将预测对象随时间推移而形成的数据序列视为1个随机序列，以时间序列的自相关分析为基础，用一定的数学模型近似描述这个序列。这个模型一旦被识别后就可以依据时间序列的过去值及现在值来预测未来值[14]。ARIMA由自回归(auto-regressive,AR)模型、移动平均(moving average,MA)模型、差分(integrated,I)法结合而成。

AR描述当前值与历史值之间的关系，用变量自身的历史时间数据对自身进行预测。首先，需要确定1个阶数p，以表示用几期的历史值来预测当前值。p阶自回归模型的计算式定义为：

(1)

式中:yt为t时刻的值，t取0,1,2,...；μ为常数项；γi为第i阶的自相关系数；εt为t时刻的误差；p为自回归过程阶数。

MA关注自回归模型中的误差项的累加，q阶自回归过程的计算式定义如下。

(2)

式中:q为自回归移动的平均过程阶数；θi为第i阶模型的估计系数。

移动平均法能有效地消除预测中的随机波动。结合AR和MA，可获得ARMA(p,q)。计算式定义如下。

(3)

再结合差分法，即可获得ARIMA(p,q,d)。其中，d是需要对数据进行差分的阶数，即差分项，一般取1阶或者2阶即可。p和q可通过自相关函数(autocorrelation function,ACF)和偏自相关函数(partial autocorrelation function,PACF)预估，再通过贝叶斯信息准则(Bayesian information criterion,BIC)进行最后确认。

BIC是衡量统计模型拟合优良性的1种标准，计算式如下。

BBIC=k× ln(n)-2ln(L)

(4)

式中:k为模型参数个数；n为样本数量；L为似然函数。

1.3 聚类算法选择

ARIMA算法的学习样本是每个时间单元内、能反映当前特征情况的、有代表性的数值。由于工业环境复杂恶劣，采集的数据可能包含一定的噪声，所以需要选择1种既能够去除异常数据又可以提取特征值的算法。监督学习算法虽然可以实现上述目标，但需要大量的对于数据的打标签工作。标签值也来自于人的经验，在实际工业庞大数据量的情况下，隶属于监督学习的算法费时费力且精度不高，不适合本文的场景。以聚类算法为代表的无监督学习算法可根据数据特征自主判断异常值和特征值，无需人工参与，工作量小，精确度也较监督学习算法高。虽然计算复杂度较高，但是本文场景的时间单元内只有30个样本值，所以计算速度并不低。因此，聚类算法是本文场景的第一选择。聚类算法有很多，考虑到本文聚类的目的是通过聚2类获得大多数样本的聚类中心、去除噪声等小数量样本，所以选择基于质心的聚类算法(如k-means和k-medoids算法)较为合适。

k-means算法中选取的中心点为当前类中所有点的重心，而k-medoids算法选取的中心点为当前cluster中存在的1点，准则函数是当前cluster中所有其他点到该中心点的距离之和最小[15]。正因为中心点较之均值更不容易被“噪声”和“极端异常值”影响，所以本文选择k-medoids算法进行数据处理。

k-medoids算法的数据处理流程如下。

①在总体n个样本点中任意选取k个点作为medoids。

②按照与medoids最近的原则，将剩余的(n-k)个点分配到当前最佳的medoids代表的类中。

③对于第i类中除对应medoids点外的所有其他点，按顺序计算当其为新的medoids时准则函数的值，遍历所有可能，选取准则函数最小时对应的点作为新的medoids。

④重复步骤②～步骤③的过程，直到所有的medoids点不再发生变化或已达到设定的最大迭代次数。

⑤产出最终确定的k个类。

1.4 预测研究过程

①生成样本。

选择控制器每秒的CPU电压、I/O管脚电压、电源电压3个硬件故障特征参数，合成1个3维向量。选取每30 s为1个时间单位，其30个向量使用k-medoids算法进行聚2类操作，获得中心向量。获取n个连续时序的中心向量，每个硬件故障特征拥有n个顺序时序样本。

②建立ARIMA。

对每组硬件故障特征样本进行平稳性分析，并作d阶差分，判断其是否平稳，从而确定d值。

观察序列的ACF、PACF，当其均呈衰减正弦波并趋向于零，表现为拖尾性时，则根据其拖尾起始位置确认p和q的取值范围。利用BIC获取最优的p值和q值，从而建立ARIMA。

③优化ARIMA。

对ARIMA的残差进行白噪声检验，如果检验不通过则调整p值和q值，最终建立通过白噪声检验的ARIMA。

④故障预测。

为每个硬件故障特征建立1个ARIMA。利用ARIMA(p,d,q)进行预测，将预测结果与实际样本进行比较，观察预测曲线的拟合度，以及预测故障时间与真实控制器故障时间的一致程度。

2 研究实践与结果分析

2.1 搭建实践环境

本文选择中控ECS 700控制器作为故障预测对象。在控制器的集成电路中嵌入1颗老化的电阻，让其长期运行。由于电阻老化，控制器会在几天之内迅速故障，可大幅度缩短试验所需时间。

编写程序，让控制器以1次/秒的频率通过工业控制网上报前述的3个特征参数给控制通信主机，即CPU电压、I/O管脚电压、电源电压这3个特征参数。主机端存储实时数据到Mysql数据库。研究人员通过工业信息网从数据库获得样本，使用Python3.5进行研究实践，分析预测结果。

试验环境如图1所示。

图1 试验环境

2.2 采集样本

本文从某天13∶10开始采集样本，至第2天8∶15控制器开始出现故障，即控制器上报的CPU电压、I/O管脚电压均超出阈值，至第2天13∶09控制器上报的电源电压也出现故障。

采集样本数为86 400个(24 h)，选择前57 600个(16 h)样本作为学习样本，后28 800个(8 h)样本进行预测结果对比验证。

每个样本都是由每秒的“CPU电压(C)、I/O电压(I)、电源电压(S)”组成向量V=[C,I,S]。

CPU电压的正常波动范围是1.7～1.9 V。I/O管脚电压的正常波动范围是3.1～3.5 V。电源电压的正常波动范围是18～35 V。

对采集的样本进行观察，样本第19小时15分1秒的实时向量为[2.08,3.59,33.34]，即CPU电压为2.08 V、I/O管脚电压为3.59 V，都超过了正常波动范围，可判断为控制器硬件故障；第21小时1分18秒的实时向量为[2.39,3.79,35.8]，即电源电压为35.8 V，也超过了正常波动范围。

2.3 生成时序模型样本

本文选择30 s为1个聚类单位时间，包含30个特征向量，使用k-medoids算法聚2类。以聚集样本多的那类中心点作为本单元的特征向量。不同单位时间内30个特征样本聚2类结果如图2所示。

图2 不同单位时间内30个特征样本聚2类结果

如在某个单位时间内，使用Python3.5实现k-medoids算法对30个特征样本进行聚2类，得到图2(a)。图2(a)中，左下角有27个样本，而右上角只有3个样本。按照取样本多的那类原则，取左下角那一族样本的中心向量为该时间单元的特征向量，表征该时间单元3个参数(即CPU电压、I/O管脚电压、电源电压)的特征。

再如某个单位时间内，同样进行聚2类操作，得到如图2(b)。图2(b)中，左下角有12个样本，而右上角有18个样本。按照取样本多的那类原则，取右上角那族样本的中心向量为该时间单元的特征向量。

如上所述，利用k-medoids算法，对86 400个(24 h)样本进行计算，获得2 880个3维向量。这些向量中，已经除去了有明显噪声的数据，可以反映每30 s时间单位内(即CPU电压、I/O管脚电压、电源电压)的特征。

从2 880个3维向量，分别获得如图3所示的3个特征的顺序时间序列曲线。

图3中，横实线表示特征值的正常范围。

图3 顺序时间序列曲线

2.4 确定ARIMA的参数

本文为3个特征分别建立3个ARIMA。本小节详细论述了CPU电压的ARIMA建模过程，而I/O管脚电压、电源电压的ARIMA建模过程与之相同，在此不作赘述。

2.4.1 参数d

首先检测CPU电压时序样本的稳定性。

取CPU电压时序数据前1 920个数据(前16 h)作为训练样本。如果样本是非稳定时序，则ARIMA预测的结果不准确。利用Statsmodels算法库进行单位根检验(augmented dickey-fuller test,ADF)，监测结果表明单位检测统计量为0.811，明显大于0.05，说明存在单位根。所以CPU电压时序样本是非稳定的，需要进行差分处理。

利用pandas算法库进行1阶差分操作，再次对1阶差分结果进行稳定性检测，得到单位检测统计量为6.93e-28，已经远小于0.05。因此，可认为其是稳定序列，则参数d=1。

2.4.2 阶数p和阶数q

使用CPU电压1阶差分后的样本，绘制CPU电压ACF图及PACF图分别如图4、图5所示。

图4 CPU电压ACF图

图5 CPU电压PACF图

由图4、图5可知，ACF在1阶以后出现了明显的截尾，而PACF在15阶以后出现拖尾情况。则初设p=0、q=1。

利用BIC对模型ARIMA(p,1,q)进行试验，从而获得最佳p值和q值的组合。试验后发现，p=0、q=1时，BIC函数最小，所以认为ARIMA(0,1,1)是最优参数。

因为ARIMA是在假设随机干扰项是1个白噪声的基础上进行的，所以模型的残差需要验证是1个白噪声序列。如果不能验证为白噪声序列，则说明残差中还有有用的信息，需要重新调整p和q的值。

使用1 920个原始样本(16 h)训练模型对训练结果的残差进行ADF，获得单位检测统计量为0.0。这说明残差不具有相关性，ARIMA(0,1,1)训练后的模型可用于预测。

2.5 预测结果验证

使用建立好的CPU电压ARMIA，进行8 h(24 h中前16 h的数据为训练样本)的预测，获得960个预测结果，并与真实样本(24 h中后8 h)进行比较。CPU电压ARMIA预测值与实际值对比如图6所示。

图6 CPU电压ARMIA预测值与实际值对比图

由图6可知，预测值(虚线)与真实值在5 h内较为贴近，拟合效果较好，从第6 h开始偏差较大。这是因为控制器内部硬件故障导致模型突变，先前使用的训练样本和模型参数已经不足以实现较精准预测。

样本显示第二天的8∶15∶31—8∶16∶00这个聚类区间，CPU电压为19.16 V，高于正常范围；预测结果显示第二天的8∶16∶31—8∶17∶00期间CPU电压为19.02 V，高于正常范围。该结果落后真实情况2个聚类时间区间。

由于预测的结果是1个30 s内的中心值，因此对于实际的故障预测效果进行再次验证。从预测结果来看，预测第二天的8∶16∶31—8∶17∶00这个30 s区间出现超过阈值的CPU电压数值，即从第二天的8∶16∶31开始CPU电压不正常、控制器出现故障。真实情况是从第2天8∶15∶00开始，CPU电压出现在阈值之外的情况，从此CPU电压一直异常，即控制器从第二天8∶15∶00开始故障。

所以利用对CPU电压的预测结果进行控制器故障预测，会晚3个聚类时间区间。其中，有1个区间的误差是由于样本采集方式导致的。因为ARIMA的样本是1个30 s区间的中心值，如果CPU电压数值异常(超出正常范围)不是该区间内的主要状态，则该区间的样本值并不会是CPU电压异常值，预测的结果也会是正常值，导致预测结果产生误差。

I/O管脚电压和电源电压预测值与实际值对比分别如图7、图8所示。

图7 I/O管脚电压预测值与实际值对比图

图8 电源电压预测值与实际值对比图

由图7、图8可知，I/O管脚电压和电源电压的预测结果在前期与其样本值的拟合度较高，后期由于硬件故障而超出模型的预测范围，导致预测精度有所下降。

通过细致的观察可以发现，I/O管脚电压的ARIMA预测硬件出现故障的时间节点是在第二天的8∶17∶31，而实际I/O管脚电压出现电压值突变的时间节点为第二天的8∶16∶10。I/O管脚电压的ARIMA预测硬件出现故障的时间比实际出现故障的时间晚了近3个聚类时间区间。

同样的，电源电压的ARIMA预测硬件出现故障的时间节点是在第二天的9∶01∶31，而实际电源电压出现电压值突变的时间节点为第二天的8∶58∶40。电源电压的ARIMA预测硬件出现故障的时间比实际出现故障的时间晚了近4个聚类时间区间。

3 结论

本文探索了1种预测工业控制器故障的方法，即分别对CPU电压、I/O管脚电压、电源电压这3个特征建立ARIMA，并验证了该方法的可行性。模型训练样本使用k-medoids算法。以30 s为1个时间单元进行聚2类操作。研究结果表明，每个特征的ARIMA预测结果与真实结果非常接近，预测特征值出现异常只比真实情况晚3～4个聚类时间区间，即时延在2 min以内。因此，本文的研究结果为预测控制器故障提供了思路。

尽管本文提出的预测方法对于控制器运行故障的中短期预测达到了良好的效果，但对于长期预测的效果不理想，即预测2 d后的故障准确率逐渐下降。这既是ARIMA算法本身造成的，又和每个预测步长的设定有关。然而，模型可根据设定的时间步长实现连续多步预测，并不断更新预测结果，实现长时间预测。同时，本文目前设定的预测步长为30 s，后期工作中在保证预测准确性的情况下将增加预测步长，从而扩大有效预测结果的时间范围。同时，后期将结合灰色预测算法，实现控制器故障的中长期预测。