Bearing Dynamics Modeling Based on the Virtual State-Space and Hammerstein–Wiener Model

This study investigates a novel approach for assessing the health status of rotating machinery transmission systems by analyzing the dynamic degradation of bearings. The proposed method generates multi-dimensional data by creating virtual states and constructs a multi-dimensional model using virtual state-space in conjunction with mechanism model analysis. Innovatively, the Hammerstein–Wiener (HW) modeling technique from control theory is applied to identify these dynamic multi-dimensional models. The modeling experiments are performed, focusing on the model’s input and output types, the selection of nonlinear module estimators, the configuration of linear module transfer functions, and condition transfer. Dynamic degradation response signals are generated, and the method is validated using four widely recognized databases consisting of accurate measurement signals collected by vibration sensors. Experimental results demonstrated that the model achieved a modeling accuracy of 99% for multiple bearings under various conditions. The effectiveness of this dynamic modeling method is further confirmed through comparative experimental data and signal images. This approach offers a novel reference for evaluating the health status of transmission systems.


Introduction
The mechanical transmission system plays an indispensable key role in transportation.For example, over 45% of maritime machinery failures are caused by transmission systems.Therefore, developing intelligent operation and maintenance technologies for transmission systems is imperative.Prognostics and Health Management (PHM) has emerged as a valuable strategy for preventing failures and reducing costs, garnering increasing research interest.Bearing failure, which accounts for over 50% of mechanical failures, can be a significant research object.This research aims to provide new methods of fault modeling and ideas for developing health management technology for transmission systems.
In PHM, fault modeling techniques are categorized into mechanism-based and datadriven methods.Mechanism-based methods employ mathematical and physical models, as illustrated in method (a) in Figure 1 [1].Bearing degradation is typically characterized using vibration and temperature sensors, with vibration signals from acceleration sensors being particularly prevalent due to their sensitivity and rapid response.A significant amount of research in bearing PHM has focused on mechanism modeling [2][3][4].Notable advances include dynamic models based on Hertz contact theory and finite element methods for detecting defects in angular contact ball bearings [5], as well as more complex models incorporating local geometrical imperfections and varying race-supporting conditions [6].Qian et al. [7] developed a model-based fault diagnosis technology for rotor-bearing systems, effectively utilizing residual generation techniques for fault detection and localization.
Despite these advancements, the application of mechanism-based methods in industrial settings is limited due to noisy and complex operating environments.Consequently, data-driven methods are increasingly prominent in bearing PHM research.These methods leverage machine learning techniques to capture dynamic degradation characteristics from measurement data and produce high-precision models.These techniques include shallow and deep learning techniques, as depicted in methods (b) and (c) in Figure 1.Since Janssens et al. first applied Convolutional Neural Networks (CNNs) to bearing fault diagnosis in 2016 [8], various enhancements have been proposed, such as 2D-CNN, multiscale CNN, and adaptive CNN [9][10][11].Other techniques like Long Short-Term Memory (LSTM), Deep Belief Networks (DBN), Recurrent Neural Networks (RNN), and Deep Deterministic Policy Gradient (DDPG) algorithms have also demonstrated efficacy in fault diagnosis and modeling [12][13][14][15][16].For instance, Tian et al. proposed a CNN-LSTM fault diagnosis model optimized using Hybrid Particle Swarm Optimization (HPSO), achieving a classification accuracy of 99.2% [17].However, data-driven modeling methods face challenges, such as the need for extensive training datasets and a lack of interpretability, particularly in hyperparameter optimization.To address these challenges, this study proposes a bearing dynamics modeling method utilizing virtual state-space and Hammerstein-Wiener (HW) models.This approach generates the necessary data for modeling by constructing multi-dimensional virtual states.It enhances model interpretability through the HW model's capacity to capture nonlinear system characteristics, thereby improving modeling accuracy.Although HW models are predominantly studied in the control field [18,19], their robust data fitting capabilities for nonlinear systems suggest the potential for high-precision dynamic modeling of bearing degradation and anomalies.Additionally, the HW model's modular structure offers enhanced configurability and extensibility.For example, Wang et al. integrated HW models with neural networks, achieving lower average absolute percentage errors (0.0223) compared to Gated Recurrent Units (GRU) and Recurrent Neural Networks (RNN) on simulation data [18].Zhong et al. [19] explored parameter estimation for HW time-varying systems using a learning recognition algorithm with a forgetting factor.Based on these advancements, this study combines the HW model with virtual state-space for fault modeling in PHM, validating it through dynamic modeling of bearings and providing a novel approach for assessing transmission system health.The innovations of this study include the following: 1.
This study uses multi-dimensional virtual state-spaces to generate sufficient data, improving the accuracy of model identification.

