A Method for Structure Breaking Point Detection in Engine Oil Pressure Data

Abstract

In this paper, a heavy-duty loader operated in an underground mine is discussed. Due to extremely harsh operational conditions, an important maintenance problem is related to engine oil pressure. We have found that when the degradation process appears, the nature of variation of pressure engine oil changes. Following this observation, we have proposed a data analysis procedure for the structure break point detection. It is based on specific data pre-processing and further statistical analysis. The idea of the paper is to transform the data into a nearly monotonic function that describes the variation of machine condition or in the statistical language—change of the regime inside the process. To achieve that goal we proposed an original data processing procedure. The dataset analyzed in the paper covers one month of observation. We have received confirmation that during that period, maintenance service has been done. The purpose of our research was to remove ambiguity related to direct oil pressure analysis and visualize oil pressure variation in the diagnostic context. As a fleet of machines in the considered company covers more than 1000 loaders/trucks/drilling machines, the importance of this approach is serious from a practical point of view. We believe that it could be also an inspiration for other researchers working with industrial data.

1. Introduction

Maintenance procedures for a fleet of the load-haul dump (LHD) mobile machines operated in underground mine are critically important and challenging. To reach the required level of efficiency, modern maintenance is supported by on-board monitoring systems and advanced analytics. The monitoring system is acquiring data in every 1 s (or more, depending on the type of variable) and measured values are stored in local memory. After the shift, when the machine is going back to the so-called machine chamber, the data are automatically transferred via WIFI to the server on the surface (data storage and processing center). Specific analysis is automatically performed and some reports are available via the website. Unfortunately, the variables are very specific and exhibit different characteristics, thus it is required to apply different algorithms to their analysis. Till now, simple comparison with threshold is used, for example, for low frequency temperature data. In this paper, we present an example related to the engine oil pressure measurement. A base for all analysis is the need to detect the structure break point (called also the anomaly or the structure/regime change point), interpret data automatically, and provide simple information to appropriate staff. Engine oil pressure is very different than mentioned temperature data. It is relatively high-frequency variable. It varies from zero to ca. 700 kPa, but we have found that distribution of the process is more informative than just simple particular sample values. It is hard to make any conclusion from the raw signal. If engine oil pressure is too high, there is an automatic safety valve that releases excess oil from the system. A much worse case is when oil pressure is too small—such a situation may lead to the damage of the engine. However, as we mentioned, the simple diagnosis via comparison with the threshold level is impossible. Thus, in this paper, we propose a statistical procedure for the detection of the so-called structure break point in the highly variable random data.

Our strategy is based on the statistical approach that minimizes the impact of local changes in the examined data. We consider the machine as a time-varying system, so we decide not to analyze the data sample by sample but to estimate its characteristics in the work shifts perspective. A daily operation in the considered mine is divided in four shifts, each 6 h. It is due to mining law regulation, specific processes required for mining production and extremely harsh condition underground. Even if some experts may notice changes in the raw signal, it is very hard to provide objective rules to detect anomalies in practice. Due to the mentioned specific operation of the machine in the mine, firstly we segment the data according to the work shifts. We received advice, that operation of LHD during the first and second work shifts is very different than during the third and fourth work shifts (during the night). The 1st work shift in the day is between 6 a.m.–12 p.m.; the 2nd work shift—12 p.m.–6 p.m.; the third work shift—6 p.m.–0 a.m. and the four work shift—0 a.m.–6 a.m. It appears that shifts are different, indeed, however, from our perspective, it is not critical. From the statistical point of view changes in the data related to change of condition are much stronger than changes related to different work shifts. Our approach allow us to process data to the form that contains information about the condition of the system, not about its daily operation parameters variation as speed, load, etc. From a highly time-varying process related to oil pressure variation, we receive a nearly monotonic function describing the change of condition in the system. Such data transformation is readable for maintenance staff and may be used as a basis for maintenance decisions. It has great practical importance as the population of the LHD fleet exceeds 1000 machines.

A novelty of the approach is related to the appropriate processing of raw data. In the data analysis, we use the fact that the empirical distribution of the engine oil pressure time series may change in the subsequent work shifts, which is strongly related to the degradation process of the machine. Thus, we examine the distance measures between distributions of the data corresponding to the work shifts. Here we utilize the divergence functions that measure the distance between two probability distributions expressed in the means of their probability density functions (pdfs). Based on this characteristic we apply the classical test statistics (i.e., Student’s t and Wilcoxon) used in the problem of testing if two samples have equal means. The final result is the detection of the structure break point which corresponds to the moment of the repair registered by the maintenance staff in the considered mine. From the signal processing perspective, the novelty of the approach is related to statistical processing of long term, historical industrial data to achieve unequivocal information about the change of condition. We propose analyzing data in segments to cancel out local disturbances and extract patterns seen as the probability density function. Then, the difference between segments is estimated by the distance between distributions and such a new feature is used to establish the structural break point in the process. We propose the original data processing procedure that covers pre-processing step (data validation, missing values handling, re-sampling, reshaping data to associate them with specific nature of the process), data analysis (probability density function of data instead of the sample by sample analysis, measures of distance in multidimensional space and finally structural break point detection—i.e., the identification of moment when the process change its nature) and the presentation of the results (transforming obtained information back to real samples domain to visualize change moment). The general concept of the procedure is demonstrated in Figure 1. Namely, the raw data is subjected to a three-stage procedure (successive blocks of the diagram), as a result of which we detect the structural break point. The steps highlighted in each of three stages are described in detail in the following sections, namely, the data pre-processing is presented in Section 2.4, the analysis is described in Section 3, and the visualization is presented in Section 4.

