systems
Article
Classifying and Optimizing Spiral Seed Self-Servo Writer
Parameters in Manufacturing Process Using Artificial
Intelligence Techniques
Chaweng Sapapporn 1,2 , Soontaree Seangsri 3 and Jiraphon Srisertpol 1, *
1
2
3
*
Citation: Sapapporn, C.; Seangsri, S.;
Srisertpol, J. Classifying and
Optimizing Spiral Seed Self-Servo
Writer Parameters in Manufacturing
Process Using Artificial Intelligence
Techniques. Systems 2023, 11, 268.
https://doi.org/10.3390/
systems11060268
Academic Editor: Fernando De la
Prieta Pintado
School of Mechanical Engineering, Suranaree University of Technology, Nakhon Ratchasima 30000, Thailand;
chawengs@hotmail.com
Western Digital Storage Technologies (Thailand) Ltd., Bang Pa-in Industrial Estate,
Phra Nakhon Si Ayutthaya 13160, Thailand
School of Mechatronics Engineering, Rajamangala University of Technology Tawan-ok, Chonburi 20110, Thailand;
soontaree.g@gmail.com
Correspondence: jiraphon@sut.ac.th; Tel.: +66-4422-4412
Abstract: This paper presents a technique for evaluating the performance of high-precision machines
and classifying machine conditions in terms of test capability, such as hard disk drive (HDD) signal
writing machines. In general, position errors generated during the signal writing process must be
minimized to ensure high-quality writing. Position errors refer to deviations in the signal writing
process and can be caused by several factors, such as deviations in the performance of the positioner
that result in a position error signal exceeding its control limit. The proportion-al-integral-derivative
(PID) controller must be optimized to minimize position errors. In model-based controller tuning,
an accurate mathematical model is essential. The first step utilizes system identification methods,
including adaptive weight least squares and peak detection, to create a partition resonance frequency
model. This mathematical model is used to determine the open-loop stability, which involves
achieving gain and phase margin at a specific crossover frequency, and the closed-loop dynamic
response, which involves minimizing the discrete Fourier transform (DFT) of the position error
signal. The DFT of the position error signal in each harmonic can be represented as a resonance
peak in the transfer function model. The DFT and other combinations of operating parameters are
analyzed and used as machine learning features. The ANN classifier was also effective in categorizing
the performance of signal writing machines into four classes: 0 (healthy machine), 1 (sensor fault),
2 (loose pushpin), and 3 (tunable machine). The results showed that the classification performance
was sufficient to separate class 1 and 2 for the maintenance process and class 3 for further optimization
achieved using the mathematical model.
Keywords: dynamic responses; system identification; deep learning; classification; HDD signal
writing machine
Received: 26 April 2023
Revised: 19 May 2023
Accepted: 22 May 2023
Published: 25 May 2023
Copyright: © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).
1. Introduction
Thailand is widely recognized as a global hub for hard disk drive (HDD) production.
The manufacturing process of HDDs involves the utilization of advanced technology at
various stages. Among these stages, the HDD testing process holds significant importance
in instilling confidence in the products being sold. The signal writing process, which relies
on the spiral seed self-servo writing machine, necessitates precise position control. In order
to achieve high-quality position control, a proportional-integral-derivative (PID) controller
must effectively handle diverse load characteristics. Presently, a single set of PID controller
gains is used to control multiple positioners. Before deployment in the manufacturing
process, the PID gains are manually tuned by expert designers to meet time and frequency
domain criteria during bench testing in a clean room environment.
Systems 2023, 11, 268. https://doi.org/10.3390/systems11060268
https://www.mdpi.com/journal/systems
Systems 2023, 11, 268
2 of 33
During operation, the signal writing machine operates across a wide range of frequencies within cycle-time limits. Any deviation in the positioner can impact the quality of the
signal writing process, resulting in the position error signal (PES) exceeding control limits.
The reference signal is generated during the HDD manufacturing process. High-precision
machines are employed to regulate the movement position of the read–write head. However, over time, these machines experience wear and deterioration in their parts, leading to
writing signal errors that exceed the specified range. Consequently, this causes wastage
in the production process. In severe cases of machine part damage, the machine needs to
be stopped for repairs. Optimizing the PID gain is crucial to minimize the PES. However,
based on repair history, some machines have significant mechanical defects that cannot
be rectified by optimizing the controller gain. These machines should undergo repairs
instead of being included in the tunable class. It is necessary to separate them from the
tunable class.
Currently, there are two approaches to address these issues. In cases where the signal
writing quality does not meet the control specifications and defective units are identified, a
high-strength permanent magnet is employed to erase the incorrectly written spiral signal.
In cases where the machinery stops, the repair process is carried out in a step-by-step
manner, beginning with the replacement of components according to predefined steps.
We conducted a study and classified the signal writing machines into four classes:
0 (healthy machines), 1 (encoder sensor degradation), 2 (pushpin degradation), and 3 (tunable
machines). Machines in class 0 function well in manufacturing and produce high-quality
signal writing. Machines in classes 1 and 2 should be repaired by replacing the degraded
components. Machines in class 3 do not exhibit any mechanical degradation, but the PES
exceeds the control limit. These machines should be categorized for controller optimization.
In model-based controller tuning, accurately estimating the mathematical model up to
high frequencies is crucial to account for the effect of each resonance. The mathematical
model can be obtained from the frequency response function (FRF), which can be gathered
and used with system identification techniques to acquire the model. Controller tuning can
then be implemented through model-based optimization to meet both time and frequency
domain criteria.
To make informed decisions for each signal writing machine class, machine learning
is employed to classify machine performance. Operating parameters during the signal
writing process that demonstrate a high correlation with the dynamic response of each
machine characteristic are considered as feature inputs for the machine learning model. The
signal writing machines in class 3, which exhibit less deviation from nominal performance,
are selected for further controller optimization.
The approaches used to perform system identification and machine learning classification were reviewed.
E.C. Levy [1] introduced complex curve fitting in the frequency domain. Levy’s
method identifies FRF in polynomial form in an ordinary least squares (OLS) sense and was
later referenced by several other least squares methods. J.R. Marti [2] presented the Bode
technique to estimate the model parameter of transmission lines. The transfer function
model was estimated exclusively using the magnitude response, wherein pole and zero are
the real part within the complex number. Jiraphon Srisertpol and Autsadayut Rodpai [3]
proposed the mathematical model for linear viscoelastic materials. They created several
weighting functions through frequency modification. Srisertpol and Rodpai define stable
systems as those with a negative pole at the left half-plane of an s-domain. The system
with the minimum norm error was selected as the candidate model in each iteration loop.
Abdullah Al Mamun, T.H. Lee and T.S. Low [4] proposed frequency domain identification
as a transfer function model for disk drive actuators. This approach uses the error residue
matrix from the prior cycle as a weighting function for the current calculation cycle. The
results of this technique indicate whether the model is stable or unstable by using the
left half-plane pole. A stable model is required for accurate controller tuning. Further
research is ongoing for frequencies below 50 Hz, due to limitations in data collection.
Systems 2023, 11, 268
3 of 33
Moreover, some resonance peaks above 4 kHz have not yet been identified, and more
iterations may be required. Taku Noda [5] suggests the identification of a multiphase
network equivalent for electromagnetic transient calculations using a partitioned frequency
response that involves iterative re-weighting. The weighting function is the absolute error
from the previous cycle. The partition is set manually at an interval of 1 kHz, which refers
to the closet equidistant point. This method performs well in high resonance mode, with a
smooth gain at low frequencies. In addition, Eduardo Salvador [6] used several adaptive
errors to modify weighting functions. This method works well when the numerator (zero)
equals the denominator (pole). Takashi Yamaguchi, Mitsuo Hirata and Chee Khiang
Pang [7] proposed HDD actuator modeling with high-speed, precision motion control.
This approach utilizes a combination of rigid body plus delay time and a multiplicative
resonance model. A weighted least squares (WLS) technique (fitsys in MATLAB® ) was
used to validate the stable model. The weighting function was set with the Rigid Body
baseline, with gain reduction in high frequency. Chen Peng and Liang Yanbing [8] presented
frequency domain identification in the context of a fast-steering mirror. The model was
fitted by a gradient search type based on the Levenberg–Marquardt algorithm, and the
error model was improved by adding zero to second-order standard forms. Tomaz K.
and Damir [9] presented an analytical method for estimation of a five-parameter model
in second-order with zero plus time delay. The proposed method uses information in
characteristic areas from arbitrary change in the steady state of the process to calculate
parameters of model. The results show good fitting of both time and frequency response.
To classify machine performance, various classification techniques have been reviewed.
Salem M. et al. [10] presented frequency response analysis (FRA) for three-phase
star and delta induction motors to perform pattern recognition and fault analysis. They
analyzed statistical parameters such as the correlation coefficient (CC), absolute sum of
logarithm error (ASLE), standard deviation (SD), and mean square error (MSE) for healthy
motors in both star and delta connections, as well as motors with short circuit faults (SC)
and open circuit faults (OC). The results demonstrated that these statistical parameters
can be used for pattern recognition. I Abdul et al. [11] presented a fault diagnosis method
for wind turbines using a decision tree classification algorithm. One advantage of this
approach is that it does not require transformation or normalization of the input features,
making it easier to trace the root cause of faults in the wind turbine. The resulting decision
tree classifier can be used by an engineer to design simulation scenarios that replicate the
faults and propose mitigating actions in supervisory control and data acquisition (SCADA).
Masayuti S. and Prabhas C. [12] used a support vector machine (SVM) to classify abnormal
assembly processes in hard disk drive manufacturing. They achieved 100% accuracy in
classifying good and abnormal assembly drives using a training set of 500 drives. They
identified that SVM can be used in manufacturing and has advantages over conventional
methods such as the trapping method because it is easier to develop and does not require
knowledge of specific abnormalities to develop a correct program. Aida et al. [13] presented
a fault detection and classification method for power transmission systems using the knearest neighbor (k-NN) classifier. The system was classified into healthy and 10 faulty
types with an accuracy of about 98%. Yordanos D. et al. [14] presented application of
machine learning for fault classification and location in a radial distribution grid. This
was employed to extract useful features from the three-phase current signal. Standard
statistical techniques are then applied to discrete wavelet transform (DWT) coefficients to
extract the useful features. Multilayer perceptron (MLP) and extreme learning machine
(ELM) were used to classify “no fault” and 10 fault types. The results showed both methods
performed very efficient for the classification and location of faults, and that ELM is faster
to train. Thanasak W. et al. [15] presented fault detection and identification (FDI) of mount
head damage in head gimbal assembly (HGA) manufacturing to improve slider loss defect
(SLD) using an artificial neural network (ANN). The mount heads were classified into
healthy and three fault types, which were dependent on damage level. The input feature is
an image of mount head each group. The performance on classifying each damage type
Systems 2023, 11, 268
4 of 33
was 94.3%, which was better than the result obtained from using voltage and vacuum
methods. Prathan C. et al. [16] developed fault detection and diagnosis techniques for
linear bearing in an auto core adhesion mounting (ACAM) machine using vibration signal
analysis and motor current analysis. The vibration signal was transformed to crest factor
and fast Fourier transform (FFT), while motor current was analyzed using analysis of
variance (ANOVA). The researchers were able to separate healthy linear motor from five
fault types with good classification accuracy. Prathan C. et al. [17] used ANN to segregate
healthy bearings from five defect types using the spectrum of FFT of vibration during
linear motor movement, motor current, and crest factor as input parameters. The results
showed an accuracy of up to 93% using the triple parameters. In subsequent work, Prathan
C. et al. [18] used ANN to classify anomalies in a high-speed auto core adhesion mounting
machine, followed by using a PI servo to compensate the machine back to an acceptable
condition. They achieved an accuracy of 100% in the classification via the confusion
matrix. They also demonstrated two types of gain scheduling, namely, discrete gain and
continuous gain, which resulted in 86% and 93% error reduction, respectively. Ahmed R.
Nasser et al. [19]. presented intelligent FD and identification on an analogue circuit using a
fuzzy-logic classifier. Testing was done on a low-pass filter analogue circuit. An average
98% F score accuracy was shown for this classifier. Yun Peng Zue et al. [20] presented
a study on the condition monitoring of a roto shaft system using frequency response
function (FRF). The researchers conducted two experimental studies to evaluate the system
using different methods. The first study focused on rub-impact faults, while the second
study examined misalignment. Three methods were compared in the analysis: the first
method utilized harmonic data, the second employed harmonic excitation, and the last
method utilized potential from experimental data. The results demonstrated that the second
method was better at detecting faults than the other methods. Naderi E. and Khorasani
K. [21] presented a data-driven approach for fault detection, isolation, and estimation of
aircraft gas turbine engine actuators and sensors. The frequency response of the gas turbine
engine was considered and estimated using Markov parameter estimation. These data
were then utilized for the direct design and implementation of fault detection, isolation,
and estimation filters. Samree J. et al. [22] presented transient analysis method using a
high-pass filter circuit in a high-voltage system. Faults were detected very quickly—within
1 ms. Abdelila E. [23] presented a fault-detection method for photovoltaic (PV) panels using
k-means clustering. The input data were collected from thermal images, and an image
processing technique was applied to convert the RGB color space to HSV. The number of
clusters, k, was determined using the elbow method and the average silhouette method.
The results demonstrated that the k-means algorithm successfully detected the location of
faults on PV panels. Christoph K. et al. [24] presented vibration classification on a motor
using machine learning. Several acceleration sensors were evaluated, and the vibration at
each motor speed was used as the label for the machine learning model. Smoothing and
differencing were used for feature extraction. Performance based on decision tree, gradient
boost, and focus cluster algorithms were presented in terms of accuracy and computing
time. A. M. Umbrajkaar et al. [25] presented vibration analysis on rotating machines using
a combined approach of SVM and ANN. The data from triaxial acceleration sensors were
collected in three groups: healthy, parallel misalignment, and angular misalignment. The
data were transformed using DWT before passing to the machine learning model.
Based on our review of various system identification techniques, we have concluded
that the WLS technique is appropriate for signal writing machines because it does not
re-quire an initial gradient search parameter. However, WLS requires an appropriate
weighting function. The Autsadayut [3] method works well under stable system conditions
and with linear viscoelastic materials. The Noda [5] method works well for electrical plant
identification when WLS reveals equidistant points in the frequency response. Either way,
signal writing machines represent a resonance mode with an inequal frequency band. Some
resonance modes tend to misidentify frequencies or require multiple iterations to meet
criteria. We used [3,5] the peak partition technique to model a signal writing machine.
Systems 2023, 11, 268
5 of 33
On several classification techniques, we decided to use an ANN classifier as our machine learning classifier. We used each machine group’s performance to represent the label
output. The operation parameters that show a high correlation to the dynamics behavior of
each transfer function will be used as feature inputs during the signal write process.
The paper is organized as follows: Section 2 describes the materials and methods
used to analyze and estimate mathematical models in transfer function form, including
time and frequency domains. The machine classification using an artificial neural network
(ANN) structure is also presented in this section. Section 3 describes the results of the
mathematical model obtained from system identification techniques and the results of
machine classification using ANN. Section 4 discusses the validation results for each model
structure and the accuracy of machine classification. Finally, Section 5 offers conclusions
and alludes to further work.
2. Materials and Methods
This section describes the operation of a spiral seed self-servo track writing machine for
closed-loop position control during HDD testing to collect various data from the machine
during testing. It also discusses the structure of a mathematical model in terms of a transfer
function, the procedure for mathematical model estimation using system identification
methods, and the analysis of the response of the test machine in the time and frequency
domain for each machine characteristic. Additionally, it demonstrates the classification
procedure of the ANN technique for testing machines.
2.1. Position Control
The spiral seed self-servo track writing machine, as shown in Figure 1, is used for
writing a reference signal in a spiral pattern. The machine consists of 16 signal write
slots, with one industrial computer controlling 4 signal write slots. It runs on multithread
software called WinSTW. Each slot consists of a drive figure clamping, an electric control
board, and a positioner known as MicroE. The positioner model is the PA2000, which
is integrated with a voice-coil motor (VCM) and an absolute encoder, and the amplifier
model is SA200. The maximum stroke range is 40 degrees, with an encoder resolution of
4.68 nanoradians. The rotary shaft of the MicroE is assembled with a customized positioner
arm that has an adjustable pushpin height to support multiple product specifications.
Figure 1. Spiral seed self-servo track writing machine.
The position control system shown in Figure 2 consists of a plant that is engaged with
complex loads such as hard disk drives (HDDs). The controller is a PID controller with
low-pass and notch filters.
Systems 2023, 11, 268
6 of 33
Figure 2. Position control system.
To obtain frequency response data, MicroE Systems Motion Control software version
6.3 was used, in frequency analysis mode. The input swept-sine signal was set to 0.4 Vpp,
with frequency range 40–5000 Hz. The sampling time was set to 39.6 usec for data collection.
Bode plots were generated by injecting the sine input signal at the junction point be-tween
the controller and amplifier. We then collected data from the encoder, where the position
changed, as illustrated in in Figure 3.
Figure 3. Frequency response function data collection.
Please note that the phase response was displayed as a phase margin and was wrapped
at positive(+)180 to negative(−)180 degrees.
The frequency response data were saved in
−
a bod file and imported into the system identification software. As shown in Figure 4,
the phase display data from the MicroE Systems Motion Control software needs to be
unwrapped and corrected as a phase response.
The correct phase response with unwrap version is shown in Figure 5, below.
2.2. Model Structure
The frequency response of the high-precision machine in Figure 5 can be used for
system identification. The model structure was defined by a combination of low- and
high-frequency models.
−
Systems 2023, 11, 268
7 of 33
Figure 4. Mechanical plus servo amplification frequency response from MicroE Systems
Control software.
Figure 5. Mechanical plus servo amplification frequency response.
The model structure for the low-frequency range can be obtained by using the response
data from the first response frequency up to 500 Hz. When we converted the frequency axis
from Hz to rad/sec, the magnitude response exhibited a slope of −40 dB/decade, with a
small gain increase around 50 to 70 Hz. These results motivated us to create a second-order
standard form model, as shown in (1), where ω0 and ζ 0 are the natural frequency and
damping ratio, respectively. Kdc represents the DC gain of the model, which combines the
plant gain and amplifier gain.
G (s)Second Order Std.Form =
Kdc ω02
s2 + 2ζ 0 ω0 s + ω02
𝐾
𝐺(𝑠)
.
=
−
(1)
𝜔
𝐾 𝜔
𝑠 + 2𝜁 𝜔 𝑠 + 𝜔
𝜁
Systems 2023, 11, 268
8 of 33
The model structure for the high-frequency range can be obtained by using the frequency response data above 500 Hz, represented by a ∏ type resonance model as shown in
(2), where ωnr and ζ nr are the natural frequencies and damping ratios of the anti-resonance
modes, respectively. ωdr and ζ dr are the natural frequencies and damping ratios of the
resonance modes, respectively.
N
G (s) Resonance =
2
s2 + 2ζ n ωnr s + ωnr
∏ 2
2
r =1 s + 2ζ r ωdr s + ωdr
(2)
The plant transfer function is shown in (3).
G (s) Plant = G (s)Second Order Std.Form × G (s) Resonance
(3)
2.3. System Identification Algorithm
To find model parameters and order in high-precision machines, OLS and WLS were
used, with an iterative weighting function. We considered using the general transfer
function G (s) in (4).
N (s)
(4)
G (s) =
D (s)
N (s) and D (s) represent the numerator and denominator of the transfer function,
respectively. They can be expressed in polynomial form, as shown in (5).
G (s) =
a0 + a1 ( s ) + a2 ( s )2 + . . . + a n ( s ) n
(5)
1 + b1 (s) + b2 (s)2 + . . . + bm (s)m
The error ε(s) of actual frequency response data F (s) and model G (s) are represented
in (6).
ε(s) = F (s) − G (s)
(6)
We replaced G (s) with the numerator and denominator.
ε(s) = F (s) −
N (s)
D (s)
(7)
Then, we multiplied both sides of (7) by D (s).
ε(s) D (s) = F (s) D (s) − N (s)
(8)
When the error related to the frequency response data and model was zero, the left
side of (8) was also zero. When we replaced N (s) and D (s) with the polynomials from (5),
the result was (9).
0 = F (s) 1 + b1 (s) + b2 (s)2 + . . . + bm(s)m
(9)
− a0 + a1 ( s ) + a2 ( s )2 + . . . + a n ( s ) n
Next, we rearranged (9) to (10).
F (s)
= a n ( s ) n + . . . + a2 ( s )2 + a1 ( s ) + a0
− F (s)bm (s)m − F (s)b2 (s)2 − . . . − F (s)b1 (s)
(10)
Then, we expressed (10) in matrix form as (11), to solve the least squares solution.
Ax = b
A is the coefficient matrix, as shown in (12).
s0k − F (sk )sm
A = snk . . . s1k
k
(11)
...
− F (sk )sk
(12)
Systems 2023, 11, 268
9 of 33
Since s0k = 1, Equation (12) is modified to (13),
1 − F (sk )sm
A = snk . . . s1k
k
...
− F (sk )sk
(13)
x is a coefficient vector of the numerator and denominator, as shown in (14).
x = an
...
a1
a0
bm
...
b1
T
(14)
b is the actual frequency response vector, as shown in (15)
b = F ( s )1
F ( s )2
F ( s )3
...
F (s)k
T
(15)
m and n represent the number of poles and zeros, respectively. k is the measurement point.
The OLS solution is shown in (16).
−1
AT b
(16)
−1
A T Wb
(17)
x = ( AT A)
The WLS solution is shown in (17).
x = ( A T WA)
W is the weighting function [5], which was changed adaptively according to the
solution of least squares in the previous step. j is the iteration step and the initial W is the
identity matrix, which is the OLS solution.
W = (W ) ( j − 1 ) × ε ( j − 1 )
(18)
To evaluate model accuracy, we had to consider the error of transfer function in (7)
from each data point. Then, the norm value [3] was calculated, as shown in (19).
kεk =
q
| ε 1 |2 + | ε 2 |2 + . . . + | ε k |2
(19)
2.4. Second-Order Standard Form Model Estimation
The second-order standard form model parameters in (1) are identified by (16) and
(17). The model parameters are shown in (20). The fitting result is illustrated in Figure 6.
The magnitude response of the model fit well to the actual magnitude response at low
frequency when the sample point of phase shift was not less than −180 degrees. The
minimum model error of model was 808. The best weighting function is shown in Figure 7.
The weighting function showed gain reduction at high frequency.
G (s)Second Order Std.Form =
75 × 3692
s2 + 2 × 0.22 × 369s + 3692
(20)
−−
Systems 2023, 11, 268
𝐺(𝑠)
𝐺(𝑠)
.
.
75
75×
×369
369
=
= 𝑠 + 2 × 0.22 × 369𝑠 + 369
𝑠 + 2 × 0.22 × 369𝑠 + 369
10 of 33
Figure 6. Second-order standard model.
Figure 7. Weighting function of second-order standard model.
There were delays related to phase response at frequencies exceeding 1 kHz. Thus, we
added a transportation delay to (1). We multiplied delay up to second-order standard form
to improve phase errors at high frequency, while maintaining magnitude response. The
model structure was modified as shown in (21).
G (s)Second Order Std.Form = e(−td ×s) ×
Kdc ω02
s2 + 2ζ 0 ω0 s + ω02
(21)
MicroE Systems Motion Control software was used in Time Response Analysis mode.
The output of the time response data was exported to a plot, as shown in Figure 8, with
the delay indicating the time between the two samples. We estimated time delays from
39.6 usec to 79.4 usec. Based on a Golden Section Method search (GSM), the time delay was
about 65 usec. The model response and actual response are shown in Figure 9. The model
Systems 2023, 11, 268
𝐺(𝑠)
𝐺(𝑠)
. .
==
𝑒 (𝑒
𝐾 𝐾𝜔 𝜔
× × 𝑠 + 2𝜁 𝜔 𝑠 + 𝜔
𝑠 + 2𝜁 𝜔 𝑠 + 𝜔
( × ×
) )
11 of 33
parameters are illustrated in (22). The minimum model error of second-order standard
form plus delay time was reduced to 277.
G (s)Second Order Std.Form w_Delay = e(−65×10
−6 s )
×
s2
75 × 3692
+ 2 × 0.22 × 369s + 3692
(22)
Figure 8. Time response.
Figure 9. Second-order standard form plus delay time.
2.5. Resonance Model Estimation
The resonance model in (2) was separately identified by using a resonance peak
command (findpeaks in MATLAB® ). The first data point from the first group was selected
from the end of the second-order standard form data.
Then, the first resonance frequency and center frequency from each group were
partitioned and indicated in dash line, as shown in Figure 10.
Systems 2023, 11, 268
12 of 33
Figure 10. Resonance peak grouping.
Next, (16) and (17) were used to estimate the model parameter, under the following
specific conditions:
The higher magnitude peak was identified first.
The maximum order in each of part was limited to the number of peaks, multiplied
by 2.
One peak equates 2 orders.
The minimum order from each part was set as a second-order.
The result of pole and zero must be in the left half-plane of the s-domain.
The result of pole and zero must be within a fitting frequency range.
The result must be a complex conjugate of pole and zero. This condition was set to
maintain phase response from a roughly phased slope of the second-order standard form
plus delay time model, from (21).
The last was around 4000 Hz. The fitting sequence for identifying the resonance
groups followed a higher-to-lower peak magnitude, as indicated by the numbers at the top
of each peak in Figure 11. The first resonance group identified was around 1200 Hz, while
the last resonance group was around 4000 Hz.
Figure 11. Resonance fitting sequence.
Systems 2023, 11, 268
13 of 33
The fitting result of the resonance model is shown in Figure 12. The model parameters
are illustrated in (23) to (28). Several resonance peaks can be identified as having a secondorder complex conjugate of pole and zero. The phase response was maintained at 0 degrees.
Figure 12. Resonance fitting result.
Resonance model group 1.
𝑠 2+ 566𝑠 + 5.6𝑒 7
𝐺𝑓
= s + 566s + 5.6e
G f(𝑠)
1 ( s ) =𝑠 2+ 570𝑠 + 5.9𝑒 7
s + 570s + 5.9e
(23)
Resonance model group 2.
𝑠 + 197𝑠 + 1.3𝑒
𝐺𝑓 (𝑠) = s2 + 197s + 1.3e8
G f 2 (s) =𝑠 2+ 221𝑠 + 1.3𝑒 7
s + 221s + 1.3e
(24)
Resonance model group 3.
G f 3 (s) =
s6 + 3876s5 + 9.7e8 s4 + 2.2e12 s3 + 3.0e17 s2 + 3.0e20 s + 3.1e25
s6 + 3561s5 + 9.8e8 s4 + 2.0e12 s3 + 3.1e13 s2 + 2.8e20 s + 3.2e25
(25)
Resonance model group 4.
s2 + 1194s + 5.0e8
s2 + 1005s + 5.0e8
(26)
s4 + 1744s3 + 1.8e9 s2 + 1.5e12 s + 7.9e17
s4 + 1366s3 + 1.7e9 s2 + 1.2e12 s + 7.5e17
(27)
G f 4 (s) =
Resonance model group 5.
G f 5 (s) =
Resonance model group 6.
G f 6 (s) =
s2 + 921s + 6.3e8
s2 + 898s + 6.4e8
(28)
The combination of second-order standard form plus delay time with a resonance
model is shown in Figure 13. The model error was 114.
𝐺𝑓 (𝑠) =
Systems 2023, 11, 268
𝑠 + 921𝑠 + 6.3𝑒
𝑠 + 898𝑠 + 6.4𝑒
14 of 33
Figure 13. Combination of second-order standard form plus delay with resonance model.
2.6. Resonance Gain Adjustment
The result of the second-order standard form plus delay time combination and the
resonance model in Figure 13 reveals a small gain, which decreased at high frequency.
We improved the model by maintaining the phase response and fine-tuning the small
resonance gain. (2) was modified to (29), where kres was the resonance gain. Because onedimensional gain searching was required, we used the GSM method. We set the minimum
𝑘
and maximum resonance gain at 0.1 dB and 2 dB, respectively.
The best resonance gain was
adjusted to 1.5 dB. The model error was reduced to 93. The comparison of actual response
and model response is shown in Figure 14.
N
2
s2 + 2ζ n ωnr s + ωnr
2
2
r =1 s + 2ζ r ωdr s + ωdr
G (s) Resonance = kres × ∏
(29)
𝑠 +
2𝜁at𝜔
𝑠 + 𝜔order, with a
The final model parameter is shown in (30). The model
was
the 22nd
𝐺(𝑠)
=𝑘 ×
small delay time.
𝑠 + 2𝜁 𝜔 𝑠 + 𝜔
G (s) = e(−65×10
−6 s )
×
×1.5
75×3692
s2 +2×0.22×369s+3692
s18 + 8549s17 + 4.1e16 s16
+ 3.0e13 s15 + 7.0e18 s14
+4.4e22 s13 + 6.6e27 s12 + 3.4e31 s11 + 3.8e36 s10
+1.6e40 s9 + 1.3e45 s8 + 4.2e48 s7 + 2.9e53 s6 + 6.3e56 s5
+3.6e61 s4 +4.9e64 s3 +2.3e69 s2 +1.5e72 s+5.4e76
s18 +7578s17 +4.1e9 s16 +2.7e13 s15 +6.9e18 s14
+3.9e22 s13 + 6.6e27 s12 + 3.1e31 s11 + 3.7e36 s10
+1.4e40 s9 + 1.3e45 s8 + 3.8e48 s7 + 2.9e53 s6 + 5.9e56 s5
+3.6e61 s4 + 4.7e63 s3 + 2.2e69 s2 + 1.4e72 s + 5.8e76
(30)
Systems 2023, 11, 268
15 of 33
Figure 14. Resonance gain adjustment.
As the signal writing machine utilizes a digital PID controller, it is necessary to convert
the plant model in (30) from continuous time to discrete time. This conversion can be
achieved using the C2D (continuous to discrete) method. In this case, we employ the zeroorder hold (ZOH) method to transform the model. The resulting transformed equation is
75 × 369
provided in (31).
×
𝐺 (𝑠) = 𝑒
×
𝑠 + 22 × 0.22 × 369𝑠 + 369
𝑠 + 3.0𝑒 𝑠 + 7.0𝑒 𝑠
18
17
16 3.4𝑒
+4.4𝑒 z 𝑠− 17.4z
+ 6.6𝑒+ 103.8z
𝑠 +
𝑠15 ++1192z
3.8𝑒14 𝑠
− 409.4z
10 𝑠 + 6.3𝑒 𝑠
+1.6𝑒 −2713z
𝑠 +131.3𝑒
𝑠 12+−4.2𝑒
2.9𝑒
+ 4992z
7584z11𝑠 ++
9634z
4 z94.9𝑒
6 − 2355z
𝑠 +
𝑠8 −
+7067z
2.3𝑒7 +𝑠4490z
+ 1.5𝑒
𝑠 +5 5.4𝑒
+ 9300z
× 1.5 +3.6𝑒 −1.0e
(31)
4
3
2
𝑠×1.5+ 7578𝑠
+
4.1𝑒
𝑠
+
2.7𝑒
𝑠
+
6.9𝑒
𝑠
+998.7z −331z +80.9z −13.1z+1.1
18
17
16
15
14
+71.6z −𝑠283.3z
+827z 𝑠
+3.9𝑒 z 𝑠−11.9z+136.6𝑒
+ 3.1𝑒
+ 3.7𝑒 𝑠
−1886z + 3478z12 − 5292z11 + 7205z10
+1.4𝑒 𝑠 +91.3𝑒 𝑠8 + 3.8𝑒7 𝑠 + 2.9𝑒
𝑠 + 5.9𝑒 𝑠
−7205z + 6506z − 4945z + 3142z6 − 1648z5
+3.6𝑒 𝑠 + 44.7𝑒 𝑠 3+ 2.2𝑒 2 𝑠 + 1.4𝑒 𝑠 + 5.8𝑒
+0.01165z+0.003311
G (z) = z−2 × 0.001013z
𝑠 z+
2 −8549𝑠
1.944z+0.9337+ 4.1𝑒
+698.3z − 231.1z + 56.5z − 9.14z + 0.7
The flowchart of the system identification method, including the adaptive weight
least squares and peak detector for the partition resonance frequency model, is shown in
Figure 15.
2.7. Frequency and Time Response Analysis
In the manufacturing
process,
even though
same type of signal writing machine is
0.001013𝑧
+ 0.01165𝑧
+the
0.003311
𝐺 (𝑧)differences
= 𝑧 × in its characteristics can result in varying dynamic response performance.
used,
𝑧 − 1.944𝑧 + 0.9337
Figures 16–20 show the comparison of the mechanical FRF and PES of three classes of
𝑧 − 17.4 𝑧 + 103.8 𝑧 − 409.4 𝑧 + 1192 𝑧
signal writing machines, with class 0 being defined as a healthy machine that produces a
−2713
+ 4992 𝑧 − 7584 𝑧 + 9634 𝑧
good HDD during the signal writing𝑧operation.
−1.0𝑒 𝑧 + 9300 𝑧 − 7067 𝑧 + 4490 𝑧 − 2355 𝑧
+998.7
𝑧 − 331 𝑧 + 80.9 𝑧 − 13.1 𝑧 + 1.1
× 1.5
𝑧 − 11.9 𝑧 + 71.6 𝑧 − 283.3 𝑧 + 827 𝑧
−1886 𝑧 + 3478 𝑧 − 5292 𝑧 + 7205 𝑧
−7205 𝑧 + 6506 𝑧 − 4945 𝑧 + 3142 𝑧 − 1648 𝑧
+698.3 𝑧 − 231.1 𝑧 + 56.5 𝑧 − 9.14 𝑧 + 0.7
Systems 2023, 11, 268
16 of 33
Figure 15. Flowchart of system identification.
(a)
(b)
Figure 16. Comparison of (a) FRF; (b) PES of machine class 1, sensor degradation.
Systems 2023, 11, 268
17 of 33
Figure 17. Sensor voltage comparison: (a) healthy machine; (b) sensor degradation.
(a)
(b)
Figure 18. Comparison of (a) FRF and (b) PES of machine class 2, pushpin degradation.
Figure 18. Comparison of (a) FRF and (b) PES of machine class 2, pushpin degradation.
(a)
(b)
Figure 19. Pushpin comparison: (a) healthy machine; (b) pushpin degradation.
Systems 2023, 11, 268
18 of 33
(a)
(b)
Figure 20. Comparison of (a) FRF and (b) PES of machine class 3, Tunable machine.
Class 1 in Figure 16a shows DC gain increase from 75 to 83 in magnitude plot compared
to the healthy machine. The (32) and (33) show second-order standard form model of a
healthy model and class 1 model, respectively. Figure 16b PES shows high variations across
HSA seek from OD to ID as indicated by the dash line. According to repair history, the
sensor voltage drops from more than 1.5 V to 0.5 V, as shown in the red line in Figure 17.
This machine is categorized as having a sensor fault and cannot be compensated due to the
actual degradation of the sensor. These machines will be shut down for repair.
75 × 369
75 × 3692
=
𝐺(𝑠)
.
G (s)Second Order Std.Form
= 2𝑠 + 2 × 0.22 × 369𝑠 + 369
s + 2 × 0.22 × 369s + 3692
83 × 426
𝐺(𝑠)
=
𝑠 +2×
× 2426𝑠 + 426
830.25
× 426
G (s)sensor f ault = 2
s + 2 × 0.25 × 426s + 4262
(32)
(33)
Class 2 in Figure 18a shows that the first resonance frequency has shifted from around
1230 Hz to 665 Hz, resulting in an increased first resonance magnitude. The natural
frequency of the second-order standard form model has moved to around 110 Hz, which
is close to the HDD spindle motor speed at 120 Hz (7200 rpm). The (34) and (35) show
first resonance model of a healthy model and class 2 machine model, respectively. The
PES in Figure 18b exhibits slow response and large oscillations when seeking starts from
the OD as indicated by the dash line. DFT is calculated using PES at a servo sample from
200 to 400, which is in the OD zone, and the DFT control band is from 600 Hz to 1500 Hz.
This causes the DFT to exceed the control limit. The machine has a loose pushpin, which
is a mechanical degradation and cannot be compensated. It will be shut down for repair.
Figure 19 shows a comparison of a healthy pushpin and a loose pushpin.
𝑠 + 566𝑠 + 5.6𝑒
𝐺𝑓 (𝑠)
=2
7
570𝑠
+ 5.9𝑒
s 𝑠+ +
566s
+ 5.6e
G f 1 (s)healthy = 2
(34)
s + 𝑠570s
+ 5.9e+7 1.6𝑒
+ 280𝑠
𝐺𝑓 (𝑠)
=
𝑠 + 266𝑠 + 1.7𝑒
s2 + 280s + 1.6e7
G f 1 (s) push pin loose = 2
(35)
s + 266s + 1.7e7
Class 3 in Figure 20a has a comparable FRF, with the first resonance frequency slightly
higher than that of the healthy machine. Upon closer examination in Figure 21, the gain
crossover frequency is slightly lower than the healthy machine as indicated by the arrow
dash line, and the phase margin tends to be higher, leading to low closed-loop bandwidth
and a larger damping ratio. When using the same controller gain, the response of this class
becomes slower than the healthy machine. The (36) and (37) show first resonance model
Systems 2023, 11, 268
19 of 33
of the healthy model and class 3 model, respectively. The PES in Figure 20b shows slow
response when compared to the healthy machine as indicated by the dash line. This class
exhibits a slight performance deviation and can be optimized by adjusting the controller
gain to meet frequency domain stability criteria and reduce
PES. + 5.6𝑒
𝑠 + 566𝑠
𝐺𝑓 (𝑠)
G f 1 (s)healthy
𝐺𝑓 (𝑠)
G f 1 (s)tunable
=
7
+ 570𝑠
s2 + 566s𝑠+ 5.6e
= 2
s + 570s + 5.9e7
+ 5.9𝑒
𝑠 + 575𝑠 + 4.7𝑒
+ 4.9𝑒
= + 4.7e7
s2 + 575s
𝑠 + 631𝑠
= 2
s + 631s + 4.9e7
(36)
(37)
Figure 21. Comparison of zero dB gain cross frequency of machine class 3.
The DFT results of each machine class have been compared with those of the healthy
machine, as shown in Figure 22. The healthy class shows minimum amplitude for each
harmonic. In the sensor fault class (red line), a high amplitude is observed around 1100
Hz, which is considered system noise since it is beyond the closed-loop design bandwidth.
In the pushpin-loose class (yellow line), a high amplitude is observed around 700 Hz,
which is considered the effect of the first resonance. In the tunable class (green marker),
a high amplitude is observed at low frequencies, which is considered a result of the slow
system response.
2.8. Tester Classification: ANN Structure
From the deviation among the three machine classes, in order to make the best decisions for each machine class, it is necessary to separate them based on their deviations.
Class 1 and class 2 are classified as machines that require maintenance, while class 3 is
identified as a tunable machine. To address this issue, the authors propose a machine
learning method that classifies machines based on their characteristics.
Systems 2023, 11, 268
20 of 33
Figure 22. DFT comparison.
The authors utilize an artificial neural network (ANN) structure from the Keras python
package to classify machines into different classes. The ANN structure comprises 10 input nodes, representing the selected features. To achieve high accuracy, the authors initially experimented with a large number of nodes in each layer, using a configuration of
200–500–200 nodes. The rectified linear unit (ReLU) activation function is applied in these
hidden layers.
The output layer of the ANN consists of 4 nodes, representing the 4 machine classes.
The SoftMax activation function is utilized in the output layer to generate class probabilities.
Figure 23 illustrates the structure of the ANN.
Figure 23. ANN structure.
During training, the network is trained for 200 epochs with a batch size of 10. The
authors employ the Adam optimizer with a learning rate of 0.01. The categorical crossentropy loss function is used to evaluate the network’s performance.
The ReLU activation function is defined by the Formula (38), while the SoftMax
activation function is defined by the Formula (39). These activation functions play a crucial
role in introducing non-linearity and producing appropriate output representations for
classification.
ReLU (z) = max (0, z)
(38)
𝑅𝑒𝐿𝑈(𝑧) = 𝑚𝑎𝑥(0, 𝑧)
𝑆𝑜𝑓𝑡𝑀𝑎𝑥(𝑧) =
∑
𝑒
𝑒
Systems 2023, 11, 268
21 of 33
So f tMax (z)i =
e zi
zj
∑N
j =1 e
(39)
2.9. Data Collection and Preparation
In machine learning, data collection and preparation are crucial. This study collected
over 200 features from 4 classes of signal writing machines. The raw data for 486 drives
were obtained from the manufacturing log file. Class 0, representing the healthy condition,
has 198 drives in its dataset. Class 1, representing a sensor fault machine, has 106 drives.
Class 2, representing a pushpin-loose machine, has 78 drives. Finally, class 3, representing
a tunable machine, has 104 drives.
The data was first randomly indexed, and then cleaned using the modified version of
the FeatureSelector Python package [26] to support multiclass data. The cleaning process
involved removing missing data that accounted for more than 20% of each feature, as
well as single unique data and linear correlations that showed a 90% or higher correlation
to other features. The feature ranking was determined using light GBM, and Figure 24
displays the 10 most important operating parameters that correlate highly with the label.
These parameters are considered input features.
Figure 24. Ten most important operating parameters.
DFT in frequency domain is calculated from raw PES data in time domain by (40).
Only DFT amplitude was considered (41):
N −1
Xk =
2πkn/N ) − i.sin
)]
𝑋 ∑
= xn [cos
𝑥 (𝑐𝑜𝑠(2𝜋𝑘𝑛/𝑁)
− 𝑖.(2πkn/N
𝑠𝑖𝑛(2𝜋𝑘𝑛/𝑁)
(40)
n =0
X𝑋k is DFT which include information of both amplitude and phase
X𝑋n is PES value in time domain
N is number of samples
𝑁
n is current sample
𝑛
k is current frequency
𝑘
q
real𝑟𝑒𝑎𝑙(𝑋
( Xk )2 +)imaginary
( Xk )2 )
+ 𝑖𝑚𝑎𝑔𝑖𝑛𝑎𝑟𝑦(𝑋
amp𝑎𝑚𝑝
= =
N 𝑁
(41)
The 10 histograms and boxplots of 10 feature are shown in Figures 25–34.
1.
DFT_Euclidean_Low: Represents the Euclidean norm (38) of DFT harmonics at low
frequencies between 600 and 1500 Hz. The symptoms associated with this feature are
Systems 2023, 11, 268
22 of 33
observed in the analysis of frequency and time response in Section 2.7. A high DFT
value can be attributed to three symptoms:
•
•
•
(a)
High noise, as shown in Figure 16, where the model exhibits an overall increase in
magnitude, including noise.
High variation during the start of movement, as shown in Figure 18, where the model’s
= 𝑎𝑚𝑝 + ⋯ + 𝑎𝑚𝑝
first resonance shifts to 𝐷𝐹𝑇
a lower frequency.
𝐷𝐹𝑇
= 𝑎𝑚𝑝resulting
+ ⋯from
+ 𝑎𝑚𝑝
Slow response due to a low bandwidth
the reduction of complex
pole/zero cancelation, as shown in Figure 20.
(b)
Figure 25. (a) Histogram and (b) boxplot of feature DFT_Euclidean_Low.
(a)
(b)
Figure 26. (a) Histogram an (b) boxplot of feature DFT_Euclidean_High.
Systems 2023, 11, 268
23 of 33
(a)
(b)
Figure 27. (a) Histogram and (b) boxplot of feature Harmonic3.
(a)
(b)
Figure 28. (a) Histogram and (b) boxplot of feature PES_Range.
(a)
(b)
Figure 29. (a) Histogram and (b) boxplot of feature FF_Factor.
Systems 2023, 11, 268
24 of 33
(a)
(b)
Figure 30. (a) Histogram and (b) boxplot of feature PES_STD.
Figure 30. (a) Histogram and (b) boxplot of feature PES_STD.
(a)
(b)
Figure 31. (a) Histogram and (b) boxplot of feature Harmonic24.
(a)
(b)
Figure 32. (a) Histogram and (b) boxplot of feature Harmonic1.
Systems 2023, 11, 268
25 of 33
(a)
(b)
Figure 33. (a) Histogram; (b) Boxplot of feature Harmonic6.
(a)
(b)
Figure 34. (a) Histogram and (b) boxplot of feature Harmonic8.
The histogram in Figure 25a indicates significant variation in the sensor fault class,
while in Figure 25b, the median value of the healthy class is significantly lower than that
of the other classes. Outlier data points are marked with diamond symbols. This feature
plays a crucial role in distinguishing the healthy class from the other classes.
DFTLow =
2.
amp2600 + · · · + amp21500
(42)
DFT_Euclidean_High: Represents the Euclidean norm (43) of DFT harmonics at high
frequencies between 1500 and 3000 Hz. This feature is also important for distinguishing the healthy class from
the other
of this feature is similar
𝐷𝐹𝑇
= classes.
𝑎𝑚𝑝 The
+ interpretation
⋯ + 𝑎𝑚𝑝
to the first feature, as discussed in Section 2.7, but with a focus on higher frequencies
of interest. The histogram in Figure 26a indicates significant variation in the sensor
fault class, while in Figure 26b, the median value of the healthy class is significantly
lower than that of the other classes.
DFTHi =
3.
q
q
amp21500 + · · · + amp23000
(43)
Harmonic3: Represents the amplitude of the third DFT harmonic (37). Due to the
nature of a low-bandwidth system, a slow response symptom will result in a high DFT
magnitude at low frequencies, as shown in Figure 22. The histogram in Figure 27a
demonstrates narrow variation in the healthy, sensor fault, and pushpin-loose classes.
Additionally, the boxplot in Figure 27b reveals that the median of the tunable class is
higher than that of the other classes. This indicates that the Harmonic3 feature has
high potential for distinguishing the tunable class from the other classes.
Systems 2023, 11, 268
26 of 33
4.
5.
6.
7.
8.
9.
10.
PES_Range: Represents the range of PES peaks in track units from OD to ID. A high
value indicates a high spread of position error in the system or high overshoot, which
can be attributed to a low closed-loop damping ratio when the plant described in
(30) is operated as a closed-loop system. However, it should be noted that some unit
conversions and truncations of floating-point numbers in reporting may affect the
accuracy. In Figure 28a,b, there is no obvious symptom to use this feature alone for
separating machine classes. While this feature alone may not be sufficient, it can
contribute to increasing accuracy when combined with other features.
FF_Factor: Represents the feed-forward gain that is adjusted in the feed-forward
controller to minimize PES during mechanical HSA warmup. If the system exhibits a
smooth slope of a second-order system at low frequency, we can add an appropriate
amount of FF factor to the feed forward filter in order to reduce PES. However, if the
effect of the first resonance impact at low frequency is present, we cannot increase the
FF gain further to reduce PES. In the histogram shown in Figure 29a, the pushpin-loose
class exhibits a significantly lower median than the other classes, which is attributed
to the effect of the first resonance shifting to a lower frequency.
PES_STD: Represents the standard deviation (44) of PES in the DFT calculation within
a certain range. A high value of PES_STD can be interpreted as a high settling time,
indicating a low closed-loop system bandwidth. In Figure 30a, the data from the
tunable and some sensor fault classes are separated from the healthy and pushpinloose classes. Additionally, in Figure 30b, the boxplot demonstrates high potential to
distinguish the tunable class from the other classes.
q
PES_STD = mean( abs( PES − PES.mean)2 )
(44)
Harmonic15: Represents the amplitude of the 15th DFT harmonic (41). The characteristic of a low closed-loop bandwidth system is a slow response, resulting in a high
DFT magnitude at low frequencies, as shown in Figure 22. In Figure 31a, a slightly
higher variation can be observed in the sensor fault class. However, in Figure 31b, the
boxplot demonstrates a high median value of the 15th harmonic for the tunable class,
indicating its potential to be separated from the other classes.
Harmonic1: Represents the amplitude of the first DFT harmonic (41), similar to other
low DFT harmonic features. The characteristic of a low closed-loop bandwidth system
is a slow response. In Figure 32a, the data of the tunable class is separated from the
other classes. In Figure 32b, the boxplot shows that the median of the tunable class
is noticeably distinct from the other classes. This indicates high potential to use this
feature for separating machine classes.
Harmonic6: Represents the amplitude of the sixth DFT harmonic (41). Although there
is no clear variation observed in the histogram shown in Figure 33a, the boxplot in
Figure 33b suggests that the data can be separated into two distinct groups. The first
group consists of the healthy and sensor fault classes, while the second group includes
the pushpin-loose and tunable classes.
Harmonic8: Represents the amplitude of the eighth DFT harmonic (41). In Figure 22,
high peaks can be observed around 1100 Hz and 2200 Hz, indicating the presence of
resonance and an increase in noise. These peaks correspond to a frequency response
where the gain increases across the entire frequency range. Figure 34a shows significant variation within the sensor fault class, suggesting the occurrence of faults.
Figure 34b highlights the high magnitude of the eighth harmonic, indicating its
potential to distinguish sensor fault classes from other groups.
Due to the 10 operating parameters being in different scales and units, the data were
normalized using (45) and (46). This normalization process ensures that the data have
a zero mean and a unit L2 norm. The normalization was performed using the normalize
Python package. Additionally, the output labels were encoded using the one-hot format, a
commonly employed technique for categorical variables in machine learning.
Systems 2023, 11, 268
27 of 33
To assess the model’s performance, the dataset was divided into three groups using
the Python train_test_split package. The split ratio used allocated 20% of the data for
testing, while the remaining 80% was further divided, with 20% for validation.
Due to the nature of the manufacturing process, the dataset contains fewer samples
from defective machines compared to healthy machines. To address this class imbalance,
the “stratify = y” argument was utilized during the data-splitting process. This ensured
that the data was partitioned in a manner that maintained an equal ratio of samples from
each class across the three groups (train, validate, and test). By stratifying the data based
on the output labels, the resulting datasets had a representative distribution of samples
from each class, facilitating more accurate model training and evaluation.
The final distribution of the data after splitting was as follows:
Train: 60% (291 drives);
Validate: 20% (97 drives);
Test: 20% (98 drives);
√
kxk = xT x
(45)
xnormalize =
xi
kxk
(46)
3. Results
This section describes the experiment result of system identification each model
structure and performance of machine classification by ANN.
3.1. System Identification Result
The signal writing machine is modeled, and the norm error is shown in Table 1. All
models are in a stable condition. The second-order standard form plus delay time model
showed a norm error of 808, but not all resonance peaks were taken into account. The final
model, which used the combined second-order standard form and resonance model under
stable conditions, showed a minimum error of 93. This was a good candidate model for
further controller optimization in tunable machine class because it was able to identify the
most important peaks.
Table 1. Model summary.
Model Type
Second-order standard form
Second-order standard form plus delay time
Second-order standard form plus delay time and resonance
Second-order standard form plus delay time and resonance
with gain (final)
Order
Norm Error
2
2
20
808
206
114
20
93
3.2. Tester Classification Result
Table 2 summarizes the time and frequency analysis of each machine class. Classes
1 and 2 will require repair. Class 3 will require optimization of the controller using the
mathematical model obtained from the system identification process.
Systems 2023, 11, 268
28 of 33
Table 2. Machine class analysis.
Machine Class
Class 0 Healthy
Class 1 Sensor fault
Class 2 Pushpin loose
Class 3 Tunable
Frequency
Represented by blue line.
Open-loop frequency
meets stability criteria
The magnitude response is
higher than that of the
normal machine
The first resonance
frequency shifts from 1230
Hz to around 665 Hz
The gain crossover
frequency and phase
margin do not meet
stability criteria
Time
Action
Small PES at measurement
area OD to ID
No action
High PES oscillation across
the head stack seeks from
OD to ID
In the OD zone, the PES
oscillation is bigger than
that of the healthy machine
The PES decreases more
slowly than that of the
healthy machine
Repair
Repair
Optimize controller
The confusion matrices for the training, validation, testing and all datasets are shown
in Figure 35. The diagonal elements of the confusion matrices represent successful classifications, while the off-diagonal elements represent unsuccessful ones. The numbers in the
matrices correspond to the data used in the training, validation, and testing processes. It
can be observed that the detection of all machine classes has an accuracy of 100%.
Figure 35. ANN confusion matrix.
To ensure the model’s performance, k-fold cross-validation was performed with 5 folds.
The training and evaluation dataset ratio is 1/k fold which is ratio 80:20: 80% for training
and 20% for evaluation. The accuracies of the validation dataset observed are 100% and
above 98%, as shown in Table 3. The cross-validation results indicate that the model is
reliable even if another train–test dataset is chosen. The loss function in each epoch is
showed in Figure 36.
Table 3. K-fold cross-validation.
Fold
1
2
3
4
5
Minimum Accuracy
Maximum Accuracy
Minimum Loss
Loss
×− 10−3
×− 10−8
0.8
3.4
5.4 ×− 10−8
8.0 ×− 10−6
51 × 10−3
−
3.4 × 10−8
−
Accuracy (%)
100
100
100
100
98.7
98.7
100
−
Systems 2023, 11, 268
29 of 33
−
Figure 36. Model loss.
The high accuracy shown in Table 3 indicates that it is possible to reduce the number
of features and neural sizes. We selected the top five features in Figure 24 and performed
hyperparameter tuning by reducing the number of nodes in the middle-hidden layer using
the keras_tuner python package. The range of node in the hidden layer was set from 32 to
512 nodes with a step size of 32. The model loss in Figure 36 remained flat after 100 epochs,
so we reduced the number of epochs from 200 to 150 to speed up the training process.
The optimization results yielded 384 nodes as shown in Figure 37 with 100% accuracy as
shown in Figure 38. This allowed us to reduce the number of parameters in the ANN
from 203,704 to 156,188 parameters. This resulted in an increase in prediction speed.
Figure 37. Loss function.
Systems 2023, 11, 268
30 of 33
Figure 38. Overall confusion matrix of tester classification with hyperparameter tuning.
The performance of the artificial neural network (ANN) in Figure 38 demonstrates high
accuracy. However, in order to explore models with lower complexity, we decided to select
the linear discriminant model. Using the same five most important features and the same
split ratio, the results are shown in Figure 39. Although misclassifications are observed
in the healthy, sensor fault, and tunable classes, the pushpin-loose class exhibits higher
accuracy. This model can serve as a baseline for further improvements in future models.
Figure 39. Linear discriminant analysis.
4. Discussion
A mathematical model of a signal writing machine has been successfully developed in a stable condition. The model consists of a second-order system with a small
delay connected in series with a resonance model. The parameters of the model were
determined using the adaptive weight least squares method and a peak detection for
resonance partition.
The first resonance peak was accurately identified as it is more critical and has a
significant impact on the system’s closed-loop response. If the first resonance frequency
falls within the DFT checking band, there is a high risk of failing to meet the DFT limit.
The model that exhibits minimal deviation from the baseline can be utilized for future
controller optimization to minimize PES and ensure stability.
However, due to the complexity of system identification and the time it takes, and
in some cases, mechanical degradation that cannot be compensated for by adjusting the
controller, machine classification is proposed.
The features of the machine learning model were derived from the DFT values in each
harmonic, which correlate with the resonance model obtained from system identification.
Systems 2023, 11, 268
31 of 33
The performance of the ANN classifier was promising, with 100% accuracy. The
cross-validation shows minimum accuracy at 98% with low variation, indicating its ability
to effectively distinguish between healthy and other classes of signal writing machines.
Hyperparameter tuning can be employed to reduce the neural size while maintaining
high accuracy.
The tunable class has significant potential for optimizing controller gain, wherein
sensor fault and pushpin-loose classes can be used for maintenance planning.
5. Conclusions
This research presents a classification system for machine performance using system
identification analysis and an artificial neural network (ANN). The machines are categorized into four classes base on symptom position error and defect in machine component:
0 (healthy machine), 1 (sensor fault), 2 (loose pushpin), and 3 (tunable machine). The transfer function model parameters are analyzed to determine their correlation with position
error response in both the time and frequency domains. Parameters from both domains
that show high correlation are used as feature inputs for the ANN, while the machine class
serves as the label.
The position error signal (PES) in the time domain is transformed to the frequency
domain using DFT at a specific OD seek location. Class 0 represents a healthy machine that
exhibits a small PES and no abnormalities in its components. All DFT harmonics are within
the pass criterion.
Symptoms of class 1 indicate low voltage in the encoder sensor. The PES shows high
variation across the stroke seek OD to ID, indicating that noise affects the response. The
system identification shows that the system gain increases, leading to increased noise
signals. The DFT shows several high harmonic peaks, resulting in HDD failing the DFT
criterion. This machine class requires repair by changing the positioner, which is unable
to compensate.
Symptoms of class 2 indicate a loose pushpin. The PES shows oscillation at the
OD zone. The system identification shows that the first resonance model changes from
1230 Hz to 665 Hz, which is within the closed-loop bandwidth, leading to an increase in
first resonance response. The DFT shows a high harmonic at the first resonance frequency,
resulting in HDD failing the DFT criterion. This machine class requires repair by changing
the pushpin, which is unable to compensate.
Class 3 PES exhibits smoother symptoms than classes 1 and 2. However, the PES
shows a slower decrease than a healthy machine. The system identification shows that the
first resonance model is closed to the healthy class. However, the gain cross-frequency is
lower than that of a healthy machine, and the phase margin tends to be larger than that of
a healthy machine, leading to a slower system response than the healthy class. The DFT
shows a high magnitude at the first harmonic, resulting in HDD failing the DFT criterion.
There are no obvious component parts degrading. Through rough experimentation, the
response of this machine class can be improved by optimizing the controller gain.
The DFT and PES from specific frequency ranges that correlate with the characteristics of the mathematical model are key parameters for classifying machine class. Ten
and five features are evaluated. Several training epochs are performed, and the ANN
classification results show high accuracy in all confusion matrices.
Machines in classes 1 and 2 are recommended for repair, whereas those in the tunable
machine class have the potential for PES improvement by adjusting the controller gain.
The mathematical model obtained from the system identification method will be used for
further controller optimization.
The classification method utilizing machine learning has demonstrated high accuracy,
greatly aiding users in identifying faulty parts instead of following a sequential fixing
process. Furthermore, we can leverage the model information from the tunable group to optimize the controller, rather than depending solely on a single controller gain. This approach
presents potential advantages in terms of efficiency and performance optimization.
Systems 2023, 11, 268
32 of 33
Author Contributions: Conceptualization, C.S. and J.S.; methodology, C.S. and J.S.; software, C.S.;
validation, C.S. and S.S.; formal analysis, C.S.; investigation, S.S. and J.S.; resources, C.S.; data
curation, C.S.; writing—original draft presentation, C.S. and S.S.; writing—review and editing, J.S;
visualization, C.S.; supervision, J.S.; project administration, J.S.; funding acquisition, J.S. All authors
have read and agreed to the published version of the manuscript.
Funding: This research was funded by the National Research Council of Thailand (NRCT), grant
N41A640442.
Data Availability Statement: Not applicable.
Acknowledgments: The authors are indebted to Suranaree University of Technology (SUT) and
Western Digital Storage Technologies (Thailand) Ltd. for their generosity and valuable comments.
Conflicts of Interest: The authors declare no conflict of interest.
References
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
Levy, E.C. Complex Curve Fitting. IRE Trans. Autom. Control 1959, 4, 37–43. [CrossRef]
Marti, J.R. Accurate Modelling of Frequency-Dependent Transmission Lines in Electromagnetic Transient Simulation. IEEE Trans.
Power Appar. 1982, PAS-101, 147–157. [CrossRef]
Jiraphon, S.; Autsadayut, R. Mathematical model of linear viscoelastic materials using weighting least square method. In
Proceedings of the 26th IASTED International Conference on Modelling, Identification, and Control, Innsbruck, Austria,
12–14 January 2007; pp. 452–457.
Mamun, A.A.; Lee, T.H.; Low, T.S. Frequency domain identification of transfer function model of a disk drive actuator. Mechatronics
2002, 12, 563–574. [CrossRef]
Taku, N. Identification of a Multiphase Network Equivalent for Electromagnetic Transient Calculations. IEEE Trans. Power Deliv.
2005, 20, 1134–1588. [CrossRef]
Eduardo, S.B.C.; Jose, A.G.R.; Bjorn, G. Rational Fitting Techniques for the Modeling of Electric Power Components and Systems Using
MATLAB Environment; Intech: London, UK, 2017; pp. 9–12. [CrossRef]
Takashi, Y.; Mitsou, H.; Chee, K.P. Chapter 2 System Modeling and Identification. In High-Speed Precision Motion Control, 1st ed.;
CRC Press Taylor & Francis Group, LLC: Boca Raton, FL, USA, 2012; pp. 24–32.
Chen, P.; Liang, Y. Model in Frequency-Domain Identification of a Fast-Steering Mirror System Based on Levenberg-Marquardt
Algorithm. In Proceedings of the 2017 2nd International Conference on Cybernetics, Robotics and Control (CRC), Chengdu,
China, 21–23 July 2017; pp. 199–202. [CrossRef]
Tomaz, K.; Damir, V. A simple analytical method for estimation of the five-parameter model: Second order with zero plus time
delay. Mathematics 2021, 9, 1707. [CrossRef]
Salem, M.A.A.; Zulkurnain, A.M.; Ali, A.S.; Zulkarnain, A.N.; Ahmed, A.A.; Mohd, F.M.Y.; Mohame, I.M.; Ahmed, A.S.;
Hammam, A.T. Frequency Response Analysis for Three-Phase Star and Delta Induction Motors: Pattern Recognition and Fault
Analysis Using Statistical Indicators. Machines 2023, 11, 106. [CrossRef]
Imad, A.; Vasilis, D.; Charilaos, M.; Konstantinos, T.; Eleni, C.; Nikolaos, D.; Keith, W.; Eoghan, A.M. Fault diagnosis of wind
turbine structures using decision tree learning algorithms with big data. In Proceedings of the 2018 European Safety and
Reliability Conference, Trondheim, Norway, 17–21 June 2018. [CrossRef]
Masayuti, S.; Prabhas, C. Abnormality Detection in Hard Disk Drive Assembly Process Using Support Vector Machine. In
Proceedings of the 2018 15th International Conference on Electrical Engineering/Electronics, Computer, Telecommunications and
Information Technology, Chiang Rai, Thailand, 18–21 July 2018. [CrossRef]
Aida, A.M.; Haidar, S.; Teymoor, G. k-NN based fault detection and classification methods for power transmission systems.
Prot. Control. Mod. Power Syst. 2017, 2, 32. [CrossRef]
Yordanos, D.M.; Yin-Der, L.; Jing, W.S.; Cheng-Cchien, K. Application of Machine Learning for Fault Classification and Location
in a Radial Distribution Grid. Appl. Sci. 2020, 10, 4965. [CrossRef]
Wanglomklang, T.; Chommaungpuck, P.; Chamniprasart, K.; Srisertpol, J. Using fault detection and classification techniques for
machine breakdown reduction of the HGA process caused by the slider loss defect. Manuf. Rev. 2022, 9, 21. [CrossRef]
Chommuangpuck, P.; Wanglomklang, T.; Srisertpol, J. Fault detection and diagnosis of linear bearing in auto core adhesion
mounting machines based on condition monitoring. Syst. Sci. Control. Eng. 2021, 9, 290–303. [CrossRef]
Prathan, C.; Siwanu, L.; Jiraphon, S. Fault Detection of Linear Bearing in Auto Core Adhesion Mounting Machine using Artificial
Neural Network. WSEAS Trans. Syst. Control 2019, 14, 31–42.
Phathan, C.; Thanasak, W.; Suradet, T.; Jiraphon, S. Fault Tolerant Control Based on an Observer on PI Servo Design for a
High-Speed Automation Machine. Machines 2020, 8, 22. [CrossRef]
Ahmed, R.N.; Ahmad, T.A.; Amjad, J.H.; Ammar, K.A.; Ibraheem, K.I. Intelligent fault detection and identification approach for
analog electronic circuits based on fuzzy logic classifier. Electronics 2021, 10, 2888. [CrossRef]
Zhu, Y.-P.; Zhao, Y.-L.; Lang, Z.-Q.; Liu, Z.-P.; Liu, Y. Online Rotor Systems Condition Monitoring Using Nonlinear Output
Frequency Response Functions Under Harmonic Excitations. IEEE Tran. Ind. Inform. 2022, 18, 6798–6808. [CrossRef]
Systems 2023, 11, 268
21.
22.
23.
24.
25.
26.
33 of 33
Naderi, E.E.; Khorasani, K. Data-driven fault detection, isolation and estimation of aircraft gas turbine engine actuator and sensors.
In Proceedings of the 2017 IEEE 30th Canadian Conference on Electrical and Computer Engineering, Windsor, ON, Canada,
30 April–3 May 2017. [CrossRef]
Samreen, J.; Dongyu, L.; Abhisek, U. Transient analysis method using high pass filter circuit in VSC interfaced multi-terminal DC
system. Electr. Power Syst. Res. 2023, 216, 109062. [CrossRef]
Et-Taleby, A.; Boussetta, M.; Benslimane, M. Faults Detection for Photovoltaic Field Based on K-Means, Elbow, and Average
Silhouette Techniques through the Segmentation of a Thermal Image. Hindawi Int. J. Photoenergy 2020, 2020, 6617597. [CrossRef]
Christoph, K.; Micha, K.; Michael, G.; Pascal, S. Classification of motor vibration with machine learning methods and simulating
the vibration using statistical models. In Proceedings of the 3rd Computer Science & Software Engineering Student Workshop,
Kryvyi Rih, Ukraine, 27 November 2020; Volume 2832.
Umbrajkaar, A.M.; Krishnamoorthy, A.; Dhumale, R.B. Vibration Analysis of Shaft Misalignment Using Machine Learning
Approach under Variable Load Conditions. Hindawi Shock. Vib. 2020, 2020, 1650270. [CrossRef]
GitHub. Available online: https://github.com/WillKoehrsen/feature-selector (accessed on 15 March 2023).
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.