2.
This study constructs multi-dimensional models of the dynamic system by utilizing the virtual state-spaces, which allows for high-precision modeling, mainly through the dimension with a more vital characterization ability of the dynamic system.

3.
This study employs the HW model in the control field to effectively mitigate the challenges posed by the nonlinear part of the system in dynamic modeling.
Figure 2 illustrates the research structure, comprising four components: (a) methods for obtaining high-quality signals and selecting valid data; (b) extraction of virtual states for fault modeling through dynamic analysis; (c) application of virtual state-space to bearing fault modeling and HW method for model identification; and (d) verification experiments of the proposed method.The paper is organized as follows: Section 2 covers theoretical foundations, Section 3 describes the database and data processing, Section 4 presents experimental results for modeling configuration, Section 5 validates the optimal configuration and method with various datasets, and Section 6 concludes the study.

Bearing 5-DoF Model for Dynamics Analysis
This study uses bearings as a case study to investigate a new fault modeling method.The foundational analysis employs the 5-DoF model, central to subsequent research.This model, shown in Figure 3, is mathematically represented by Equations ( 1)-( 5) [20]: where m ir , m or , m r , respectively, represent the mass of the inner race, outer race, and the unit resonator; k ir , k or , k r , c ir , c or , c r represent their stiffness and damping, respectively; f x , f y represent the nonlinear contact force of the bearing in the horizontal and vertical direction, respectively; x ir , y ir , x or , y or , y r mean the two inner and outer race DoF, and measured vibration response, respectively.Additionally, g represents the gravitational acceleration, and F load denotes the vertical external force applied to the inner race of the bearing [21].The system is simplified using the lumped parameter method by treating all balls as a single entity.The contact force f j between the j-th ball and the raceway is split into horizontal and vertical components ( f x j and f y j ).Then, through Newton's second law, the vertical acceleration a y and can be given by Equation ( 6) [3]: where n b denotes the number of balls in the bearing and m j the mass of the j-th ball.From the Hertz contact theory, the contact force f j and its vertical component f y j between the j-th ball and the raceway are given by Equations ( 7) and (8), respectively: where k b represents the stiffness of the balls.γ j is a switch function indicating the ball is within the load region, thereby enabling deformation due to the contact force f j .The total deformation δ j of the j-th ball is related to the relative displacement between inner and outer races, its angular position ϕ j , and the bearing clearance.γ j and δ j are calculated by Equations ( 9) and (10), respectively [3]: where c stands for the clearance.Many essential state quantities in the dynamic model are difficult to measure in practice, such as (y ir and y or in Figure 3, which makes it difficult to establish an accurate physical model.Therefore, this study will explore constructing a high-precision fault model based on dynamic analysis from the perspective of acceleration response and multi-dimensional data.