1.1. Brief State of the Art

LHD machines are commonly used in the mining industry and are critically important for the production process. There are many papers related to various aspects of LHDs. Most of them are focused on reliability analysis [1,2,3,4,5,6], prediction and assessment of machines breakdowns [7], machine performance measures [8,9,10,11,12], utilisation analysis, optimisation, production analysis [11,13,14,15], risk evaluation [16], residual life estimation [17], etc. Just a few works have mentioned condition monitoring [18,19,20,21,22], technical and operational aspects related to LHD machines [23,24,25,26] and real measurement of output torque, identification of operation regimes, etc. [27,28,29,30,31,32]. LHD machines are equipped with on-board monitoring systems. They are used for control as well as for maintenance purposes. As one receives real data from the machine, it opens many opportunities for their analysis.

As it was mentioned, from the mathematical/statistical perspective, the considered problem is related to the segmentation of the data which, is strongly connected with the structure break point detection. In recent years, substantial works on segmentation methods for different applications appeared in the literature. A few interesting applications include condition monitoring [33,34] (where structural break detection method based on the adaptive regression splines technique has been proposed to recognize a change of operational regime in copper ore crusher vibration and local maxima method has been proposed in the time–frequency domain for spike detection in bearings vibration, respectively), biomedical signals (e.g., electrocardiogram) [35,36,37,38,39] (where hidden Markov models, moving average and Savitzky-Golay filter, cepstral analysis, wavelet transform, envelope-based segmentation, etc., have been used), speech analysis [40,41,42] (where nonlinear speech analysis based on the microcanonical multiscale formalism, adaptation of Appel and Brandt algorithm, innovation (Shur) adaptive filter have been discussed), econometrics [43,44] (where the regime switching model is applied), and seismic signals [45,46,47,48,49] (where, among others, cumulative sum of Gaussian probability density functions and Markov regime-switching models as well as the empirical second moment of given raw signal have been proposed for seismic signal segmentation in order to extract seismic events).

Many segmentation methods are based on simple statistics in time domain, the cumulative squared data or empirical second moment [50,51,52,53]. However, one can also find methods based on the representation of the data in different domains, such as time-frequency [34,54]. See also the effective segmentation methods used in the physical sciences, like the methods based on the so-called recurrence statistic [55,56,57,58,59]. The method proposed in this paper is based on the relatively simple statistics, however, they are applied to the characteristics of the data describing the distribution in the work shift.

The problem discussed in the paper could be seen as a process diagnosis (fault detection in the process). Various approaches to fault detection and diagnosis in process data have been developed over decades. They can be divided into three main categories: data-driven approaches, deep-knowledge-based approaches, and analytical-model-based approaches [60]. An analytical-model- and deep-knowledge-based methods rely heavily on fundamental knowledge of the process [60]. Unfortunately, in our data, it is very difficult to provide a model of the process as a pressure variation is depending on the behavior of the operator, the environment, performed task, etc. It leads to—using data science language—highly non-Gaussian, nonlinear, non-stationary data. A data-driven methods specifically refer to methods that rely purely on operational data without using process knowledge. Data-driven approaches learn from history and place no requirements on models or expert knowledge [60].

Data-driven process monitoring or statistical process monitoring applies multivariate statistics and machine learning methods to fault detection and diagnosis for industrial process [61]. Various variants of PCA as recursive PCA (RPCA), dynamic PCA (DPCA), and kernel PCA (KPCA), have been used for process diagnosis. The main idea is that structure of PCA will change when an abnormal situation will appear in the data. For data-driven approaches, methods based on the hidden Markov model, Ddynamic neural network, kernel independent component analysis, Gaussian mixture model, hidden semi-Markov model have been applied. A deep review has been provided in [60].

In the context of real-time analysis one can consider two approaches: decision “sample by sample” or decision “segment by segment”. In the first case, we can compare the amplitude of the incoming sample with the threshold or we can test if the sample belongs to the same distribution as previous samples. One may search for pre-defined events (zero value, flat line, min/max, outliers detection, etc.) [62]. It is very difficult to apply this approach in our case. In the second-mentioned technique, one may play with more reliable statistical approaches, instead of samples the statistics or the distributions may be compared. We are doing this in our research. We test if the mean value estimated from at least two new samples is different than the historical one.

It should be also mentioned that there is a class of data-mining based solutions, which could be adopted here, for example-time-series clustering, see review [63], anomaly/novelty detection in time series, [64], see review [65,66], process mining [67]. We consider checking this in near future, however, we believe that the proposed techniques are simple and quick and do not require advanced training algorithms as most of the data mining techniques.

As an interesting approach one may find in [68] where historical SCADA data for the normal condition from wind turbines were used to train this multi-layer network model layer-wise to extract the relationships between SCADA variable trained for the healthy case, the model will not be able to reconstruct real data if input will be abnormal. Thus, residual data will be higher, and easy to notice an anomaly. As mentioned it used an advanced neural network-based model rather than simple mean testing as we proposed.

