Computational Effective Fault Detection by Means of Signature Functions.

The paper presents a computationally effective method for fault detection. A system's responses are measured under healthy and ill conditions. These signals are used to calculate so-called signature functions that create a signal space. The current system's response is projected into this space. The signal location in this space easily allows to determine the fault. No classifier such as a neural network, hidden Markov models, etc. is required. The advantage of this proposed method is its efficiency, as computing projections amount to calculating dot products. Therefore, this method is suitable for real-time embedded systems due to its simplicity and undemanding processing capabilities which permit the use of low-cost hardware and allow rapid implementation. The approach performs well for systems that can be considered linear and stationary. The communication presents an application, whereby an industrial process of moulding is supervised. The machine is composed of forms (dies) whose alignment must be precisely set and maintained during the work. Typically, the process is stopped periodically to manually control the alignment. The applied algorithm allows on-line monitoring of the device by analysing the acceleration signal from a sensor mounted on a die. This enables to detect failures at an early stage thus prolonging the machine's life.

Nomenclature

As scientific texts use different notations for vectors, matrices, random variables, derivatives with respect to vectors, Table 1 provides the convention applied in this paper.

Table 1

The list of symbols and notations used in this paper.

Symbol	Description
σ	standard deviation
σ2	variance, second central moment
E{ ⋅ }	expectation; average value; expected value; first moment
x, X	scalar values
f(t)	one dimensional signal
〈f(t), g(t)〉	dot (inner or scalar) product of two signals; 〈f(t),g(t)〉=∫T_0T_1f(t)g(t)dt
r × c	size of a matrix: r rows and c columns
A	matrix—bold, capital letter
x	column vector—bold, small letter;x = [x1, x2, …, xm]T
f(t)	column vector of one dimensional signals—bold, small letter;f(t) = [f1(t), f2(t), …, fm(t)]T
∂f(x)∂x	derivative of function f(x) with respect to column vector x; denominator layout; ∂f(x)∂x=[∂f(x)∂x_1,∂f(x)∂x_2,…,∂f(x)∂x_m]^T

Introduction

There are many applications that require on-line fault diagnosis, e.g. to monitor the quality of manufactured products or to improve the reliability and safety of a system.

The literature provides many sophisticated solutions that are, in most cases, computationally demanding, hence problematic to implement in embedded devices where processing power is restricted. Moreover, many scenarios require a strictly limited processing time. Computationally effective solutions can be implemented using simple and low-cost hardware which also cuts development time and overall cost.

The article presents a method for detecting faults by analysing the time responses. The solution was applied in monitoring a moulding process which is commonly used in industry to manufacture objects from pliable material. A moulding device is typically composed of two forms whose inner shape determine the outline of a manufactured object. At the beginning of the process, the dies are shifted towards and pressed against each other with huge force to assure no inbetween gaps. The raw material is injected into the forms which are then detached to eject the solidified object. The proper alignment of the dies is crucial to the process. When the contact planes of the dies are skewed, the impact causes uneven wear-off of the forms and their untimely deterioration. Therefore, the production process must be stopped every now and again to measure the forms’ alignment and recalibrate if need be. A reliable and computationally inexpensive method for automatic detection of the dies’ adjustment is desirable in industry.

The suggested approach makes use of an acceleration signal from the sensor to detect a faulty moulding device. In the application, the response signal can be considered as stationary and can also be analysed solely in the time domain by using the dot (inner, scalar) product whose computational cost is proportional to the number of samples. There is no need to employ transforms such as wavelet or Fourier, as in other approaches. This greatly reduces computation complexity and allows to use simpler and lower-priced hardware. This method does not require classifiers such as a neural network, support vector machine (SVM), etc.

The subject of controlling the dies of a moulding device is not widely described in the literature. However, similar challenges arise in different applications. The excellent paper [1] shows an interesting application whereby a stamping process is monitored for a missing part in the production line. The tonnage signal is analysed with help of the recurring plot (RP) technique [2] which is used to capture a fine deviation of a non-stationary signal. The computation of the presented two dimensional RP matrix is fairly burdensome. The classification is determined on the base of the difference between the calculated RP matrix and a reference one. The overall computation cost is proportional to the squared number of analysed samples (being around several hundreds in the considered application).