Virtual State-Space for Dynamics Modeling
This study's theoretical foundation is based on the concept of virtual states and the virtual state-space model for generating multi-order data and exploring system dynamic characteristics from multiple dimensions [22].The structural and compositional expressions of the virtual state-space align with the general state-space formulation, as illustrated in Equations ( 11) and (12): where the state matrix A, input matrix B, output matrix C, and feed-forward matrix D retain their standard definitions and functions.The primary distinction lies in the virtual states within the virtual state vector x v .These virtual states possess undefined but real physical significance, which can be represented by one or more derivatives of the actual state.Consequently, the variables in the virtual output vector y v and virtual input vector u v can also be obtained by derivation of actual measurable variables.
Based on the relationships between various actual state variables in the dynamic model, y ir , y or , and y r can be selected as the actual state variables in the initial bearing model (0thorder).The derivatives of these states retain specific physical significance, fulfilling the requirements of the theory of virtual states.Additionally, the vertical acceleration signal ÿr is a measured physical variable associated with the state variables y ir , y or and y r .Therefore, ÿr can be used as an in-or output of the initial model and expressed as a y = ÿr .As an example of a three-input single-output model, Equations ( 13) and ( 14) can be used to construct the initial and virtual state-space models for the bearing: Higher-order virtual state-spaces can be derived from the initial state-space using different input forms, as depicted in Figure 4. Multi-order virtual state-space models can also be developed through high-order derivatives of initial inputs and outputs, providing a data foundation for model identification.In virtual state-space equations consisting of various orders, variations in the state matrix A n reflect the degradation of bearing dynamics in multiple dimensions, with specific parameters that may contain information regarding bearing aging and faults.Thus, virtual state-spaces serve two primary purposes: generating extensive multi-dimensional data from limited data and exploring system dynamic degradation from various dimensions to uncover implicit dynamic characteristics.
Based on the inference of generating multi-dimensional data, multiple-order derivatives of the acceleration signal a y can serve as in-and outputs for the virtual state-space model.Hence, the in-and output vectors can include one or more of the acceleration signal a y and its multiple-order derivatives a (i) y , where i = 1, 2, 3, . . ., n, (n ∈ N + ).Although the time function of a y is assumed to be unknown during model identification, high-frequency sensors provide the minimum interval between sampling points.This approach allows for the use of finite difference rather than differentiation as an alternative method to generating multi-dimensional data.Low-order derivatives (e.g., velocity and acceleration) of physical variables (e.g., displacement) are typically well-defined and measurable, making them suitable for input or output in virtual state-space models.However, higher-order derivatives are more challenging to obtain.Therefore, finite difference can be employed in place of differentiation to calculate derivatives in specific scenarios, such as sampling at very high frequency.The minimal time intervals (∆t ≤ 0.0001 s) between points t 0 and t 0 + ∆t make it justifiable to replace differentiation with finite differences, as represented in Equation ( 15): Simultaneously, filtering the generated data can help mitigate the adverse effects of using finite differences to calculate derivatives.The entire workflow for utilizing virtual state-space in bearing fault modeling is depicted in Figure 5.

Hammerstein-Wiener Method for Model Identification
Based on the framework of multi-dimensional virtual state-space, this study proceeds to identify dynamic systems with known inputs and outputs.Given the complexities of bearings in real-world applications and the inherent nonlinearity of their acceleration signal fluctuations, the HW modeling method is employed.While HW models have demonstrated significant success in various nonlinear systems and control applications [23], their application in PHM is relatively underexplored.Thus, a key innovation of this study is integrating the HW model with virtual state-space, which facilitates high-precision model identification for multi-dimensional bearing dynamics and explores new methods for assessing the health of transmission systems.
The HW model can serve as a black-box for estimating linear models, with its accuracy enhanced by incorporating nonlinearities in the inputs and outputs.It can also function as a grey-box model that captures physical process characteristics.As illustrated in Figure 6, the HW model simplifies a complex nonlinear dynamic system into two static nonlinear blocks, f and h, in series with a linear block L, where u(t) and y(t) represent the input and output data, respectively.The the linear block's input and output are denoted as w(t) and x(t).The model structure is mathematically described by Equations ( 16)-( 18): where α, β, and ) are parameter vectors of memoryless nonlinearities f and h, and the linear time-invariant system L. The HW model represents the dynamics of the polynomial-based linear block L within the system using the discrete transfer function H(s) in Equation ( 19) [24]: where N Z and N p represent the number of zero poles, respectively, z n , z * n , p k and p * k denote the values of zeros and poles along with their complex conjugates, and B represents the scaling coefficient.For instance, in a two-input-two-output system, the specific Hammerstein nonlinearity f and Wiener nonlinearity h are given by Equations ( 20) and ( 21) [25]: where f 1 and h 1 represent the input and output nonlinearities corresponding to the first input u(t, 1), while f 2 and h 2 correspond to those associated with the second input u(t, 2).This block structure provides strong data fitting capabilities, reconfigurability, and parameter interpretability.The HW model does not require complete identification of the nonlinear process; instead, it uses nonlinear functions f and h to process the linear system's input and output signals, effectively capturing the system's nonlinear characteristics.In the modeling experiment, attention is given to the stability, accuracy, and efficiency of configuring the nonlinear functions f and h.Firstly, the nonlinear functions must avoid introducing mathematical errors to ensure process stability.Secondly, high accuracy is crucial to reflect the model's fitting capability.Lastly, the modeling process should balance stability and accuracy to achieve efficiency.The dynamic components of the linear module involve three key parameters: input delay, and number of zeros and poles.Input delay influences the system's response to external stimuli over time, while the quantity and values of zeros and poles affect the stability and performance of the dynamic system.