1.2. Structure of the Paper

The paper is organized as follows: in Section 2 we describe the machine, experiment, and the data as well as we give the step by step pre-processing procedure used for the engine oil pressure. In Section 3 we present the statistical methodology applied for the data processing, it consists of the calculation of the distance measures for empirical probability density functions corresponding to work shifts and finally the application of the Student’s t and Wilcoxon statistics for the structure break point detection. In Section 4 we demonstrate the results for the real engine oil pressure data. The last section contains the summary and concludes the paper.

2. Machine, Obds, Experiment and Data Description

2.1. Machine Description

The machine used in the experiment is a Loader LKP-0903 produced by KGHM ZANAM. Basic parameters of the truck are: length 10,600 mm, width 3150 mm, height 1750 mm total weight 28,700 kg, standard bucket 4.6 ${\mathrm{m}}^3$, tramming capacity 96 kN, power rating 181 kW, driving speed 4.5 km/h (1st gear) up to 19.0 km/h (4th gear). The photo of the machine is shown in Figure 2. An LHD loader is quite a complex machine that consists of several subsystems as the drive unit, transmission system, hydraulic system for lifting the loader bucket, etc. The heart of the machine is the combustion engine. One of the most critical parameters of engine operation is Engine Oil Pressure. For more details see [69]. The loader is used to transport copper ore from the so-called mining face (mining front, where extraction of ore is performed) to the screen (reloading point, from loader to continuous belt conveyor system for transporting bulk material for long distances up to mining shaft and then to the surface). The mining process is complicated and consists of several steps. It implicates the way of the usage of LHD machines. Daily operation is divided into four shifts, but due to several factors (i.e., blasting procedures), the LHD cannot operate continuously. If LHD is not in operation, the monitoring system will produce NaNs (Not a Number values). It may happen that most of the data for a given shift are NaNs. During the weekend LHD is not used, either. All these specific cases lead to dedicated data processing and analysis methodology related to particular data stream from SCADA.

2.2. On Board Diagnostic System Description

To maintain such types of machines a list of critical parameters has been defined to be measured on the machine. The list of parameters is case-depended. In

Table A1. A list of monitored parameters for LHD machines.

Variable Name	Description
	Temperatures
’ENGCOOLT’	temperature of the cooling liquid of the internal combustion engine
’GROILT’	oil temperature of transmission and torque
’HYDOILT’	hydraulic oil temperature
	Pressures
’ENGOILP’	oil pressure of the internal combustion engine
’BREAKP’	breaking pressure
’GROILP’	transmission oil pressure
’HYDOILP’	pressure in the hydraulic system
	Others
’ENGRPM’	engine speed
’FUELUS’	instant fuel consumption
’SELGEAR’	direction and current gear
’SPEED’	average speed every 1s

, see Appendix A, we have present the parameters that are most important for our case study. Parameters are associated with the various components. Their dynamics, variability, etc, are different. Some of the data are sampled with 1 Hz, some others considered as low-frequency variables may be sampled with every 5 s or more. In this paper, we will use oil pressure signal only sampled every 5 s.

The data considered in this paper described one month of LHD operation. We received them from the data server as historical data.

2.3. Experiment Description

The machine considered in the paper is a regular example taken from the mine (LK3 419R machine). The data covers one month of operation. We have selected such a period because we received information that during that period some repair action has been done. Our experiment is a passive one. We just observe the operation of the machine using OBD monitoring data. According to information received from maintenance staff, turbocharger replacement has been done on 14 May 2019. The diagnostic task it to identify the replacement moment based on proposed feature obtained thorough raw pressure data processing. Once again, we would like to highlight that based on observation of raw data we cannot notice this point.

2.4. Data Pre-Processing

Raw engine oil pressure data has been presented in Figure 3 (top). As one can observe the pressure values vary from 0 up to 700 kPa. The engine oil pressure was acquired every 5 s, however, for each day the data acquisition process could start with slightly different time points (exactly at midnight or 1, 2, 3 .. seconds after midnight). The empty spaces indicate the weekends. Thus, only the working days were taken to the analysis. To be able to compare the variability of data during a single shift we re-sampled data to the same time-basis using the linear interpolation. As one may see, it is a cosmetic change—the shape of the signal has not been changed see Figure 3 (bottom).

As mentioned, we know that turbocharger has been replaced on 14 of May but it is still difficult to identify a significant change in data on that date. Similar data were analyzed in [21], where the specific representation of the oil pressure was proposed, namely, the authors reshaped the data into the two-dimensional array where the x-axis describes the number of shift or date for a given day, and on the y-axis we have the so-called local shift time (0–6 h) or time corresponding to a single day (0–24 h). In this paper, we propose to apply the local shift time representation of the one-month measurements, see Figure A1 given in Appendix B. As one can see, it is difficult to clearly indicate the structure break point (14th of May) on the map presented in Figure A1. Thus, we propose to seperately analyze also the first, second, third, and fourth shifts. Moreover, to have the same time for each shift we re-sampled the data according to the most frequent time-basis using linear interpolation. The representation of the oil pressure measurements broken down into four sub-sets (corresponding to four shifts) after re-sampling is demonstrated in Figure A2 given in Appendix B.