The comprehensive work [3] describes an application of detecting several faults in a stamping process. The sampled strain signal is transformed using wavelets to acquire a vector of coefficients. The approach necessitates a classifier to make a decision about the fault. The authors compare the effectiveness of hidden Markov models (HMM), artificial neural network (ANN), support vector machine (SVM), support vector regression (SVR) and a proprietary classifier. The optimal training of a classifier is a problem of its own. The approach is suitable to analyse non-stationary signals with a poor signal-to-noise ratio (SNR), however it is too computationally demanding to be implemented in a microcontroller.

A similar work is also reported in [4] where a number p of autoregressive coefficients need to be computed. Thus the computation cost is p (being at least 4) times the cost of calculating the dot product. Moreover, the classifier based on HMM needs to be trained.

The approach presented in [5] analyses the tonnage signal from a stamping process with the help of wavelet transform and a fairly complex technique called statistical process control (SPC).

Similar solutions are employed to prevent the damage of mechanical systems or warn of a fault, e.g. in engines [6], pumps [7], gear-boxes [8], wind turbines [9], etc.

The typical approach to fault diagnosis is composed of the following three stages:

signal transformation—the output response (or responses) of a device or system is often decomposed by Fourier, wavelet [10-12], (short time Fourier transform) STFT, Wiegner-Ville [13], EEMD (empirical mode decomposition) [7] transforms to acquire a vector of parameters;

reduction of the vector size—the vector contains many irrelevant features that can be neglected to simplify the classification problem. The following methods are commonly used in fault diagnosis: principal component analysis (PCA) [14], kernel entropy analysis (KECA) [11], kernel principal component analysis [11, 15], uncorrelated multilinear PCA, uncorrelated multilinear PCA [16] and others;

classification—the reduced vector is an input to a neural network to perform fault classification. As a neural network can produce not deterministic results [17, 18] due to overlearning and local extrema, often the support vector machine (SVM) is used, as in [10, 19]. The work [7] uses Bayesian network and shows its superiority over neural network or SVM in an application of monitoring a gear pump. Some works also use hidden Markov models, as in [3, 4].

The above-presented methodology is general, flexible and allows to analyse signals produced by non-stationary and non-linear systems. This comprehensiveness, however, incurs computational complexity that in most cases precludes implementation in embedded devices. At the same time, there are systems that can be considered stationary and linear and whose output responses for different faults are distinguishable in the time domain. Therefore, there is no need to apply signal transform which imposes a considerable computational burden. The approach suggested in this paper does not require a complex classifier, such as a neural network, HMM, SVM, etc. The system’s response is projected into the space created by the so-called signature functions, which are computed from the training responses. The location of the signal in that space determines the fault type. The properties of the presented method make the approach suitable for embedded industrial solutions where hardware cost, reliability and simplicity are important factors.

Materials and Methods

Derivation

In this section, we will develop a method for determining the system’s state basing on its time response. Thus, the responses pertaining to different conditions should be distinguishable in the time domain. This implies that the analysed system should behave as stationary and linear.

First, the responses corresponding to different conditions of the system need to be acquired. This constitutes training data for the method. Let a(t), where i ranges from 1 to A, denote the system’s responses measured in its in healthy condition. Let b(t), where i ranges from 1 to B, represent the system’s responses under the first type of malfunctioning. Let c(t), where i ranges from 1 to C, represent the system’s response under the second type of malfunctioning. By analogy, we can define other functions representing the system’s responses under other types of malfunctioning. The total number of states (healthy and malfunctioning) is denoted by ξ. The total number of time responses is

For the sake of compactness, we arrange the training signals in the following vector of functions whose dimension is Θ × 1.

On the base of the training signals, we want to calculate so-called signature functions. Let denote the signature of the healthy state. We demand, that be similar to signals a(t) (i = 1…A) and at the same time unrelated (orthogonal) to signals b(t), c(t), …. Consequently, signature function should be similar to signals b(t) and orthogonal to signals a(t), c(t), …This simplifies greatly the classification problem. The current response h(t) of the system is registered and compared to signature functions , , …Then, for instance, in the healthy state, h(t) bears a resemblance to and is orthogonal to other signature functions.

We assume that the signatures can be expressed as linear combinations of the training signals, hence can be written as where unknown vectors x, x, x, …are of dimension Θ × 1.

The challenge is to compute signature functions, which is equivalent to calculating vectors x, x, …so that the signatures are robust against noise and provide minimum classification error for unseen data. We will come back to this problem later.

Another issue is to measure the resemblance between two signals. For that purpose, various techniques are employed, e.g. the sum of squared differences, the sum of absolute differences, dynamic time warping.