Database and Data Preprocessing
The proposed method is evaluated using the IEEE PHM 2012 Challenge dataset, a widely recognized resource in the field of PHM.This dataset's rolling bearing test data are collected from the PRONOSTIA experimental platform at the FEMTO-ST research institute and have been previously used in the IEEE PHM 2012 Prognostics Challenge.The PRONOSTIA platform and accelerometers are shown in Figure 7 [26].
In the measurement part of the test bench, various high-quality sensors are used to obtain the instantaneous radial force applied to the bearing, the speed of the shaft handling the bearing, and the torque applied to the bearing.The acceleration sensors (Type DYTRAN 3035B) operate in vertical and horizontal directions at a sampling frequency of 25.6 kHz.Additionally, this sensor collects samples for 0.1 s every ten seconds.The test bearing is illustrated in Figure 8, and its parameters are listed in Table 1.Bearing accelerated degradation experiments are conducted until the vibration signal amplitude exceeds 20 g.As detailed in Table 2, these run-to-failure datasets include 17 rolling bearings tested under three different operating conditions.The operating conditions in the experiments are set to higher levels than typical industrial scenarios to expedite bearing degradation.To validate the generalizability of the proposed model, additional comparative experiments are performed using datasets from Case Western Reserve University (CWRU) [27], the Society for Machinery Failure Prevention Technology (MFPT) [28], and Paderborn University (PU) [29].Their parameters are summarized in Table 3. Figure 9 illustrates the data flow from sensor acquisition to modeling.During degradation experiments, data from the bearing's healthy state throughout most of its lifecycle are not essential for the study.Modeling these stable acceleration responses is less challenging but consumes considerable computational resources.Therefore, we focus on the aging stage during data preprocessing.Acceleration signals are filtered using a one-dimensional stationary wavelet transform with Wavelet Symlets 8 and a frequency resolution of 50 Hz.The boundaries between the healthy, aging, and failure stages of the bearing are determined using the 3σ method.The division between the healthy and aging stages is completed by calculating the kurtosis value, while the starting point of the failure stage is determined by referring to the Root Mean Square (RMS).After segmenting the data from the bearing's aging stage, finite difference is performed to generate multi-dimensional signals, filtered using a median filter with a window size of five to remove impulse noise from the finite difference process.These filtered finite difference data are used as inputs for the virtual state-space and further applied for model identification.

Experimental Results
This section presents sets of experiments to explore the effectiveness of the HW model based on virtual state-space and the impact of its different configurations.Critical configurations for HW modeling include the input and output structures of the state-space, their corresponding nonlinear estimators, and the zeros and poles of the transfer function.HW modeling is implemented using MATLAB's System Identification Toolbox, with the Normalized Root Mean Square Error (NRMSE) as the criterion for evaluating modeling accuracy.NRMSE measures the fit between the model's response and the actual data, expressed as a percentage (Modeling Accuracy = 100(1 − NRMSE)).The time required for modeling under various configurations is recorded to assess modeling efficiency.All experiments are conducted on a consistent software and hardware platform, and time measurements are averaged over three trials.The transfer function's delay in the linear subsystem is set to one.Unless specified otherwise, all experiments use 126 samples collected from Bearing 3 3 during the aging phase.