However, as was mentioned, the monitoring system produces NaNs. If in the analyzed sample a significant amount of data are NaNs, then the analysis and corresponding interpretation may not be reliable. Thus, we examined how much data (in percent) are NaNs and before the further analysis, we removed the work shifts with more than 40% of NaNs. The oil pressure after removing the days with more than 40% NaNs represented as shift by shift is demonstrated in Figure 4. The oil pressure broken down into four shifts after removing the shifts with more than 40% NaNs is presented in Figure 5. The data demonstrated in Figure 4 and Figure 5 are analyzed using the procedures described in the next section. The consecutive steps of the pre-processing scheme described above are demonstrated in Figure 6.

3. Methodology

In this section, we present the methodology applied for the real data presented in Section 2. At first, we set the notation used further in the paper. Let $m_{ij}$ denote the i-th measurement during the j-th work shift where $i=1,\dots ,N$ and $j=1,\dots ,M$, and consequently let ${\mathbf{m}}_j=(m_{1j},m_{2j},\dots ,m_{Nj})$ denote a vector of measurement corresponding to the j-th work shift. Moreover, let $p_j\left(x\right)$ indicate the theoretical probability density function of ${\mathbf{m}}_j$. To detect the moment when the character of the data changes, we consider the pdfs corresponding to subsequent work shifts, namely $p_1\left(x\right),p_2\left(x\right),\dots ,p_M\left(x\right)$.

In probability theory and statistics, to quantify the similarity of two distributions one can use the so-called divergence (or contrast) functions that measure the distance of one probability distribution to another. In general, the divergence is not a concept as strong as distance, because it does not have to be symmetric in arguments or satisfy the triangle inequality. Among various contrast functions, we distinguish one very important class of divergence coefficients, namely the class of the so-called f-divergences of the following form [70,71,72]

$$I_{f,g}(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))=g\left(\int p_{j^{\ast \ast }}\left(x\right)f\left(\frac{p_{k^{\ast }}\left(x\right)}{p_{k^{\ast \ast }}\left(x\right)}\right)dx\right),$$

(1)

where $p_{k^{\ast }}\left(x\right)$, $p_{k^{\ast \ast }}\left(x\right)$ are the probability density functions corresponding to two variables, $f\left(t\right)$ is a continuous convex real function on ${\mathbb{R}}_{+}$ and $g\left(t\right)$ is an increasing function on $\mathbb{R}$. It is important to mention that the f-divergences are always non-negative and they are equal to zero if and only if the densities $p_{k^{\ast }}\left(x\right)$ and $p_{k^{\ast \ast }}\left(x\right)$ coincide. For more properties of the divergence functions defined in Equation (1) we refer the readers to [72,73]. Depending on the choice of $f\left(t\right)$ and $g\left(t\right)$ one can obtain different forms of the contrast functions. In this paper, we consider three specific measures belonging to the class defined above.

The first measure, called the Hellinger distance, corresponds to the case of $g\left(t\right)=\sqrt{0.5t}$ and $f\left(t\right)={\left(\sqrt{t}-1\right)}^2$ in Equation (1) and it is given by the following formula [72]

$$H(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))=\sqrt{0.5\int {\left(\sqrt{p_{k^{\ast }}\left(x\right)}-\sqrt{p_{k^{\ast \ast }}\left(x\right)}\right)}^{2}dx}.$$

(2)

Let us notice that the Hellinger distance is symmetric with respect to the arguments and obeys the triangle inequality. It satisfies the property that $0\le H(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))\le 1$ with the minimum value corresponding to the case when $p_{k^{\ast }}\left(x\right)=p_{k^{\ast \ast }}\left(x\right)$ for every $x\in \mathbb{R}$, and the maximum value achieved when $p_{k^{\ast }}\left(x\right)$ is equal to zero for every x for which $p_{k^{\ast \ast }}\left(x\right)$ is nonzero and vice versa.

As the second divergence measure, we consider a modification of the Hellinger distance defined above. Namely, for $g\left(t\right)$ being an identity function and $f\left(t\right)$ defined as in Hellinger case, we obtain the so-called Jeffreys distance of the following form [74]

$$J(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))=2H^2(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))=\int {\left(\sqrt{p_{k^{\ast }}\left(x\right)}-\sqrt{p_{k^{\ast \ast }}\left(x\right)}\right)}^{2}dx.$$

(3)

Let us notice that similarly to the Hellinger distance, the Jeffreys measure is symmetric in the arguments. Moreover, it takes values between 0 and 2.

The Chernoff distance, the third example of the divergences conisidered in this paper, corresponds to the case when $g\left(t\right)=-log(-t)$ and $f\left(t\right)=-t^{1-\alpha }$ with $0<\alpha <1$ in Equation (1) and takes the following form [73,74]

$$CH(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))=-log\left(ch(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))\right),$$

(4)

where

$$ch(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))=\int {p_{k^{\ast }}\left(x\right)}^{\alpha }{p_{k^{\ast \ast }}\left(x\right)}^{1-\alpha }dx,$$

(5)