For our purpose, we employed the dot (inner, scalar) product, which is a linear operation. Let 〈g(t), h(t)〉 denote the dot product of two functions, defined here by where (T0, T1) denotes the time range of the relevant system’s response.

Let h(t) denote the system’s response. If the system is healthy then

When the system exhibits a fault of the first type, then

By analogy, in case of the second type of failure

In other words, a signal space is spanned on the signature functions. The location of the system’s response h(t) determines its state. Fig 1 shows an example of projecting an analysed signal h(t) on the plane spanned by two signature functions and . The projection of h(t) reads . This vector has the following coordinates in the considered signal space

Fig 1

An example of projecting a signal h(t) on the two-dimensional space spanned by two signals and .

The projection of h(t) reads . The coordinates of the head of this vector are denoted by P.

Signal h(t) can also be considered as a point whose coordinates in the signal plane read

Given the system’s response h(t) and calculated coordinates P(h(t)) in the signal space, the problem of state classification is then easy. Checking the location of P(h(t)) with reference to the decision boundaries is straightforward. An example criterion can be based on the maximum value of c, c, c, ….

As indicated earlier, the challenge is to determine the signature functions , , …that is to calculate vectors x, x, …of unknown coefficients—see Eqs (3–5). The resultant signatures should be resilient against noise and provide a minimum misclassification rate for unseen signals.

The task is similar to regression analysis, whereby a curve equation is calculated on the base of the training data. When the underlying physical model is unknown, the problem is difficult. The curve should provide a minimum error for unseen data. Overfitting the curve to the training data results in a poor fit to the unseen data—this is illustrated in Fig 2. This problem is discussed in detail in [20] (especially Examples 4.3 and 4.4) and in [21]. In the case of neural networks, this phenomenon is described as overlearning.

Fig 2

Illustration of overfitting.

Blue bullets—training data; red circles “unseen” (verification) data; continuous blue line—polynomial fitted accurately to 7 training points; red dashed lined—polynomial of optimal order fitted to the training data.

Calculating the signature functions in the way presented below leads to overfitting to the training data. However, this is an instructive step, as it suggests a solution. Demanding for and using the linearity of the 〈·, ·〉 operator, leads to the following matrix equation where is a Θ × 1 vector and D is the following symmetric, Θ × Θ matrix Then and

The signature function calculated in this way is perfectly fitted to the training signals. Fig 3a, 3b and 3c present example training signals a(t), b(t), c(t) and Fig 4a, 4b and 4c corresponding signatures which are noisy. The interpretation of the shape of these signals is difficult. Moreover, the signatures should assume non-zero values only where the training signals representing different states differ. For example, the training data representing healthy and ill states in Fig 3 start to differ only at around t ≈ 0.2612 s, however the overfitted signatures are of significant values already at t ≈ 0.2605 s. This causes a poor fit to the unseen data and is illustrated later in the Result section.

Fig 3

Plots of functions from the training set.

a) a(t)—signals measured in healthy state (i = 1…6); b) b(t)—first type fault (i = 1…6); c) c(t)—second type fault (i = 1…6).

Fig 4

Plots of the overfitted signature functions (calculated on the base of the training set f(t)).

a) under healthy conditions , b) under the first type fault and c) under the second type fault .

To mitigate this problem, we need to constrain the signatures to take non-zero values, solely where the signals representing different states are different.

We demand the subsequent orthogonality conditions where E{ ⋅ } denotes expectation (average value). This constraint prevents the signature functions fit the training data perfectly. Instead, the average residual error equals 0. This mitigates the influence of noise on fitting. Defining the vectors whose dimensions are Θ × 1, the conditions (23–25) can be expressed by or where E = [e, e, e, …] is of size Θ × ξ and v = [1, 0, 0, …] is of size ξ × 1. By the same vein, conditions (26–28) can be written as where v = [0, 1, 0, …]. Vectors x, x, etc. cannot simply be calculated from Eqs (33 and 34) which pose an under-determined set of equations (E is a non-square matrix). We also impose the following minimum energy constraints which have the beneficial effect of smoothing the signature functions. This phenomenon is described as “the blessing of smoothness” and improves the prediction capabilities for unseen data, as explained in [21]. This constraint plays an important role—the signatures will assume non-zero values, only where the training data corresponding to different states bear differences.

Hence, the problem is to minimize subject to Eqs (23–25). Minimizing subject to Eqs (26–28) is separable from the previous problem. We will use the method of Lagrange multipliers [22] and form the following Lagrangian where λ = [λ, λ, …] is the vector of Lagrange multipliers. The derivative should vanish implying that the minimum of was met and the conditions (23–25) satisfied.

