Anomaly Detection in a Logistic Operating System Using the Mahalanobis–Taguchi Method

: Product delivery via logistic systems is becoming more e ﬃ cient, rapidly and continuously bringing products to the customer. The continuous operation of logistic equipment, however, can lead to mechanical stoppages due to excessive use. To avoid system failures, fatigue in each part of the system should be monitored, enabling the accurate prediction of potential stoppages and thus promoting overall system e ﬃ ciency. To date, various kinds of anomaly-detection methodologies have been proposed. Among them, the Mahalanobis–Taguchi method, which simply describes the extent of a failure using the Mahalanobis distance, has been utilized to detect changes in the mechanical condition of facilities. However, the technique has not yet been applied to anomaly detection in a logistic operating system. In this paper, anomaly detection using the Mahalanobis–Taguchi method targeting the operational characteristics of a large-scale vertical transfer system is proposed and the validity of the method is discussed. The calculation used to produce proper values of the Mahalanobis distance is ﬁrst developed based on simple excitation using a shaker. Mahalanobis distances under conditions of continuous operation of the target vertical transfer system are then obtained; distances for the system in an artiﬁcially damaged condition are compared to values produced under normal conditions, and any signiﬁcant increase is used as an indicator of a problem. The applicability of the approach to a case involving continuous long-term operation is discussed using a simulation in which the target vertical transfer system is in continuous operation over a two-year period.