is called the Chernoff coefficient. In general, expect the case of $\alpha =0.5$, the Chernoff distance given in Equation (4) is not symmetric in the arguments nor does it satisfy the triangle inequality. Moreover, let us notice that since $0\le ch(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))\le 1$ we have that $0\le CH(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))\le +\infty$ with the minimum value corresponding to the instance when $p_{k^{\ast }}\left(x\right)=p_{k^{\ast \ast }}\left(x\right)$ for every $x\in \mathbb{R}$. The special case of the Chernoff coefficient given in Equation (5), i.e., the case of $\alpha =0.5$, is called the Bhattacharyya coefficient related directly to the Hellinger distance given in Equation (2), namely

$$H(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))=\sqrt{1-ch(p_{k^{\ast }}\left(x\right),p_{k^{\ast \ast }}\left(x\right))}.$$

Let us note that from the practical point of view to calculate the empirical counterparts of the divergences defined in Equations (2)–(4) there is a need to estimate the probability density functions and use them instead of the theoretical pdfs in the above definitions. In our case we will apply the empirical divergences to the measurements corresponding to the work shifts, namely to the vectors ${\mathbf{m}}_j=(m_{1j},m_{2j},\dots ,m_{Nj})$ where $j=1,\dots ,M$. For this purpose, one can use the kernel density estimator of $p_j\left(x\right)$ defined as follows [75,76]

$${\widehat{p}}_j\left(x\right)=\frac{1}{Nh}\sum _{i=1}^{N}K\left(\frac{x-m_{ij}}{h}\right)$$

(6)

for any $j=1,\dots ,M$ and $x\in \mathbb{R}$, where $K(·)$ is the non-negative kernel smoothing function, and h is the bandwidth. The choice of kernel smoothing function determines the shape of the curve used to estimate the probability density function. In our case, we take the normal kernel that is simply the standard normal pdf of the form

$$K\left(x\right)=\frac{1}{\sqrt{2\pi }}exp\left\{-\frac{x^2}{2}\right\}\mathrm{for}x\in \mathbb{R},$$

(7)

and the bandwidth is chosen using the Silverman’s rule of thumb to be optimal for estimating normal densities [77]. The procedure is implemented in many programming languages and it is available in numerous mathematical packages, e.g., the function “ksdensity” in Matlab.

Since the distribution of the subsequent samples ${\mathbf{m}}_1,{\mathbf{m}}_2,\dots ,{\mathbf{m}}_M$ changes for certain $j^{\ast }\in \{1,2,\dots ,M\}$, we can identify the moment of change by examining the empirical contrast functions describing the similarity of the pdf corresponding to the first work shift and the pdfs corresponding to all the other work shifts in the sample, namely we analyze the following vector

$${\widehat{I}}_{f,g}(p_1\left(x\right),p_1\left(x\right)),{\widehat{I}}_{f,g}(p_1\left(x\right),p_2\left(x\right)),\dots ,{\widehat{I}}_{f,g}(p_1\left(x\right),p_M\left(x\right)),$$

where the f-divergence are specified in Equations (2)–(4). We expect that the values taken by the similarity measures in the sub-samples

$${\widehat{I}}_{f,g}(p_1\left(x\right),p_1\left(x\right)),\dots ,{\widehat{I}}_{f,g}(p_1\left(x\right),p_{j^{\ast }}\left(x\right))$$

and

$${\widehat{I}}_{f,g}(p_1\left(x\right),p_{j^{\ast }+1}\left(x\right)),\dots ,{\widehat{I}}_{f,g}(p_1\left(x\right),p_M\left(x\right))$$

differ significantly. The preliminary analysis of the vector of measurements corresponding to the work shifts indicates the corresponding distributions are different, more precisely, we observe that some characteristics responsible for the location change with respect to time. We refer the reader to the next section for more details. Thus, to determine the structure break point in the data, we decided to apply the test statistics used commonly in the problem of testing if two samples have equal means. However, we apply this methodology not to the raw vector of observations, but to the empirical convergence distance described above. The procedure used to detect the structure break point $j^{\ast }$ in the convergence distances corresponding to the work shifts is described below.

Let us consider M independent random variables $I_{f,g}(p_1\left(x\right),p_1\left(x\right))$, $I_{f,g}(p_1\left(x\right),p_2\left(x\right))$, …, $I_{f,g}(p_1\left(x\right),p_M\left(x\right))$ with the cumulative distribution functions denoted as $F_{I,1}$, $F_{I,2}$, …, $F_{I,M}$ and the following means ${\mu }_{I,1},{\mu }_{I,2},\dots ,{\mu }_{I,M}$. Now, let us examine the location testing problem based on the expected value with the null and alternative hypotheses defined as follows

$$\begin{array}{cc}\hfill & H_0: {\mu }_{I,1}={\mu }_{I,2}=\dots ={\mu }_{I,M}\hfill \\ \hfill & H_1: \exists M^{\ast }\in \{2,\dots ,M-1\}\mathrm{such} \mathrm{that}{\mu }_{I,1}=\dots ={\mu }_{I,M^{\ast }}\ne {\mu }_{I,M^{\ast }+1}=\dots ={\mu }_{I,M}.\hfill \end{array}$$

Let us notice that the alternative hypothesis can be also presented as

$$H_1=\bigcup _{M^{\ast \ast }=2}^{M-1}{H}_{1,M^{\ast \ast }},$$

where $M^{\ast \ast }$ is a fixed value and $H_{1,M^{\ast \ast }}$ is the alternative hypothesis corresponding to the two-sample test comparing means in two populations, i.e.,