The derivatives of conditions (35 and 36), as can be shown by direct computation, read

Differentiating conditions (23–25) yields

Hence

Thus, to calculate x, the set of two matrix Eqs (33 and 44) needs to be solved. For the sake of readability, we rewrite these equations where E is a non-square matrix of size Θ × ξ; D is a square matrix of size Θ × Θ. x cannot be determined directly from Eq (46) by inverting the non-square matrix E. Instead, solving for x in Eq (47) yields

Substituting this to Eq (46) and calculating λ leads to

Plugging λ back to Eq (47) finally results in

The function can now be calculated from Eq (3) or from

By forming an analogous Lagrangian we can calculate x and consequently and so forth.

Variance

Estimating whether the system is healthy or malfunctioning on the base of its response h(t) by the following scalar products is encumbered with an error. The sources of this error are measurement noise, stochasticity (non-repeatability) of the monitored process and differences in the training responses a(t), b(t), etc. Thus c, c and so on can be treated as stochastic variables, for which we will calculate the variance. On the base of conditions (23–25), (26–28) we can calculate the unbiased estimator of the said variance where, for the sake of simplicity, the result of is understood as the following vector

Parameters , , … will serve as measures for the goodness of fit to unseen data.

Results

Application

The proposed method was applied to monitor a bottle moulding process. A simplified diagram of the device is presented in Fig 5.

Fig 5

A simplified diagram of a moulding process being monitored.

“A” stands for an accelerometer, “PC” for personal computer where the signals from the measurement card were stored and analysed.

The system is composed of a static and moveable form (die) where the latter is shifted by an actuator, till both dies meet. The forms are pressed against each other with strong force. Then, a preformed plastic bottle is blown and expanded into the inner shape of the forms.

The fastening screws provide proper alignment of the forms. Over time and due to high pressure, the fastening can become loose which results in an early degradation of the forms—an exaggerated situation is shown in Fig 6.

Fig 6

A sketch showing an exaggerated situation of misaligned forms.

The second type of fault that occurred, albeit less important, was the loose fastening of the actuator to the body which caused additional vibration due to resonance.

Measurements

An analog piezoelectric accelerometer was placed on the moveable form and the acceleration signal along the axis of movement (horizontally) was registered. The parameters of the accelerometer were as follows:

3 dB bandwidth: 1 Hz to 10 kHz,

acceleration range: ±100 g,

sensitivity: 50 mV/g,

where g denotes the gravity of Earth, i.e. g ≈ 9.81 m/s2. The analog signal from the accelerometer was amplified, filtered and sampled at 65536 Hz by an acquisition card. Due to external noise, namely other machines in the factory producing vibrations, the accuracy of the measurements was degraded. The resultant noise was measured during the machine idle time and standard deviation equalled σ = 1.04 m/s2. The noise histogram is presented in Fig 7.

Fig 7

Histogram of the measurement noise (acceleration in a horizontal direction).

Signature functions

We measured 60 system’s responses—20 for each state (healthy, misaligned forms, loose fastening of the actuator). We divided this set into two parts: the training set containing 18 functions and the verification set containing 42 functions. Fig 3a, 3b and 3c show training signals measured when the system was calibrated, the forms misaligned and the actuator’s fastening screws became loose, respectively. The signature functions , and of the healthy system and under two types of malfunctioning respectively, were calculated on the base of the training set f(t) and are presented in Fig 8. The verification set f(t) is used to assess the results.

Fig 8

Plots of signature functions.

a) under healthy conditions—; b) under the first type of fault—; c) under the second type of fault—.

Table 2 shows the results of the standard deviations of c, c and c (see Eqs (54 and 55)) calculated for the verification set of functions f(t).

Table 2

Standard deviations of c, c and c (for the verification set f(t)) are calculated as σ, σ and σ—see Eqs (54 and 55).

σa	σb	σc
0.146	0.134	0.056

For each function from the verification set f(t), the following projections (see Eq (17)) into the signal space were calculated

Points P(f(t)) for the functions from the verification set f(t) are visualized in Fig 9. The device’s state can simply be determined by analysing the location of the P(h(t)) point whereby h(t) is the current response. Perpendicular planes can be drawn to denote the boundaries for different conditions, which in the considered case, are easily distinguishable.

Fig 9

Projection of points P(f(t)) for the functions from the verification set f(t) on 3 planes.