Model In-and Output
Initially, the composition of the model's inputs and outputs needs to be determined.Single-input and single-output models are limited in their configurability and cannot simultaneously consider operating conditions and measurements, so they are not included in this study.Additionally, due to potential discrepancies between differentiation and finite difference data affecting HW model stability, particularly in complex Multi-Input-Multi-Output (MIMO) scenarios, the study focuses exclusively on Multi-Input-Single-Output (MISO) models.
Two primary configurations for the model inputs are examined.The first type includes a y and its multiple-order derivatives: u = [a y , ȧy, . . ., ay (n) ], where n ∈ N + .The second type incorporates operating conditions to the first type, such as rotational speed n r and torque T N , in addition to the derivatives: u = [n r , T N , a y , ȧy, . . ., ay (n) ].The study constructs a Multi-Input-Single-Output-Hammerstein-Wiener (MISO-HW) model based on the statespace representations in Equations ( 13) and ( 14) and evaluates the impact of four basic input-output compositions on modeling accuracy, as shown in Table 4.The linear module's transfer function is configured with three poles and zero zeros.The experiments use the same nonlinear estimator for input and output modules, and various nonlinear estimators are compared for modeling accuracy, as shown in Table 5. Results indicate that the modeling accuracy is significantly low when the input data's order is lower than the output data's.Therefore, Type 1 and Type 2 are not suitable for HW modeling.In Type 4, when using Unit Gain as the nonlinear estimator, the model accuracy reaches 98.70%.However, due to mathematical constraints, this configuration results in computational errors when validating modeling on other bearings.Similar errors occur in Type 2 as well because using finite differences instead of differentiation during data generation results in errors in the calculation process of a particular sample, leading to poor model stability.Therefore, we choose Type 3 as the input configuration, where the model input can be multiple variables of higher order than the output while disregarding operating conditions.
With the input data composition established, the study explores the HW model's capacity to capture dynamic characteristics across different dimensions by varying input orders.It is assessed by modeling accuracy.With the experimental configuration and sample data constant, Unit Gain (UG) and One-Dimensional Polynomial (1-DP) are utilized as nonlinear input and output modules, and first through tenth-order derivatives of the acceleration signal serve as inputs.The original acceleration signal is used as the output.The resulting experimental data are presented in Table 6.The highest modeling accuracy is achieved when using three inputs with second-order data.However, increasing the order of inputs beyond the second order does not improve the model fitting results.Instead, it leads to decreased modeling efficiency due to increased computation time.

Nonlinear Module Estimator
Since the input-output configuration is determined, this set of experiments focus on selecting appropriate in-and output estimators to match the bearing's nonlinear characteristics according to model fitting accuracy.The experiments use the optimal input-output configuration from Table 6, varying only the estimators.Due to numerous combinations, the same nonlinear estimator is initially applied to both input and output modules, as shown in the first group of experiments in Table 7. Results reveal that 1-DP and UG estimators achieve higher modeling accuracy, exceeding 87%.These findings suggest that these estimators better align with the nonlinear characteristics of the MISO-HW model's input and output modules.
Consequently, UG is employed for input and output modules, as demonstrated in the second and third groups of experiments shown in Table 7.With the UG-1-DP combination, the modeling accuracy achieves a peak of 97.8659%, while the modeling time is 181.21 s.This time is slightly longer than the 134.87 s required for the 1-DP-UG combination, which, despite a lower accuracy improvement of 18.8475%, demonstrates superior efficiency.Further tests are conducted to assess the performance of the 1-DP combined with other estimators.Although the SN-1-DP combination yielded an accuracy of 96.0861%, its modeling time is significantly longer at 445.57 s, 2.46 times that of the UG-1-DP combination.This comparison indicates that the MISO-HW model performs optimally, with UG and 1-DP as input and output nonlinear estimators.