Introduction
Recent improvements in technologies related to IoT (Internet of Things) have made real-time monitoring of various types of structures possible. In the field of structural health monitoring, sensor-embedding and wireless-transfer technologies have been developed to monitor the structural integrity of buildings and civil engineering facilities. In the logistics industry, as product delivery has become more efficient and less costly, the working time required for inspection and repair necessitated by unexpected problems can be a major issue, especially for an automatic logistic system operating in a twenty-four-hour active environment. In a conventional management system involving mechanical devices, the majority of maintenance work is done by maintenance workers based on their experience and intuition related to anomaly detection, the timing of which is typically driven by a stoppage of the operation due to some mechanical problem that has already occurred. This can lead to non-negligible economic losses. Thus, it is important to detect a deterioration or change in the mechanical state of a logistic operating system early and accurately in order to predict the time of potential failure and to avoid stoppages of the system to the extent possible. To this end, various damage-detection algorithms [1] and their application to mechanical structures have been studied. A number of studies have been conducted on damage-detection methods for rotating equipment such as bearings [2] and gearboxes [3], the failure of which can cause an entire system to be shut down. When detecting damage of machinery, the use of physical data such as the mechanical ones including vibration and the sound, electrical, and thermal ones can be effective as a basic measure [4]. Using acoustic emissions [5] from damaged parts can also provide an essential clue in damage detection. In addition to such basic measures, a variety of methods for detecting damage features by using machine learning such as neural networks [6,7] has been proposed in recent years. Zhang et al. validated the applicability of SIDL (shift-invariant dictionary learning), which is an extension of regular dictionary learning and makes each atom replicated and shifted at an arbitrary position within the signal [8]. Li et al. have suggested a useful method for nonnegative matrix factorization (NMF) to detect damage to a gearbox [9]. This investigation described a feature extraction and feature selection scheme for hybrid fault diagnosis of a gearbox based on S transform, nonnegative matrix factorization, and mutual information. Assaad et al. also proposed a method for detecting abrasion in gearboxes [10], in which a model-based technique for detecting wear in a multistage planetary gearbox used by lifting cranes was proposed, and established a vibration signal model which dealt with cyclostationary and autoregressive models. In each of these cases, estimating changes in mechanical conditions based on a large amount of accumulated physical data is necessary for effective anomaly detection. Agrawal et al. investigated previously proposed anomaly-detection methods based on data mining [11], grouping the techniques into three categories: classification, clustering, and a hybrid method combining the two. According to this scheme, classification corresponds to neural networks, support vector machines, or genetic algorithms used to classify an unknown condition as normal or anomalous by using a classifier learned with a large amount of physical data collected under normal and anomalous conditions. Presenting research examples of classification problems addressed by a machine-learning-based model, Iannace et al. have proposed a noise-prediction method using random forest regression [12] and a noise-detection method [13] for HVAC (Heating, Ventilation, and Air Conditioning) noises using recursive partitioning. However, such a method is difficult to apply in the absence of sufficient anomaly data or if various kinds of anomalous states can occur, since in both cases the opportunity for machine learning will be limited. On the other hand, clustering, which includes the K-nearest neighbor method, is based on the idea of classifying already known data into groups and of identifying any outliers beyond the range of the already categorized groups as anomalies. For example, the clack-detection methods based on clustering algorithms [14,15] have been investigated.
While various methods have been proposed for outlier detection, using the Mahalanobis-Taguchi (MT) method [16,17], by which the operating state of mechanical equipment can be represented with a single scalar value, offers a powerful approach to anomaly detection. The MT method was originally developed in quality engineering and was applied in various engineering fields due to its versatility. In particular, although various studies have been conducted in the field of anomaly detection, Soylemezoglu et al. have produced a system using the MT method to detect bearing damage [18]. Shakya et al. have also proposed a method for monitoring the health of bearings online and for detecting their deterioration [19]. Balsamo et al. verified a method for detecting the state of damage to a structure using the Mahalanobis distance [20], targeting simple systems such as frame structures. Herein, first of all, the feature of the vibration caused in the mechanical system can be described by various kinds of feature values, e.g., the power spectral density obtained by the Wavelet analysis [21] and its statistical values such as the kurtosis and skewness [22]. However, if the changes of the mechanical condition are complicatedly reflected to the power spectral density of the vibration, it is quite difficult to detect the complex change of the frequency characteristics. In the MT method, the Mahalanobis distance, which describes how the anomaly vibration data differ from the normal one, can be calculated based on the multidimensional feature vectors related to each of the vibration data. Therefore, this method has an advantage to describe the multidimensional values such as the power spectral density with a single value. The multidimensional feature vector can be easily composed of the multiband power spectrum with optimal bandwidth, which can be calculated from the power spectral density obtained by Fourier or Wavelet transform. Herein, we can also reduce the degree of freedom of the feature value by exchanging the narrowband frequency characteristics into wideband ones such as the one-third octave band levels.
Although various anomaly-detection methods have been proposed as mentioned above, a few studies have applied them to logistic equipment such as horizontal belt conveyors. The on-line inspection method of conveyor belts by using machine vision has been investigated [23,24]. The longitudinal tear of the conveyor belt is investigated for detection by a multispectral visual method [25]. The remaining lifetimes of steel cord conveyor belts were observed by using magnetic sensors [26]. Although there are some research cases for these belt conveyors, almost all of them are detected by the sensors other than those by vibration. Furthermore, there are few cases related to the anomaly detection of cargo lifters that vertically convey the logistic products. In any case, the techniques for the anomaly detection of logistics systems have not been yet fully studied and have been still in development. One reason for the lack of research cases in this field is that the mechanical characteristics of these operating devices can be rather complicated; such equipment is often composed of a wide variety of machine elements. Notably, however, since the structure of logistic equipment can include relatively flexible structures such as thin beams or columns, its operating state may be sensitively reflected in the vibration characteristics of the structure. A detection method involving vibration would be difficult to apply using supervised machine learning because of the following reasons. Firstly, it is difficult to recover the labels for specific malfunctions. Secondly, the vibration characteristics change dynamically as a function of the operating situation. Thirdly, a variety of unknown anomaly states in the mechanical system have the possibility to happen, and it is difficult to cover all the possible cases in advance by the detection algorithm. Although it is difficult to collect labeled data for the supervised classification, it is possible to detect an acceleration in vibrations with a sensor and to use the measured acceleration directly; differences between the equipment's normal and anomaly states can be observed in its frequency characteristics as well as in the absolute value of the raw acceleration level.
As described, the accuracy of anomaly detection in a logistic operating system using the MT method can be substantially improved by considering not only the absolute vibration level that occurs during the operation of the equipment but also transitions in the frequency characteristics. Then, the MT method makes it possible to clearly and simply detect changes in state, as transitions in the vibration state can be expressed as a single scalar value. In this study, an anomaly-detection method using the MT approach and targeting the operational characteristics of a large-scale vertical transfer unit (hereinafter, a "lifter") is developed. Firstly, the applicability of the method for steady-state vibrations generated from a simple shaker was examined as a basic study. Then, after that, based on this finding, the properties of the lifter in normal and artificially damaged conditions were investigated through the vibration measurement results obtained by setting the vibration sensor on the surface of the lifter. From the calculation results of the Mahalanobis distances for fragment data acquired under three types of anomaly conditions, it was suggested that the operating state of the lifter can be estimated by using the differences of the Mahalanobis distances. Finally, the validity of the proposed method was discussed from the simulation results in which the lifter was continuously operated for a period of approximately two years and the mechanical condition of the lifter's guide wheel changed gradually from its normal state to the anomaly condition.