Big green spheres—projections of points P(a(t)); red cubes—projections of P(b(t)); small blue spheres—projections of P(c(t)).

The algorithm using 32-bit fixed-point arithmetic was implemented on a simple ARM7 microcontroller (AT91SAM7x [23]) clocked at 50 MHz. The dot product of the signature functions , and with the device response h(t) required ca 0.83 ms of the microcontroller’s time. There were only around 2750 fixed point multiplications and additions.

Overfitting

The problem of overfitting was mentioned earlier during the justification of the applied assumptions. As indicated before, the overfitted signature functions are calculated (on the base of the training set f(t)) from the following equations

Solving for the unknown coefficients x, x, x we obtain and finally the overfitted signatures equal where is the signature under the system’s healthy condition; is the signature of the first type fault; and —of the second type fault. The plots of the overfitted signatures were presented in Fig 4. The overfitted functions are less smooth. The cause-observation relationship is obfuscated. For instance, at t ≈ 0.2605 s the signature assumes significant non-zero values, whereas the healthy a(t) and faulty b(t) responses differ significantly only after t > 0.2612 s—compare Figs 8 and 4.

Table 3 shows the results of the standard deviations of c, c and c (see Eqs (54 and 55)) calculated for the verification set of functions. The results are worse than in the previous case (whereby , , signatures were used)—cf. Table 2.

Table 3

Standard deviations of c, c and c (for the functions from the verification set f(t) and overfitted signatures) are calculated as σ, σ and σ—see Eqs (54 and 55).

σoa	σob	σoc
0.185	0.227	0.023

By analogy, we define the projections (see Eq (17)) for the overfitted signatures

Points P(f(t)) for the functions from the verification set f(t) are shown in Fig 10. The P(a(t)) points are well fitted for the healthy state, however in the malfunctioning states the fit is poor and characterized by large variances—compare Tables 2 and 3. In this case, the decision boundaries cannot be represented by simple perpendicular planes. The points are more dispersed. Overfitting resulted in poor prediction for the verification responses.

Fig 10

Projection of points P(f(t)) for functions from the verification set f(t) on 3 planes.

Big green spheres—projections of points P(a(t)); red cubes—projections of P(b(t)); small blue spheres—projections of P(c(t)).

Discussion

A diagnosing method for an industrial moulding process of was presented, whereby the device was under three conditions: healthy, first and second type of fault. The solution proved to be effective in the considered application, where the device could be modelled as a stationary and linear system.

The solution can be adapted for different applications and expanded to diagnose multiple faults by creating corresponding signature functions, providing the stationary and linearity requirements are met.

The strong point of this approach is its simplicity for determining the system condition. The system response is projected into a signal space to determine its state. Thus, there is no need for a classifier such a neural network, hidden Markov models (HMM), support vector machine (SVM) for which the training phase can be a complex process.

Moreover, the underlying dot product operation is comparatively undemanding, as opposed to wavelet or Fourier transforms employed by other approaches. Therefore, a simple and low-cost processing unit can be applied to monitor an industrial process.

The presented solution can perform poorly when the responses are periodic with a greatly varying period or phase and also when the analysed system cannot be treated as stationary. In this case, the presented method can use results from the short time Fourier (STFT), wavelet transform, etc. to create multidimensional signature functions and benefit from a simplified classification. For example, employing STFT, 3 dimensional signature functions (amplitude, frequency and time) can be calculated. Such signatures can be robust against the above mentioned obstacles.

Another limitation occurs when one set of training functions is a scaled version of another function set. If this is the case, then an inconsistent system of equations is obtained, e.g. the following equations cannot be simultaneously satisfied: where k is some constant different from 0. To detect these situations, the linear dependence of the training functions needs to be examined. The dimension of the signal space (i.e. number of signature functions) must be less than the total number of the system’s states. The problem needs to be slightly reformulated, however can be handled by the presented method.

A non-linear relationship between the faults and the output responses can compromise the solution, as the coexistence of faults does not correspond to the superposition of the responses under pertaining faulty conditions. The challenge of multiple fault diagnosis is difficult in general [24], however arises considerably less often than detecting a single failure.

Measurement data.

MATLAB files. Sampling rate: fs = 65536 Hz; a.mat—6 training functions for healthy state; b.mat—6 training functions for misaligned dies; c.mat—6 training functions for resonance; av.mat—14 training functions for healthy state; bv.mat—14 training functions for misaligned dies; cv.mat—14 training functions for resonance.

(ZIP)

Click here for additional data file.