Zeros and Poles of Transfer Function
After identifying the optimal configuration for the nonlinear modules, further experimentation focuses on optimizing the model by adjusting the linear module settings.Various experiments are conducted to explore the impact of altering the number of zeros and poles in the transfer function of the HW model's linear module.As per the requirement that a physically feasible control system must have more poles than zeros, the experimental design is outlined in Table 8.All configurations, except for the number of zeros and poles, use previously determined optimal values.In the first test group, increasing the number of poles improves modeling accuracy at first, reaching its peak with three poles before fluctuating downwards.Models with more than one zero exhibit instability, and all models with more than three zeros could not be constructed stably.Thus, it is inferred that the HW model is more stable when configured with no zeros.Both linear modules with three or eight poles achieve fitting accuracy exceeding 97.5% with a negligible computation time difference of 4.68 s.Therefore, the configuration with no zeros and three poles can be used for the highest modeling accuracy.

Analysis
Analysis of the four sets of experimental results reveals significant differences in how various parameters affect model accuracy.Firstly, the setting of zeros and poles in the linear module directly influences modeling stability.Based on the results in Table 8, configuring the transfer function of the linear module with no zeros allows the MISO-HW model to more accurately capture the linear characteristics of bearing dynamics.Secondly, the input-output configuration is crucial for achieving high model precision.Thirdly, among nonlinear modules, the choice of estimator impacts modeling accuracy, though to a lesser extent than the input-output configuration.Fourthly, the input data dimension has a relatively minor effect on accuracy, ranging from 2.7%.The experiments in Table 6 show that accuracy did not improve with an increasing number of inputs.However, suppose the modeling accuracy increases with the number of inputs.In that case, it suggests that higher dimensions of data and model lead to high accuracy and better performance, which is different from the expected conclusion.This set of experiments focuses on the model dimension, aiming to identify model dimensions that enhance modeling performance.The additional input between different experiments is higher-dimensional data generated based on the virtual state theory.Thus, increasing the number of inputs expands the model dimension and reconstructs the mapping relationship instead of simply adding data.For example, the modeling accuracy declines when threeorder derivative data are utilized, as shown in Table 6.This decline is attributed to the relatively simple system structure, where three dimensions suffice to effectively describe the dynamic characteristics of the bearing's steady-state operation.Capturing the system's dynamic feature in higher dimensions is more challenging, resulting in less contribution of higher dimensions to the modeling performance, thereby reducing the accuracy of the entire model.
The correctness of experimental results in Table 7 can be supported by the analysis of the bearing 5-DoF model.First, the 1-DP in the nonlinearity block of the HW model can be expressed by Equation ( 22): where c(n p ) denotes the polynomial coefficients for the n p -th term.Second, according to the Hertz contact theory, Equations ( 6)-( 8) reveal that the acceleration a y is positively correlated with the 1-DP of deformation (δ 1.5 j ) and UG composed of other parameters (k b , γ j , sin ϕ j ).By observing the similarity of these functional forms, we can conclude that it is mathematically interpretable to generate a high-precision bearing model using the 1-DP as a nonlinear module.

Validation of Optimized Configuration
Summarizing the comprehensive experimental results, the optimal configuration of the MISO-HW bearing model includes the transfer function with no zeros and three poles in the linear module, estimators employing UG and 1-DP in input and output nonlinear ] and y = a (n) y as input-output type.Under this configuration, samples from Bearing 3_3 are used to construct 0th to ninth-order virtual state-space models.Herein, the original acceleration signal and its 0th to ninth-order derivatives are used as the outputs for the ten virtual state-space models.
After excluding outliers using the three-standard deviation method, the precision of the MISO-HW model constructed for Bearing 3 3 consistently exceeds 90%.As detailed in Table 9 corresponding to Figure 10, the average precision across all state orders surpasses 93.72%, with relatively small standard deviations.Notably, for virtual state-space models with orders less than three, the modeling precision exceeds 98.29%, with a standard deviation within 3.64.This suggests that lower-order virtual state-space models may better capture the system's dynamics.In conclusion, this modeling approach proves to be effective for Bearing 3 3 .As a supplementary analysis, Figure 11 compares the original signal (in red) from the virtual state-space model for a single sample of Bearing 3 3 with the HW model output signal (in green).The right side of the figure provides an enlarged comparison of six sampling points, where the high modeling accuracy results in a significant overlap between the two curves.Therefore, additional comparison plots with greater precision are presented in the Appendix A, as shown in Figure A1.The nearly identical trends and values of both signals further validate the efficacy of this modeling approach.