$$H_{1,M^{\ast \ast }}:{\mu }_1={\mu }_{I,1}=\dots ={\mu }_{I,M^{\ast \ast }}\mathrm{and}{\mu }_2={\mu }_{I,M^{\ast \ast }+1}=\dots ={\mu }_{I,M}.$$

Now, to verify $H_0$ against $H_1$ we can use the maximum-based statistic of the following form

$$S_M\left(ϵ\right)=\underset{\lfloor ϵM\rfloor \le M^{\ast \ast }\le \lfloor (1-ϵ)M\rfloor }{max}\left|T\left(M^{\ast \ast }\right)\right|,$$

(8)

where $T\left(M^{\ast \ast }\right)$ can be any statistic for testing $H_1$ against $H_{1,M^{\ast }}$ and $ϵ\in (0,0.5)$ is used to guarantee that there are at least $\lfloor nϵ\rfloor$ elements in both sub-samples. Let us notice that the choice of $ϵ$ is crucial since a small value allows us to detect the very early or very late change in the distribution, but at the same time, it involves considering a small size sample which hinders the proper statistical inference. Moreover, it is important to mention that the distribution of the max-type statistic given in Equation (8) under $H_0$ does not depend on the analyzed random sample distribution. In our method, we want to locate the moment of the most significant change in data, i.e., we want to identify the value of $M^{\ast \ast }$ corresponding to the max-type statistic given in Equation (8).

For our purpose, we use two exemplary statistics from two-sample tests for equal means in independent groups, namely the one corresponding to the parametric Student’s t-test and the one corresponding to the non-parametric Wilcoxon test. For our data, the Student’s t statistic $T_S\left(M^{\ast \ast }\right)$ takes the following form

$$T_S\left(M^{\ast \ast }\right)=\frac{(\overline{I_1}-\overline{I_2})}{I_{std}\sqrt{1/M^{\ast \ast }+1/(M-M^{\ast \ast })}},$$

(9)

where

$$I_{std}^2=\frac{(M^{\ast \ast }-1)I_{1,std}^2+(M-M^{\ast \ast }-1)I_{2,std}^2}{M-2},$$

and $\overline{I_1}$, $\overline{I_2}$ and $I_{1,std}^2$, $I_{2,std}^2$ are the empirical means and empirical variances corresponding to the sub-samples

$$I_1=({\widehat{I}}_{f,g}(p_1\left(x\right),p_1\left(x\right)),\dots ,{\widehat{I}}_{f,g}(p_1\left(x\right),p_{M^{\ast \ast }}\left(x\right)))$$

and

$$I_1=({\widehat{I}}_{f,g}(p_1\left(x\right),p_{M^{\ast \ast }+1}\left(x\right)),\dots ,{\widehat{I}}_{f,g}(p_1\left(x\right),p_M\left(x\right))),$$

respectively. In turn, the Wilcoxon statistic denoted by $T_W\left(M^{\ast \ast }\right)$ is given by the following formula

$$T_W\left(M^{\ast \ast }\right)=\frac{W\left(M^{\ast \ast }\right)-\mathrm{E}W\left(M^{\ast \ast }\right)}{\sqrt{\mathrm{Var}\left(W\left(M^{\ast \ast }\right)\right)}},$$

(10)

where $W\left(M^{\ast \ast }\right)={\sum }_{i=1}^{M^{\ast \ast }}{R}_i$, $\mathrm{E}\left[W\left(M^{\ast \ast }\right)\right]=M^{\ast \ast }(M+1)/2$, $\mathrm{Var}\left[W\left(M^{\ast \ast }\right)\right]=M^{\ast \ast }(M-M^{\ast \ast })(M+1)/12$ and $R_i$ denotes the number of elements in a combined vector $I=({\widehat{I}}_{f,g}(p_1\left(x\right),p_1\left(x\right)),\dots ,{\widehat{I}}_{f,g}(p_1\left(x\right),p_M\left(x\right)))$ which are smaller or equal to $i-th$ element of $I_1$. For more information regarding the Student’s t-test and the Wilcoxon test we refer the readers to [78,79,80,81,82,83]. The successive stages used in the data analysis are presented as the diagram in Figure 7. Similar methodology, also based on the above-mentioned tests, is applied by the authors in [84].

4. Real Data Analysis

In this section, we analyze the real data introduced in Section 2 using the methodology described above. At the same time, we consider the data presented shift by shift in Figure 4, and the data divided into the separate classes corresponding to four work shifts presented in Figure 5. To indicate the difference in the distributions, we calculate the empirical probability density functions for the subsequent columns of the data matrices given in Figure 4 and Figure 5. The corresponding plots are presented in Figure 8 and Figure 9 as the two-dimensional graphs where the probability density functions are evaluated using the kernel density estimation as described in the previous section. As one can notice, the shape of the empirical pdfs changes at a certain point. As a consequence, the empirical pdfs take smaller values with greater probability.