Method
To determine the state of vibration based on the MT method, the Mahalanobis distance (hereinafter, "MD") is calculated. The unit space generated from the measured data in the normal operational state is used to evaluate the MD measured in the anomaly state. The detailed scheme is described below. Figure 1 shows the components of the automatic storage system targeted for anomaly detection in this study. The general function of the system is to store products scheduled for shipping; storage is automatically managed in the rack section using the shuttle vehicles shown in the figure. To put the products into and to take them out of storage, a conveyor for horizontal movement and a lifter for vertical movement are used. In many cases, the racks and conveyors in the system can be replaced immediately if one or more of their parts fail. Repairing/replacing the lifter, however, is not so easy. The lifter is essential for vertically conveying products and for distributing them to the various storage layers. It is expected that it will deteriorate more quickly than the other system components because of following reasons: Firstly, the lifter operates continuously for 24 h every day. Secondly, it frequently repeats drastic acceleration/deceleration. Thirdly, loading and unloading products on the carriage of the lifter are continuously repeated. If a problem occurs in the lifter, because of following reasons, it can have a relatively large impact on the movement of products to the shelves in each storage layer. Firstly, because the mechanical elements in the lifter are closely related and joined with one another, it often takes time to replace them in the event of a failure. Secondly, the repair of the lifter requires much time compared to that of the horizontal belt conveyor because the lifter is located in quite a narrow space surrounded by tall partition walls and, therefore, it needs to be repaired by a few people. As a consequence, effective anomaly detection for this part of the system is particularly important; if potential failures in the lifter can be predicted in advance, it is possible to minimize the economic loss due to equipment stoppage. In this study, we developed a prediction technique that can detect lifter failure in advance. As shown in Figure 1, as the lifter moves up and down, the carriage on which the product is placed moves up and down with it; in the operation of the lifter, vibration propagates to the entire body, including the steel columns supporting the lifter and other individual parts. With the proposed method, the transmitted vibration is detected by an accelerometer placed at optimal points on the body and changes in the vibration characteristics are effectively captured. Details of the method are described in the next section.