Verification with Bearings under Other Conditions
Following these initial findings, it is crucial to validate the proposed method across different bearings and operating conditions.With the same modeling configuration, experiments are conducted on 1700 samples from 17 bearings under all working conditions listed in Table 2. Figure 12 presents the experimental results for two bearings under each condition, with additional numerical results available in Table A1 in the Appendix A. The results indicate that all modeling accuracies exceed 90%, with an average accuracy of 96.30%.For instance, Figure 13 illustrates the modeling results for two bearings under varying operating conditions.Box plots, obtained after removing outliers, show that the average modeling accuracy for all bearing samples exceeds 95% under the NRMSE standard, with most samples achieving over 92% accuracy.Bearings under Condition 2 demonstrate an average accuracy exceeding 99% for virtual state-space models up to the fifth order.Thus, the MISO-HW model based on virtual state-space is highly effective for bearing fault modeling and exhibits robust generalization capabilities.

Verification with Other Datasets
Comparative experiments using the CWRU, MFPT, and PU datasets are presented in Figure 14.Two types of bearing fault samples from each dataset are modeled with state orders ranging from 0 to 9. The number of sampling points for each bearing sample ranges from 100,000 to 120,000, with relevant parameters detailed in Table 3.The modeling accuracy for the six bearings across the three datasets consistently exceeds 90% in all state orders, with average accuracies of 94.44%, 95.34%, 97.28%, 97.76%, 96.96%, and 98.18%, respectively.The location of the bearing fault does not significantly impact the modeling accuracy.The modeling accuracy of the two CWRU samples (12.0 kHz) is lower than that of other bearing datasets (48.8 kHz and 64.0 kHz), presumably due to differences in data quality.The sampling frequency of these two bearings is lower than that of the other samples, and given that the data acquisition predates the other datasets by more than a decade, there may be gaps in equipment performance.Nonetheless, the high modeling accuracy of bearing samples from different datasets verifies the effectiveness of the proposed method and demonstrates its universality and generalizability.

Discussion
The first three orders depicted in Figure 12 reveal that the average modeling accuracy for the three bearings under Condition 3 is 96.84%, 98.57%, and 97.32%, respectively.For the 17 bearings overall, the average modeling accuracy exceeds 97.26%.These results further substantiate the effectiveness of this method in object transfer.Additionally, when assessing modeling accuracy across various working conditions, the average accuracy for the first three orders are 96.66%,98.25%, and 97.57%, respectively.The MISO-HW modeling method, optimized with bearing samples from Condition 1, demonstrates robust performance across bearings in Conditions 2 and 3, underscoring the efficacy of this approach in transferring across working conditions.
Considering model dimensionality in Figure 14, the highest accuracy for each bearing is predominantly achieved at state orders 4 and 6.However, in Figure 12, the peak accuracy is concentrated at state orders 0 to 2. The data in Figure 12 are derived from healthy bearings in the aging stage, while Figure 14 reflects data from faulty bearings.This distinction suggests that the bearings in Figure 14 exhibit more significant dynamic nonlinearity and signal complexity.Thus, a higher-dimensional model is necessary to accurately capture these bearing faults' dynamic characteristics, reinforcing the research value and motivation behind the proposed method.