Since in Figure 8 and Figure 9 one can see the difference in the shape of the probability density function, to identify the structure break point we use the characteristics based on the distance between the empirical pdfs, namely the Hellinger distance, the Jeffreys distance and the Chernoff distance with $\alpha =0.5$ introduced in Section 3. At first, we plot the maps presenting the matrices of the values taken by the above characteristics where each two empirical pdfs are compared to each other. For the Hellinger distance, the matrices are given in Figure 10 and Figure 11. As one can see they are symmetric to the diagonal. Besides, the calculated Helligner distance divides the work shifts into two groups, between which the values taken by the distance measure differ significantly. The analogous plots corresponding to the Jeffreys and Chernoff distances are given in Appendix B, see Figure A3, Figure A4, Figure A7 and Figure A8.

Now, we apply the Student’s t-based and Wilcoxon-based procedures leading to the identification of the structure break point in the data to the values representing one row of the Hellinger distance matrices presented in Figure 10 and Figure 11. The selected row corresponds to the Hellinger distance between the subsequent pdfs and the base pdf chosen here as the one corresponding to the first shift in a sample. In other words, we identify the structure break point based on the comparison of each empirical pdf with the empirical pdf of the first shift. The values taken by the distance measure and the identified structure break point are presented in Figure 12 and Figure 13 which correspond to the case of the data presented shift by shift and each work shift treated separately, respectively. As one can see, in most of the cases (except the fourth work shift in Figure 13) the procedures based on the Student’s t-test and the Wilcoxon test are consistent and indicate the same structure break point. We mention here that the shift marked on the plots is the last shift before the change in distribution occurs. This means that the nature of the data changes with subsequent work shifts taking place after the one marked on the plot. The similar graphs corresponding to the Jeffreys and Chernoff distances are presented in Appendix B, see Figure A5, Figure A6, Figure A9 and Figure A10. In

Table A2. The identified date corresponding to the last shift before the change in data and the values of test statistics for the data presented shift by shift.

Hellinger Distance
Student’s t	Wilcoxon
1st shift on	1st shift on
14th of May	14th of May
$\mathcal{T}\left(m^{\ast }\right)=-19.76$	$\mathcal{T}\left(m^{\ast }\right)=-7.54$
Jeffreys distance
Student’s t	Wilcoxon
1st shift on	1st shift on
14th of May	14th of May
$\mathcal{T}\left(m^{\ast }\right)=-19.72$	$\mathcal{T}\left(m^{\ast }\right)=-7.54$
Chernoff distance
Student’s t	Wilcoxon
1st shift on	1st shift on
14th of May	14th of May
$\mathcal{T}\left(m^{\ast }\right)=-19.76$	$\mathcal{T}\left(m^{\ast }\right)=-7.53$

,

Table A3. The identified dates corresponding to the last shifts before the change in data and the values of test statistics corresponding to the data representing four work shifts—Hellinger distance.

Hellinger Distance
First Shift
Student’s t	Wilcoxon
14th of May	14th of May
$\mathcal{T}\left(m^{\ast }\right)=-8.39$	$\mathcal{T}\left(m^{\ast }\right)=-3.83$
Second shift
Student’s t	Wilcoxon
13th of May	13th of May
$\mathcal{T}\left(m^{\ast }\right)=-8.53$	$\mathcal{T}\left(m^{\ast }\right)=-3.20$
Third shift
Student’s t	Wilcoxon
13th of May	13th of May
$\mathcal{T}\left(m^{\ast }\right)=-10.73$	$\mathcal{T}\left(m^{\ast }\right)=-3.97$
Fourth shift
Student’s t	Wilcoxon
14th of May	15th of May
$\mathcal{T}\left(m^{\ast }\right)=-7.83$	$\mathcal{T}\left(m^{\ast }\right)=-3.77$

,

Table A4. The identified dates corresponding to the last shifts before the change in data and the values of test statistics corresponding to the data representing four work shifts—Jeffreys distance.

Jeffreys Distance
First Shift
Student’s t	Wilcoxon
14th of May	14th of May
$\mathcal{T}\left(m^{\ast }\right)=-11.21$	$\mathcal{T}\left(m^{\ast }\right)=-3.84$
Second shift
Student’s t	Wilcoxon
13th of May	13th of May
$\mathcal{T}\left(m^{\ast }\right)=-11.93$	$\mathcal{T}\left(m^{\ast }\right)=-3.20$
Third shift
Student’s t	Wilcoxon
13th of May	13th of May
$\mathcal{T}\left(m^{\ast }\right)=-14.76$	$\mathcal{T}\left(m^{\ast }\right)=-3.97$
Fourth shift
Student’s t	Wilcoxon
14th of May	15th of May
$\mathcal{T}\left(m^{\ast }\right)=-9.68$	$\mathcal{T}\left(m^{\ast }\right)=-3.77$

and

Table A5. The identified dates corresponding to the last shifts before the change in data and the values of test statistics corresponding to the data representing four work shifts - Chernoff distance.

Chernoff Distance
First Shift
Student’s t	Wilcoxon
14th of May	14th of May
$\mathcal{T}\left(m^{\ast }\right)=-13.94$	$\mathcal{T}\left(m^{\ast }\right)=-3.84$
Second shift
Student’s t	Wilcoxon
13th of May	13th of May
$\mathcal{T}\left(m^{\ast }\right)=-11.55$	$\mathcal{T}\left(m^{\ast }\right)=-3.20$
Third shift
Student’s t	Wilcoxon
13th of May	13th of May
$\mathcal{T}\left(m^{\ast }\right)=-10.52$	$\mathcal{T}\left(m^{\ast }\right)=-3.97$
Fourth shift
Student’s t	Wilcoxon
14th of May	15th of May
$\mathcal{T}\left(m^{\ast }\right)=-9.68$	$\mathcal{T}\left(m^{\ast }\right)=3.77$