Normalization of Vibration Data for Creation of Unit Space
With the proposed method, normal data are used to generate the unit space; this unit space is then used as the criterion to evaluate the characteristic difference between normal and anomalous vibration. Table 1 summarizes the normal data before and after normalization. The physical properties of the lifter gradually deteriorate with the passage of operating time; the vibration characteristics gradually change as a consequence. The operation properties are expressed in terms of the multidimensional MD, and judgment as to whether the operation is proceeding normally or abnormally is validated by using the criterion.As shown in Figure 2a,b, the vibration waveform reflecting the movement of the carriage is acquired by an accelerometer placed in an appropriate position on the body of the lifter. The details of the positioning of the sensors are indicated in Section 4.1. For simplicity, to obtain basic knowledge regarding the effect on the value of MD of damage to each of the various parts, we limited the targeted lifter movement to a single condition: one-way movement from the bottom to the top of the lifter. (We were aware from actual measurement that the influence of the lifter's movement on its vibration characteristics differed slightly depending on whether the movement is rising or descending.) Samples of the waveform in rising were extracted, and their frequency response was obtained using FFT (Fast Fourier Transformation). As shown in Figure 2c,d, the converted results of the frequency characteristic have several peaks. These are caused by the structural resonance of each part of the lifter. It is suggested that the vibration components in the frequency band around the resonance frequencies are largely influenced by aging of the structure. Therefore, in this study, the remarkable peak frequencies that occur during normal operation are read from the spectrum and the energy levels of the signals in the 1/12 octave band range above and below these peak frequencies are calculated. These energy levels are sequentially calculated for each sample waveform generated by each rising. The calculation results for these energy levels are expressed as x in Figure 2c. For example, in the kth operation, the energy levels from x k1 to x km in m bands are calculated. The unit space of the normal state is generated using the data output from all the energy levels from the first to the kth operation, as shown in Figure 2c. The output data are summarized in section (a-1) of Table 1.

Generation of Correlation Matrix
The correlation coefficient measuring the correlation between u ij 1 (1 ≤ i ≤ k) in the j 1 th band and u ij 2 (1 ≤ i ≤ k) in the j 2 th band is calculated as Since u has already been normalized, the relationships are satisfied and r can be calculated as Finally, the correlation matrix R can be shown as The correlation matrix R will be used in the next section to obtain MD.

Calculation of MD
Normalization is performed on data group y (Table 1(b-1)), the group to be judged as either normal or anomalous, to produce data group v (Table 1(b-2)). Note that the normalization here is not based on the average and standard deviation of each band of data group y, but rather, it is based on the average and standard deviation of each band of data group x. First, the following column vector V i is composed by using the ith data value v ij (1 ≤ j ≤ m): The value of MD for the jth sample can then be calculated as [17] MD Herein, the value of MD indicates the physical distance between the mean value of the unit space created by the normal data group x and the judged data of the group y. Therefore, as the MD is larger, the frequency characteristics of the judged data y is more different from the normal data group x. Therefore, the value of MD is used as the indicator of the abnormality of the mechanical system in this study.

Quantitative Evaluation of MD
If each feature value for the normalized data used to generate the unit space follows a normal distribution, then k·MD v j 2 follows a χ 2 distribution with k degrees of freedom [16,17]. In this study, the number of bands is k = 4, as will be described later. To confirm the statistical properties such as the normality of the data, Figure 3 shows the probability density distribution of the calculation results of k·MD v j 2 in the normal case. Here, it can be seen that k·MD v j 2 follows the χ 2 distribution with four degrees of freedom. The χ 2 distribution with four degrees of freedom has a critical value of 13.3 when the probability level is set at 0.01. Therefore, in this study, the evaluated vibration is regarded as anomalous when MD v j 2 exceeds 13.3/4 = 3.32. In order to make it easier to refer to this level, the threshold value for an anomaly was set to MD 2 , obtained by simply squaring MD.

Basic Study for Calculation Scheme of MD
The purpose of this study is to determine the mechanical characteristics of the lifter operation. However, to use the MT method for the vibration waveform of the various parts of the lifter, which has a complex assembly structure, it was first necessary to confirm the applicability of the MT method in the case of a simple mechanical device. In this study, it was found that the vibration characteristics measured at each part of the lifter have several peaky frequency components that correspond to the modal components of the lifter. Therefore, a basic study was performed using a shaker (Figure 4) that was excited by a signal superimposed with pure tones at ten frequencies. The time waveform of the output vibration was then measured, and the basic validity of the method to obtain MD 2 was verified. The procedure is described in greater detail below. First, a steady-state signal with the frequency characteristics shown in Figure 5a was generated by superimposing pure tones at every 200 Hz in the frequency range from 200 to 2.2 kHz. The shaker was excited by this signal, and the acceleration waveforms of the shaker were continuously acquired over three days. The measured continuous waveforms were cut out every 6.5 s to generate 13,300 data fragments per day for a total of 39,900. Then, the unit space was generated using the data obtained on the first day and the MD 2 for the waveforms over the entire three days was calculated by using the unit space. While, as described here, the vibration frequency components of each frequency band are used to calculate the value of MD 2 , the scheme for determining the bandwidth to efficiently represent the frequency feature of vibration is discussed next through case studies based on two types of selection methods.

MD 2 Calculation Using Frequency Components in Constant Percentage Bandwidth
We sought to calculate MD 2 by first calculating the energy level in each of the 14 constant percentage bandwidths in the 1/3 octave band from 100 to 2 kHz; the value of MD 2 could then be calculated by using the 14 energy levels as the band information for k = 14 in the MT method (Figure 5a). Figure 6a shows the fluctuation of the energy level in each band over the three-day period. While the measurements were recorded under ideal state conditions with little disturbance vibration, slight fluctuations can be observed, particularly in the 200 Hz band. Figure 7 shows the MD 2 calculation results based on the data in Figure 6a. The average MD 2 for the fragment data between 1 and 13,300 on the first day is approximately 1.0. This seems natural, since the data that generated the unit space are evaluated in the unit space. On the other hand, even though the vibration was continuously output in the steady state, the MD 2 values after the second day greatly exceeded 1.0, and a stable MD 2 result was not obtained despite the vibration being essentially in steady state. This caused us to speculate that the calculated values of MD 2 are materially affected by minute time fluctuations in vibration.

MD Calculation Using Frequency Components at Around Peak Frequencies
As noted, when the method described in Section 3.1 was used, the MD 2 values increased even though the equipment was operating in a steady-state condition. To cope with this problem, an alternative method was considered. Specifically, as shown in Figure 5b, the frequency components used to calculate MD 2 were extracted as 1/12 octave bands around the peaks of the frequency characteristics at 200 Hz, 400 Hz, 600 Hz, and 800 Hz. The reason for this is that background noise components other than peak frequency are included in the data when the energy levels in all of the 1/3 octave bands are used; however, these background noise data are not related to equipment malfunction. Therefore, developing a method that uses only the data around the peak frequencies would seem appropriate.
In addition to the abovementioned treatment, the minute fluctuations in vibration observed in Figure 6a may affect the increase in the value of MD 2 . To establish the periodicity of these fluctuations, the results of FFT processing at the sampling frequency of 13,300 Hz in the 200 Hz band with the largest fluctuation among all the bands in Figure 8 was used as a reference. It was observed that the periodic components in the low-frequency range were dominant. For this reason, a high-pass filter with a cutoff frequency of 4.16 Hz in the sampling frequency of 13,300 Hz was applied to the time-series data of the energy level in each band, and the effects of the fluctuations were cut off to the extent possible. Note that the cycle in 4.16 Hz corresponds to periodic components over approximately four hours. The filtered time-series data are shown in Figure 6b. It can be seen that the minute time fluctuations have been eliminated. Figure 9 shows the calculated MD 2 results. The MD 2 values from the 1 to 13,300 time-series data on the first day have an average of approximately 1.0; the average after the second day is approximately 1.0 as well. Thus, a relatively constant MD 2 is observed.  From the results in Sections 3.1 and 3.2, when the MT method is applied to a vibration waveform having peak components, relatively stable MD 2 values of approximately 1.0 during normal operation are obtained by using the energy level in the frequency band around the peak components of the target waveform, making this the preferred alternative. Moreover, a periodic fluctuation component longer than approximately six hours should be cut off.

Validation of MD Calculation Scheme for Lifter Operation
The validity of the method described in Section 3 for detecting the failure of the lifter operation was subsequently tested. To verify the change in vibration characteristics caused by mechanical failure, three kinds of damage were artificially inflicted on parts of the lifter, and their effects on the MD 2 results were examined.

MD Calculation Using Frequency Components Around Peak Frequencies
The vibration characteristics of the lifter were carefully determined. Figure 10 shows in detail the configuration of the target lifter. The lower left of panel Figure 10a shows a pair of assembled lifters, only one of which was operated and targeted in this study. The lifter has a carriage ( Figure 10c) for loading products. The lifter has a belt to move the carriage up and down and a pulley to operate the belt. The pulleys are located at the bottom and top of the lifter, shown in Figure 10b,d, respectively. To monitor the vibration of the lifter, two accelerometers, R1 and R2 (PCB, 603C01, upper limit frequency: 10 kHz within a variance of ±3 dB), were set at the positions shown in Figure 11. As will be described later, the receiver at R1 was set to determine the mechanical conditions of the belt and guide wheel, whereas that at R2 was set to determine that of the bearings. The reason for the location of these sensors are as follows: Firstly, the vibration caused by the contact between the pulley and belt shown in the right figure of Figure 11a can be sensitively captured at the receiver R1. In addition, the vibration caused by the contact between the guide wheel on the carriage (Figure 10c) and the aluminum column easily transmits via the aluminum column and can be sensitively captured at receiver R1 as well. On the other hand, since the rotational vibration caused by the bearing is relatively weak compared to the vibration caused by the above devices that needs capturing on the surface of the bearing shown in the right figure of Figure 11b, each was attached magnetically to a steel plate. Receiver R1 was set on the steel plate surrounding the pulley, while R2 was set on the steel plate near the bearings used when the pulley at the top of the lifter is rotated. Each of the lifter parts-the belt ( Figure 12A), the guide wheel ( Figure 12B), and the bearings ( Figure 12C)-was artificially damaged. For the belt, the green coating material on the side meshing with the pulley was sanded to make the meshing of the pulley and the belt unstable. For the guide wheel, its rubber-coated surface was melted, making the carriage operation unstable. Finally, a vertical force totaling 1 kN was applied to each of the bearings to destabilize the pulley rotation.

Mechanical Effect of Artificial Damage on Vibration Characteristics of Lifter
The acceleration waveforms of the unloaded lifter under three different conditions are shown in Figure 13: operating in the normal state (Figure 13(a1)), operating with a damaged belt (Figure 13(a2)), and operating with a damaged guide wheel (Figure 13(a3)). All of these results were measured by receiver R1. In this study, the lifter carriage continuously moved up and down in the range between 0.3 and 2.8 m above the ground. Each round trip of the carriage took 8.6 s. Thus, the carriage made 418 trips per hour or 10,046 trips per day. As can be seen in Figure 14, the vibration waveforms during the rising and descending phases of operation are slightly different. In this study, only the vibration waveform for the rising phase was used. Looking at the waveforms shown in Figure 13(a1)-(a3), it can be seen that they are quite different depending on the type of damage.
The spectrograms obtained from the time waveforms are shown in Figure 13(b1)-(b3), respectively. When the carriage rises, it first starts to move, then accelerates, reaches maximum speed, then decelerates, and finally stops at the top. As can be seen from these spectrograms, the eminent frequency of the main frequency component of vibration begins to rise at the start of operation, then begins to fall as deceleration starts. Although the spectrogram of the respective conditions has slightly different features, it is difficult to determine the effect of each type of damage on the time-frequency response. Consequently, the results of FFT processing of the waveforms in Figure 13(a1)-(a3) are shown in Figure 15. As can be seen in Figure 15a, the frequency characteristics of the peaks and dips when operating in the normal condition differ from those associated with operating with a damaged belt. Figure 15b also shows some differences between operating in the normal condition and operating with a damaged guide wheel. Figure 13. Acceleration waveforms and their spectrograms for the lifter (a1,b1) operating in normal state, (a2,b2) operating with a damaged belt, and (a3,b3) operating with a damaged guide wheel. As shown in Figure 13(a1) and (b1), the waveform starts to sweep up after around 0.9 s when the lifter starts to accelerate; after that, the waveform starts to sweep down after 1.8 s; then finally, the large fluctuation of the waveform finishes at around 2.8 s when the lifter stops. Figure 16 shows the time waveforms and spectrograms obtained by R2 under the normal condition and under the condition of damaged (impeded) bearings. The differences between the normal and damaged conditions observed in the time waveform and spectrogram are relatively small compared to those shown in Figure 13. Figure 17 shows the frequency characteristics in these conditions. In these results, the characteristic difference in the peak-dip tendencies shown in Figure 15 is not observed; however, in the damaged-bearing condition, the energy level is approximately 5 dB lower than under the normal condition over the entire frequency range.   . Acceleration waveforms and their spectrograms for the lifter (a1,b1) in the normal state and (a2,b2) operating with a damaged bearing: As shown in Figure 16(a1) and (b1), the waveform starts to sweep up after around 0.9 s when the lifter starts to accelerate; after that, the waveform starts to sweep down after 1.8 s; then finally, the large fluctuation of the waveform finishes at around 2.8 s when the lifter stops. The differences between the normal and damaged conditions observed in the time waveform and spectrogram are relatively small compared to those shown in Figure 13. Figure 17. Comparison of the frequency characteristics of the vibration in the normal condition and when operating with a damaged bearing: This result shows that the characteristic difference in the peak-dip tendencies received at R1 shown in Figure 15 is not observed; however, in the damaged-bearing condition, the energy level is approximately 5 dB lower than under the normal condition over the entire frequency range.
The above results suggest that differences in the mechanical state of the equipment from the normal state can be detected from the vibration characteristics under each of the three conditions tested: a damaged belt, a damaged guide wheel, and a damaged bearing.

Calculation of MD 2 Under Normal Condition
To calculate MD 2 under the normal condition, the unloaded lifter was continuously operated in the normal state for nine days. A unit space was generated based on the fragment data (10,046 waveforms) obtained from measurements on the first day; the MD 2 of the fragment data for the remaining eight days (80,368 waveforms) in the generated unit space was then produced using the method described in Section 2.
From the frequency characteristics of the vibration data obtained from R1 during normal operation (Figure 18), the energy levels in the 1/12 octave band around the four prominent peak frequencies at 724 Hz, 1020 Hz, 1378 Hz, and 1744 Hz were obtained and used for the calculation of MD 2 . Note that the frequency components around the lowest peak at 340 Hz were excluded from the calculation since the level fluctuation in each operation was quite large. Figure 19 shows the time series of the energy levels in each of the four selected frequency bands. Since some trends in the time series of the energy levels were apparent, the periodic components of these series were checked using the results of FFT processing at the sampling frequency of 10,046 Hz. The frequency component of 724 Hz that included the largest trend component is shown in Figure 20. Since long-term periodic components were dominant in the quite low-frequency range, a high-pass filter with a cutoff frequency of 4.16 Hz in the sampling frequency of 10,046 Hz was applied to the time-series data in each band and the effects of the fluctuations were cut off. The filtered time-series data of the energy levels and the MD 2 values obtained from these data are shown in Figures 19b and 21, respectively. The trends were removed by the filtering, and the MD 2 values during the nine days of normal operation appear to be relatively constant around an average value of 0.94. Figure 18. Selection method using frequency bands to calculate the MD 2 of the lifter: From these frequency characteristics of the vibration data obtained from R1 during normal operation, the energy levels in the 1/12 octave band around the four prominent peak frequencies at 724 Hz, 1020 Hz, 1378 Hz, and 1744 Hz were obtained and used for the calculation of MD 2 . Figure 19. Time-series data of the vibration in the normal conditions (a) without high-pass filtering and (b) with high-pass filtering: Since some trends in the time series of these energy levels in the respective bands were apparent, the periodic components of these series were investigated using the results of FFT processing, and will be removed by filtering.

Calculation of MD 2 Under Damaged Conditions
The MD 2 values under the various damaged conditions were calculated based on the generated unit space and the filtered time series of the vibration under the damaged condition. Note that the time-series energy levels under each of the damaged conditions were also filtered by the same high-pass filter designed and were adopted for the normal cases in the previous section. The MD 2 results based on over 30,000 data values obtained from the continuous three-day operation of the lifter for the cases of belt and guide wheel damage are shown in Figure 22a,b, respectively. The MD 2 results from 10,000 data values obtained during a one-day operation for the case in which the bearings were loaded by a strong external force are given in Figure 22c. All of these results show an MD 2 that is substantially larger than the threshold value of 3.32. Thus, it appears that the value of MD 2 is severely affected by the various types of artificially induced damage examined in this study, at least to the extent that they were applied here.

Time-Series Simulation of MD 2 Based on Long-Term Operation
We simulated the transition of the MD 2 values based on the measured vibration data collected from the continuous operation of the lifter over a period of two years, assuming that the mechanical condition of the guide wheel will change from the normal state to the anomaly state when damaged to the same extent as described in the previous section over this period. First of all, we have two kinds of vibration data in the normal and damaged conditions that does not include the intermediate transition data in the period between the normal and anomaly terms. Therefore, the transition data were simply simulated assuming that the energy level in each frequency band linearly decays. The required time for the value of MD 2 to exceed the anomaly threshold value was estimated based on the simulated transition. In this case study, we sought to validate that high-pass filtering can remove the short-term periodic trend but can preserve the long-term transition of the MD 2 that appears as a gradual increase in the estimated values over several years. As reported earlier, the short-term periodic components that result in the unnecessary trend observed in the vibration time-series data were removed by high-pass filtering; however, we wanted to validate its usefulness in a practical case involving the long-term continuous observation of MD 2 . Figures 19a and 23 show the time-series energy levels in each frequency band under normal and anomaly (guide wheel damage) conditions, respectively. Here, the transition of the energy levels from the normal condition (Figure 19a) to the anomaly condition involving damage to the guide wheel ( Figure 23) was simulated to linearly change from one state to the other when operating the lifter continuously for 780 days. Note that, for the time-series data in the normal condition, 90 days were added before the 780 days of operation to verify that the MD 2 value in the normal condition was constant at around 1.0. The transition of the time-series energy levels during the simulated 780 days was generated as follows. First, the data for the normal condition for the eight days that were described in Section 4.3.1 were combined repeatedly. The steady-state time-series data for 780 days was then linearly decreased to smoothly connect to the energy levels under the condition of a damaged guide wheel. As an example, the data for the 1378 Hz band are shown in Figure 24. In these data, the energy level decreases monotonically by approximately 5 dB over the 780-day period. Then, as shown in Figure 25, by using the time-series data of the first five days (50,460 data values), the unit space was generated. Next, the 50,460 values of MD 2 for the five days beginning on the second day were calculated and the average value, MD 1 2 , was treated as the first MD 2 data value. This operation was performed in the same manner for all of the data simulated over 870 days, each time advancing the five-day range by one day. The resulting MD i the value of MD 2 changes as the energy level changes from the normal state to the anomaly state.
The MD 2 for the first 200 days is shown in Figure 26b. The threshold value of 3.3 is also shown.   The values of MD 2 do not appear to change from approximately 1.0 through the first 90 days of operation, as the state during this period had not changed. However, the value of MD 2 starts to rise from the 91st day due to the changing state of the device to the anomaly condition; 53 days later (i.e., 143 days from the beginning of the data), MD 2 exceeds the 3.3 threshold. Thus, based on the results of our simulation, where we assumed that, over the course of 780 days, the normal state would change to an anomaly state in which the value of MD 2 is approximately 300 (Figure 22c), we were able to detect at quite an early point-in fact, within the first 7% of the total elapsed time-that an anomaly had occurred.

Conclusions
An approach for detecting anomalies in a logistic operating system using the MT method was proposed. In this approach, the frequency characteristics of vibration acceleration are used as multivariate parameters. The approach was validated through experimental studies involving an automated cargo lifter.
The applicability of the method for steady-state vibrations generated from a simple shaker was first examined. A steady-state signal superimposed with pure tones at 11 frequencies was output from the shaker, and MD 2 values were calculated from the resulting vibration acceleration data. It was confirmed that the MD 2 values for the measured acceleration calculated using the unit space generated from the same data had an average of 1.0 when the state of the installation of the shaker was unchanged. However, it should be noted that only the frequency components in frequency ranges having predominant relative levels such as the modal peak frequencies should be used in the analysis. Furthermore, moderate cycle fluctuations should be removed from the time-series data by using high-pass filtering.
Based on the above findings, the properties of the lifter in normal and artificially damaged conditions were investigated. MD 2 values for fragment data acquired under three types of anomaly conditions were calculated based on the unit space generated from the fragment data obtained from the measurement results under the normal condition over nine days. In the normal condition, the MD 2 values showed an average of 1.0. In contrast, in the anomaly conditions, the MD 2 values substantially exceeded the threshold value of 3.3 (determined using a probability level of 0.01 in a χ 2 distribution with 4 degrees of freedom). This result suggests that the operating state of the lifter can be estimated by using the differences in MD 2 calculated from the frequency characteristics of the vibration.
Finally, the transition of the mechanical condition of the lifter was simulated over 870 days by using the collected measurement data in this paper and was used for the validation of the proposed method. Since we had two kinds of measured vibration data in the normal and damaged conditions that does not include the intermediate transition data in the period between the normal and anomaly terms, the transition data were simply simulated assuming that the energy level in each frequency band decays linearly. Based on this simulation in which the lifter was continuously operated for a period of approximately two years and the mechanical condition of the lifter's guide wheel changed gradually from its normal state to the anomaly condition, the transition of MD 2 was predicted based on the simulated transition of the fragment data. It was observed that the MD 2 value gradually increased in response to changes in the vibration characteristics over 780 days. As a result, we confirmed that the transition to the anomaly condition could be detected quite early in the process.
As a future work, the proposed method will be validated if the transition of the mechanical condition of the target lifter can be accurately detected by the method in a situation in which the lifter is continuously operated for over years. Furthermore, the applicability of the proposed method will be confirmed in more practical and complicated cases in which the lifter is operated with many loading products in a real distribution factory. After these validations are confirmed, the proposed method will make it possible to accurately monitor the transition of the mechanical conditions of the logistic systems. Therefore, the field of application of the detection method will be expanded to monitor the conditions of more general mechanical objects that provide regular mechanical movements, such as the robotics devices for production machine in the factories.