Conclusions
This study presents a bearing fault modeling method based on virtual states and the HW model.The approach involves analyzing the 5-DoF dynamic model of bearings to identify suitable virtual states and input-output configurations, which are then used to establish virtual state-space models.Bearing vibration signals collected by acceleration sensors under various working conditions are employed in the experiments.Extensive testing of HW modeling configurations result in highly accurate bearing dynamics models, validating the effectiveness of the proposed method.The key contributions of this work can be summarized as follows: 1.
The concept of virtual states is introduced, highlighting its significance in exploring system dynamic characteristics from multiple dimensions.This study employs virtual states and multi-dimensional bearing state-space models derived from measurable states in the dynamic analysis of bearings.

2.
A bearing fault modeling method is developed using virtual state-space and HW models, utilizing acceleration signals and their multi-order derivatives to construct high-precision MISO-HW models based on virtual states.
The virtual state theory offers a theoretical framework for describing the dynamic characteristics of systems across various scales and dimensions, providing new insights for system modeling, identification, and fault diagnosis.While the proposed method shows great potential in PHM, it also presents certain limitations and areas for improvement.For instance, this method is currently restricted to MISO models due to the use of finite difference rather than differentiation in data generation, which compromises data quality and fails to meet the stability requirements of MATLAB's built-in algorithms during the modeling process.Based on these findings, the following research directions are proposed: 1.
To optimize data generation and improve the quality of high-dimensional data, derivative estimators could be developed, enabling the construction of MIMO models.These MIMO models allow for the extraction of characteristic parameters that more effectively represent the dynamic degradation processes in multidimensional models.

2.
The state matrix A n in the models (refer to Figure 4) captures the relationships among states and the system's dynamic characteristics across different dimensions.This matrix allows for the extraction of characteristic parameters highly correlated with bearing aging and faults, facilitating the development of high-precision diagnostic and prognostic methods.

3.
This method can be further applied to different data sets and objects.A test bench equipped with a bearing test module and high-quality data acquisition terminals has been constructed.Experiments are planned to create a comprehensive dataset to validate the proposed method further.Furthermore, gear and rotor data will be employed to test this method following the design of a new test bench and relevant experiments.

4.
Future improvements in the accuracy of this modeling method could be achieved by incorporating varying friction states.For example, integrating a friction coefficient that adapts to different friction conditions could enhance the model's precision.Additionally, adopting an improved elastohydrodynamic lubrication modeling approach could serve as a superior alternative to the current Hertz contact theory-based modeling.

Figure 5 .
Figure 5. Construction process of virtual state and virtual state-space: (a) 5-DoF model, (b) dynamics analysis, (c) virtual state generation, (d) virtual state-space construction.

Figure 7 .
Figure 7. Test bench and sensors: (a) the PRONOSTIA platform, (b) measuring points in vertical and horizontal directions, (c) accelerometer.

Table 1 .
Parameters of bearings.Parameter Value Diameter of the outer race D or 29.1 mm Diameter of inner race D ir 22.1 mm Bearing mean diameter D m 25.6 mm Inside diameter D i 20.0 mm Outside race diameter D o 32.0 mm Diameter of rolling elements D b 3.5 mm Thickness 7.0 mm Number of rolling elements n b 13

Figure 9 .
Figure 9. Data processing (aging data refers to the bearing data in the aging stage).

Figure 12 .
Figure 12.Modeling accuracy of bearings under different conditions.

Figure 14 .
Figure 14.Modeling accuracy of bearings with faults in Outer Race (OR) and Inner Race (IR) in different datasets.

Figure A1 .
Figure A1.Signal comparison of HW model for the last sample of Bearing 3_3.

Table 3 .
Spectrum error in ablation experiment results.

Table 4 .
Configurations of the input and output.

Table 5 .
Modeling accuracy with different configurations of the input and output.

Table 6 .
Configuration experiments with different input orders.

Table 7 .
Modeling accuracy with different nonlinearity estimators.

Table 8 .
Modeling accuracy with different configurations of zeros and poles.

Table A1 .
Accuracy of bearings under different working conditions [%].