we present the summary of the results together with the values taken by the test statistics corresponding to both procedures (Student’s t-based and Wilcoxon-based) and the data presented shift by shift and each work shift treated separately. As one can see, for the data presented shift by shift, both procedures applied to all three considered measures of distance between pdfs indicate that the structure break point occurs after the 1st work shift on 14 May. On this day, from the 3rd work shift the character of the data changes. It is important to mention that the structure break point occurs between the 1st and 3rd work shift on 14 May because during the 2nd work shift the LHD machine was not operating. The same conclusions can be drawn when we analyze the results corresponding to four work shifts treated separately, besides the case of the Wilcoxon-based procedure applied to the data from the fourth work shifts where the structure break point is indicated only after 15 May. It is important to notice that the choice of divergence measure does not cause any change in the final result for the analyzed data. Although the values taken by the calculated distances differ, the overall behaviour of the functions is similar. As a consequence, we obtain exactly the same output for both procedures (based on Students t-test and on Wilcoxon test) applied to the values taken by Hellinger, Jeffreys and Chernoff distances.

5. Discussion

The developed procedure allows us to transform the one-month data from highly variable to a nearly monotonic feature. It should be said that the behaviour/shape of our feature is different than a bath-curve known from reliability theory. Here, there is a “switch” from good to bad condition. Therefore, using this data we are able just to detect the moment of change, without the possibility to detect the early stage of damage or tracking its development. It is not related to the proposed method, but to specific data.

The a priori knowledge that a repair action has been performed is necessary to validate the results. The detected date of structure break point is the same as the date of maintenance action done in the company. The extracted feature—i.e., the distance between the distributions estimated for each shift is convincing and easy to interpret. From an industrial perspective, it is essential. We believe that it has also appropriate quality from a scientific point of view as we proposed a new method for structure break point detection, i.e., the method for identification of the moment when our feature is changing (that corresponds to the change in the process). This method is based on the Student’s t and Wilcoxon statistics that are commonly applied in testing if two samples have the same means. The base of the procedure is very intuitive, we expect that the distribution of the time series is changing along with the degradation process. Thus, we proposed analyzing the characteristics of the empirical distributions of the data corresponding to the work shifts, i.e., the measures of the distance between empirical probability distribution functions. Using the characteristics, we could obtain the monotonic feature of the data that was further segmented.

The problem considered here is essential from the practical point of view. Even if various SCADA data are increasingly available, practical information for an end-user is still doubtful. It is especially the case for highly variable data acquired from time-varying systems such as mining machines, wind turbines, etc. [85,86,87]. To develop novel diagnostic procedures, historical data are used to assure the “training” process. Often the data are normalized as they depend on operating conditions [85,87]. Unfortunately, we cannot use normalization—engine oil pressure cannot be higher if the machine is more loaded. However, we have noticed that the statistical properties of the process (engine oil pressure variation) are changed when the machine changes the condition. It was the base for our approach.

6. Conclusions

As the conclusion we need to highlight few important issues:

• We have proposed a novel multistep procedure that covers pre-processing, statistical analysis, and visualization.
• The proposed procedure is the novel one. The innovation is related to the combination of the crucial steps of the proposed methodology, namely, the initial segmentation and the representation as a matrix in the work shift perspective as well as the analysis of the characteristics of the data (probability density functions) and the distance measures based on them. Finally, the last step (segmentation) is performed not for the real data, but for the distance measures of the time series’ pdfs. According to our knowledge, this approach is rarely used in real applications.
• The utilization of the distance measures of time series’ pdfs causes that we do not consider here the problem when one parameter of the data changes (like mean or variance). The examined issue is much more general. The analysis of the pdf’s changes causes the algorithm is sensitive to the dynamics of various characteristics of the data, not only the single one.
• The proposed approach is a universal one. It can be used for any cycle that corresponds to the considered phenomena (in our case it is a work shift), to any characteristics of the data (in our case it is the probability density function), and to any distance measure applied to the characteristics (in our case there are distance measures based on the probability density functions).
• The whole procedure is automatic, thus we believe it could be implemented in the monitoring systems used in the company. When the new data corresponding to the next day (four work shifts) come, then using the introduced procedure we can test if they belong to the current regime. In our case, the new sample means the data corresponding to the next day. This approach is often used in monitoring systems. Thus, in some sense, the methodology can be used in a continuous manner.
• The historical data with precise knowledge about replacement has been used for training and validation. Implementation of the proposed method as an automatic data processing procedure should not be a problem for any new machine. Small dataset from a couple of shifts from a new machine (new data set) will be enough to establish the averaged picture of the signature of good condition. If the damage will appear (change of regime), the method will be able to detect it after a few work shifts (min. 2), that is much better than the current situation. Note that a machine with such damage was able to operate for two weeks as there was not a tool to detect the problem.
• It should be highlighted, the proposed methodology has also some limitations. One of this is related to the special requirements of the data. More precisely, there is a need to consider data that could be arranged as work shifts (or any other cycles). This influences the identified structure break point corresponds to the work shift, not the real time point (like hour).