Suspension parameter identification method for rail transit vehicles using an AO-GRBF surrogate model and non-dominated sorting genetic algorithm

1.1. Novelty and contribution

Within the proposed AO-GRBF–NSGA-II framework, the principal contributions of this study are summarized as follows:

(1) To improve the stability of the GRBF surrogate model for suspension parameter identification under real operating conditions, an error-driven sampling mechanism is introduced to improve local fidelity in high-error regions. In addition, a multiobjective optimization scheme is employed to automatically adjust kernel-related parameters that determine model performance.

(2) To prevent convergence to local optima during parameter identification, a modified NSGA-II algorithm is adopted. This algorithm enables global search through population evolution combined with Pareto-based selection strategies.

(3) To demonstrate practical feasibility and deployment potential, the proposed method is validated using field-measured data obtained from rail transit operations.

The core technical challenge is to maintain accuracy and reliability when computationally intensive mechanistic models of suspension systems are replaced by data-driven surrogate models. This study addresses this challenge through residual-driven sampling, kernel self-tuning, and advanced evolutionary optimization.

1.2. Significance

Conventional approaches are constrained by lengthy offline computations, which limit the ability to detect suspension parameter variations in real time. As a result, early anomalies tend to be overlooked and the effectiveness of predictive maintenance is reduced. By contrast, the proposed method is designed to enable rapid parameter identification and timely detection of deviations, so that preventive actions can be taken before safety risks escalate. This advantage is expected to strengthen both railway safety and operational efficiency.

Beyond methodological innovation, the combined AO-GRBF and NSGA-II framework offers practical engineering benefits. First, an error-driven sampling strategy is employed to improve the robustness of the surrogate model to noise and non-stationary excitations, thereby ensuring stable performance under real operating conditions. Second, compared with repeated high-fidelity SIMPACK simulations, the surrogate-based approach significantly reduces computational cost. Consequently, efficient condition monitoring and predictive maintenance become more feasible in practical railway applications.

2.1. Mechanistic approaches

Mechanistic methods are based on nonlinear finite element (FE) or multibody dynamics (MBD) models of railway vehicles. These models are typically implemented on industrial platforms such as SIMPACK or Ansys to simulate system responses under parameter variations ^[7,8]. Suspension parameters are identified by solving an inverse problem in which discrepancies between simulated and measured responses are minimized ^[8]. For instance, Zhang et al. developed a mechanistic model for high-speed trains and determined stiffness and damping coefficients using a min-max method ^[6]. Yao et al. employed a simplified lateral dynamic model combined with NSGA for robust stability identification ^[20]. Kraft et al. formulated parameter identification as an adjoint-state inverse problem, which enabled accurate calibration ^[8]. Although these methods provide physically interpretable results, their high computational costs limit practical field applications ^[7,8,21,22].

2.2. Surrogate model-based approaches

To reduce computational costs, surrogate models such as Kriging, radial basis function (RBF), and neural networks have been employed to approximate the mapping from parameter inputs to suspension system performance responses; subsequently, these surrogates are coupled with optimization algorithms to identify parameter sets that best correspond to actual operating conditions ^[9–16].

Among these surrogates, RBF models are widely used for nonlinear black-box systems because they provide mesh-free interpolation with competitive accuracy and training cost in high-dimensional settings ^[23]. However, performance depends critically on the sampling policy: when constructed from static, one-shot designs, RBF surrogates tend to miss localized nonlinear dynamics, and accuracy degrades across the high-dimensional suspension parameter space ^[17–19]. To address these issues, adaptive RBF frameworks have been proposed. For example, Liu et al. introduced a trust-region-based adaptive RBF (TARBF) that sequentially enriches samples within moving trust regions to improve sample efficiency and local fidelity for computationally expensive constrained problems^[18]. Ye et al. presented an RBF-assisted adaptive differential evolution (RBF-ADE) with a cooperative dual-phase infill strategy to accelerate convergence and enhance surrogate accuracy in high-dimensional expensive tasks^[19]. Nevertheless, in both approaches, the basis-function class remains fixed, and kernel hyperparameters are not jointly optimized with the sampling policy. Moreover, direct integration of adaptive RBF surrogates into suspension parameter identification remains scarce, which limits their application under in-service conditions.

In addition, as suspension parameter identification is fundamentally a surrogate-based optimal-parameter search, the multiobjective evolutionary algorithm NSGA-II is commonly adopted to prevent premature convergence and maintain population diversity. For example, Zhang et al. developed a neural-network surrogate to replace computational fluid dynamics and vehicle-dynamics models under crosswind conditions and coupled it with NSGA-III to identify bogie suspension parameters ^[14]. Zhang et al. also established a deep-learning surrogate to approximate heavy-haul freight train dynamics and employed NSGA-II to estimate key longitudinal stiffness and damping coefficients ^[15]. Jiang et al. constructed surrogate models of parameterized suspension components and integrated them with NSGA-II for lightweight optimization ^[16].

2.3. Novelty of the proposed method

This study employs two key algorithms. First, a data-driven AO-GRBF surrogate is introduced to replace the computationally intensive mechanistic model of the suspension system while retaining fidelity under variable in-service conditions. Within AO-GRBF, residual-driven sampling and in-loop kernel updates operate in a closed loop to target difficult regions of the high-dimensional, nonlinear suspension parameter space. After each enrichment step, the kernel shape and regularization are updated to fit local curvature, improving local fidelity and preserving robustness as operating conditions drift.

Second, NSGA-II is employed to define a surrogate-based multiobjective optimization framework, converting parameter identification into a tractable model-optimization problem. Pareto-based global evolution maintains population diversity and conducts a reliable search across the suspension parameter space to balance accuracy and robustness, yielding estimates that best match actual operating conditions.

3.1. Construction of the mechanism- and data-driven AO-GRBF surrogate model

To expedite suspension parameter identification, a data-driven AO-GRBF surrogate model for the suspension system is constructed in this section. The approach comprises: (1) mechanism-based training dataset generation via MBD modeling; (2) construction of the data-driven AO-GRBF model with adaptive optimization techniques; and (3) model performance validation.

3.1.1. MBD model of rail transit vehicle

Scenarios involving suspension parameter degradation cannot be safely reproduced through real-vehicle testing and pose substantial safety risks. To enable simulation and obtain vehicle vibration responses under different suspension parameter configurations, a train multibody dynamics model is constructed on the SIMPACK platform. The model incorporates critical nonlinear characteristics, including axle-box suspension clearances and dry friction in friction wedges ^[24,25]. Additionally, the wheel-rail contact is represented using Hertzian theory combined with the FASTSIM algorithm ^[26,27].

The validated train multibody dynamics model ^[28,29] comprises both primary and secondary suspension systems (highlighted in red boxes in Figure 2), which are designed to control load bearing and ensure ride stability. Fourteen suspension parameters ($$ x_1^{o} $$ to $$ x_{14}^{o} $$) are included in these systems, as listed in Table 1. To avoid computational redundancy during subsequent parameter identification, the effect of each suspension parameter on train performance is assessed.

Figure 2. Dynamic modeling structure of the train and its suspension system.

Table 1

VIP values of all suspension parameters

Symbol	Parameter name	Unit	Original value	Value range	VIP($$ R_h $$)	VIP($$ R_v $$)
Bold values indicate $$VIP > 1.0$$. Underlined values indicate additionally selected parameters. VIP: Variable importance in projection.
$$ x_1^{o} $$	Lateral stiffness of the primary steel spring	MN/m	0.98	0.10–9.8	0.0252	0.15303
$$ x_2^{o} $$	Longitudinal stiffness of the primary steel spring	MN/m	0.98	0.10–9.8	0.55299	0.34093
$$ x_3^{o} $$	Vertical stiffness of the primary steel spring	MN/m	4.52	0.45–45.2	1.5818	3.17028
$$ x_4^{o} $$	Lateral stiffness of the primary axle box arm joint	MN/m	13.70	1.37–137.0	0.26975	0.03386
$$ x_5^{o} $$	Longitudinal stiffness of the primary axle box arm joint	MN/m	5.52	0.55–55.2	2.94629	0.52944
$$ x_6^{o} $$	Vertical stiffness of the primary axle box arm joint	MN/m	5.52	0.55–55.2	0.02464	0.20156
$$ x_7^{o} $$	Stiffness of the primary vertical damper joint	MN/m	0.62	0.06–6.2	0.51207	0.70251
$$ x_8^{o} $$	Damping of the primary vertical damper	MN/m	0.01	0.00–0.1	0.10559	0.35821
$$ x_9^{o} $$	Lateral stiffness of the secondary air spring	MN/m	0.17	0.02–1.7	0.74644	0.21934
$$ x_{10}^{o} $$	Longitudinal stiffness of the secondary air spring	MN/m	0.17	0.02–1.7	0.45688	0.15309
$$ x_{11}^{o} $$	Vertical stiffness of the secondary air spring	MN/m	12.32	1.23–123.2	0.18619	1.62243
$$ x_{12}^{o} $$	Stiffness of the secondary lateral damper joint	MN/m	17.20	1.72–172.0	0.00487	0.08812
$$ x_{13}^{o} $$	Stiffness of the secondary antihunting damper joint	MN/m	8.82	0.88–88.2	1.13243	0.24706
$$ x_{14}^{o} $$	Longitudinal stiffness of the secondary traction rod	MN/m	7.88	0.79–78.8	0.28475	0.30521

The variable importance in projection (VIP) method is applied to identify critical suspension parameters based on their contributions to carbody vibration acceleration through partial least squares (PLS) regression analysis^[30]. The VIP value ($$ V_j $$) of the $$ j $$-th input variable is expressed as:

(1) $$ V_j = \sqrt{\frac{14}{\sum _{a=1}^{m} r_a^2(Y;t_a)} \sum\limits_{a=1}^{m} r_a^2(Y;t_a) \cdot w_{ja}^2} $$

where $$ m $$ is the number of extracted principal components, $$ Y $$ is the output variable, $$ w_{ja} $$ denotes the weight of input variable $$ x_j $$ on principal component $$ t_a $$, and $$ r_a(Y; t_a) $$ represents the correlation coefficient between $$ Y $$ and $$ t_a $$. Parameters with $$ V_j > 1 $$ are regarded as important driving factors in the model.

To ensure accurate VIP calculation, parameter ranges are determined based on the initial suspension parameters and historical data from heavy-haul railways, as summarized in Table 1. Parameter sampling is performed using the multiobjective Latin hypercube design (MOLHD) method to reduce the risk of local optima while respecting computational resource constraints. This method employs an improved iterative local search (ILS) algorithm to optimize both sample uniformity and orthogonality. Several optimization criteria are combined to guarantee the representativeness of the samples^[30]. In total, 200 sets of suspension parameter values are generated through MOLHD, and the train multibody dynamics model is subsequently evaluated on these samples to obtain the root-mean-square values of lateral ($$ R_h $$) and vertical ($$ R_v $$) vibration acceleration for VIP computation.

Based on the VIP results in Table 1, six parameters are retained for surrogate-model-based identification to capture their influence on both lateral and vertical vibration responses. Four parameters with VIP values exceeding 1.0 are selected: vertical stiffness of the primary steel spring ($$ x_3^{o} $$; VIP $$ =1.58 $$ for $$ R_h $$, $$ 3.17 $$ for $$ R_v $$), longitudinal stiffness of the primary axle-box arm joint ($$ x_5^{o} $$; VIP $$ =2.95 $$ for $$ R_h $$), vertical stiffness of the secondary air spring ($$ x_{11}^{o} $$; VIP $$ =1.62 $$ for $$ R_v $$), and stiffness of the secondary antihunting damper joint ($$ x_{13}^{o} $$; VIP $$ =1.13 $$ for $$ R_h $$). In addition, two parameters with the highest VIP values among the remaining candidates are retained, namely stiffness of the primary vertical damper joint ($$ x_7^{o} $$) and lateral stiffness of the secondary air spring ($$ x_9^{o} $$). This selection ensures that both primary and secondary suspension dynamics are comprehensively covered in multiobjective parameter identification.

The six selected parameters are renumbered as $$ x_1 $$ to $$ x_6 $$ for surrogate model training, as listed in Table 2.

Table 2

Selected key suspension parameter ranges for surrogate model training

Symbol	Parameter name	Unit	Nominal value	Sampling range	Original symbol
$$ x_1 $$	Vertical stiffness of the primary steel spring	MN/m	4.52	2.26–9.04	$$ x_3^{o} $$
$$ x_2 $$	Longitudinal stiffness of the primary axle-box arm joint	MN/m	5.52	2.76–11.04	$$ x_5^{o} $$
$$ x_3 $$	Stiffness of the primary vertical damper joint	MN/m	0.62	0.31–1.24	$$ x_7^{o} $$
$$ x_4 $$	Lateral stiffness of the secondary air spring	MN/m	0.17	0.085–0.34	$$ x_9^{o} $$
$$ x_5 $$	Vertical stiffness of the secondary air spring	MN/m	12.32	6.16–24.64	$$ x_{11}^{o} $$
$$ x_6 $$	Stiffness of the secondary antihunting damper joint	MN/m	8.82	4.41–17.64	$$ x_{13}^{o} $$

3.2. Suspension parameter optimization based on NSGA-II

Based on the constructed AO-GRBF surrogate model, a mapping is established from the suspension parameters $$ X \in \mathbb{R}^{6 \times 1} $$ to the predicted vibration responses $$ \hat{Y}(X)=\big[\hat{Y}_h(X), \, \hat{Y}_v(X)\big]^T $$, where $$ \hat{Y}_h(X) $$ and $$ \hat{Y}_v(X) $$ denote the lateral and vertical responses, respectively. Subsequently, NSGA-II is formulated to identify the optimal parameter vector by minimizing the deviation and the correlation loss between the predicted and measured responses ^[32,33].

The multiobjective optimization problem is formulated as:

where $$ X_{\min} $$ and $$ X_{\max} $$ are bounds specified in Table 2.

The correlation coefficients and mean squared errors are computed as:

where $$ \hat{y}_{\mathrm{h/v}, i} $$ and $$ y_{\mathrm{h/v}, i} $$ are predicted and measured accelerations, $$ \mu $$ denotes mean values, and $$ T $$ is the total number of samples.

Based on the optimization objectives defined above, the NSGA-II algorithm performs parameter identification through the following steps:

(1) Population Initialization: An initial population of size $$ M $$ is randomly generated as:

where $$ \mathbf{x}^{(0)}_i \in \mathbb{R}^{6} $$ represents the $$ i $$-th individual (suspension parameter vector) in the initial generation.

(2) Evolutionary Operations: For each generation $$ t = 0, 1, 2, \dots $$, each individual $$ \mathbf{x}^{(t)}_i $$ in the current population $$ \mathbf{P}^{(t)} $$ is first evaluated with the AO-GRBF surrogate to compute the four objective functions defined in Equation (10). Next, parent pairs $$ (\mathbf{x}^{(t)}_{p_1}, \mathbf{x}^{(t)}_{p_2}) $$ are selected via tournament selection to generate the offspring population $$ \mathbf{Q}^{(t)}=\{\mathbf{x}^{o}_1, \mathbf{x}^{o}_2, \dots, \mathbf{x}^{o}_M\} $$, and simulated binary crossover (SBX) is applied to produce the offspring:

where $$ u \sim U(0, 1) $$ is a uniform random variable, $$ \eta_c $$ is the crossover distribution index controlling solution spread, and $$ \beta_q $$ is the crossover parameter determining offspring proximity to parents. Polynomial mutation is subsequently applied to each offspring $$ \mathbf{x}^{o}_j $$ with probability $$ p_m $$:

where $$ x^{o}_{j, k} $$ is the $$ k $$-th component of the $$ j $$-th offspring, $$ \delta_k $$ follows the polynomial mutation distribution with index $$ \eta_m $$, and $$ (x_{k}^{\max}, x_{k}^{\min}) $$ are the bounds for the $$ k $$-th parameter.

The parent and offspring populations are merged to form $$ \mathbf{R}^{(t)}=\mathbf{P}^{(t)} \cup \mathbf{Q}^{(t)} $$, so that $$ |\mathbf{R}^{(t)}|=2M $$. Non-dominated sorting is then applied to obtain the Pareto fronts $$ \{F_1, F_2, \dots, F_K\} $$, where $$ F_1 $$ contains the best solutions. For each front $$ F_l $$, the crowding distance of a solution $$ \mathbf{x}_i $$ is computed as:

where $$ n_f = 4 $$ is the number of objective functions, $$ f_m^{(i)} $$ is the $$ m $$-th objective value of the $$ i $$-th solution, and $$ (f_m^{\max}, f_m^{\min}) $$ are the maximum and minimum values of the $$ m $$-th objective in front $$ F_l $$. The next-generation population $$ \mathbf{P}^{(t+1)} $$ is formed by selecting the best $$ M $$ individuals based on Pareto rank and crowding distance.

(3) Termination: The algorithm terminates upon convergence, yielding the final Pareto-optimal parameter set:

The integration of AO-GRBF and NSGA-II forms a closed-loop identification process, as illustrated in Figure 3. At each iteration, deviations between the predicted $$ \hat{Y}_i $$ and measured $$ Y_i $$ responses are evaluated to update the parameters $$ X_{k+1} $$, which are then fed back to the surrogate model until convergence is reached.

Figure 3. Closed-loop diagram integrating the AO-GRBF surrogate model and NSGA-II for parameter identification. AO-GRBF: Adaptively optimized Gaussian radial basis function; NSGA-II: nondominated sorting genetic algorithm II.

4.1. On-site data collection

The experimental dataset was collected along a heavy-haul railway line in China. Two categories of field measurements were obtained: (1) track-geometry irregularities (vertical profile, alignment, cross-level, superelevation, and gauge); and (2) carbody vibration accelerations. To ensure sufficient diversity in operating conditions for model training and evaluation, the dataset comprises multiple track segments that together cover small-radius curves, tangent sections, ascending grades, and descending grades.

The track irregularity data were recorded in real time using a track inspection system equipped with inertial navigation and vibration compensation modules, as shown in Figure 4. These data were subsequently used to construct a simulation environment that replicates actual operating conditions. In addition, carbody vibration accelerations were measured by accelerometers installed near the rear wheelset of the bogie, and all acceleration signals were sampled at a frequency of 100 Hz.

Figure 4. Track inspection system adopted for field measurements (photographs by the authors).

4.2. Result analysis

To evaluate the effectiveness of the proposed AO-GRBF-based identification framework, comparative experiments were conducted against the traditional G-RBF method, as well as recent approaches based on RBFNN and RBF-HDMR ^[10,35,36], and deep learning baselines including LSTM and Transformer networks ^[37,38]. All methods were executed on the same dataset generated by the SIMPACK model and under identical NSGA-II configurations to ensure a fair comparison.

The identified suspension parameter sets are summarized in Table 4, and the corresponding objective values predicted by the surrogates are provided in Table 5. Because some surrogates exhibit limited fitting accuracy and a propensity for local optima, the predicted results may deviate from the true system behavior. To validate the identification results, the optimal parameter sets in Table 4 were applied to a high-fidelity mechanistic train model implemented in SIMPACK. Furthermore, a 200 m operating subsegment from the field runs was selected to provide a common test condition. Simulations under this condition were conducted to generate vibration responses, enabling waveform-level and metric-level comparisons across methods.

The lateral and vertical vibration accelerations before and after parameter identification are compared in Figures 5 and 6. The predicted responses (blue) show varying degrees of agreement with the measured data (black), compared with the unidentified baseline (black dashed line), thereby confirming the effectiveness of surrogate-assisted parameter identification. As shown in Figures 5A and 6A, both the AO-GRBF method (red) and G-RBF (blue) improve waveform alignment, with AO-GRBF achieving closer amplitude agreement and more accurate phase matching.

Figure 5. Comparison of lateral vibration acceleration before and after parameter identification by various methods. (A) Comparison with GRBF; (B) Comparison with RBFNN; (C) Comparison with RBF-HDMR; (D) Comparison with LSTM; (E) Comparison with Transformer. All subfigures show the measured data (black) and the unidentified baseline (black dashed) as references. The proposed method (red) and other approaches exhibit different levels of consistency with the experimental observations. GRBF: Gaussian radial basis function; RBFNN: radial basis function neural network; RBF-HDMR: radial basis function high-dimensional model representation; LSTM: long short-term memory.

Figure 6. Comparison of vertical vibration acceleration before and after parameter identification by various methods. (A) Comparison with GRBF; (B) Comparison with RBFNN; (C) Comparison with RBF-HDMR; (D) Comparison with LSTM; (E) Comparison with Transformer. All subfigures show the measured data (black) and the unidentified baseline (black dashed) as references. The proposed method (red) and other approaches show different levels of consistency with experimental observations. GRBF: Gaussian radial basis function; RBFNN: radial basis function neural network; RBF-HDMR: radial basis function high-dimensional model representation; LSTM: long short-term memory.

Further comparisons with RBFNN, RBF-HDMR, LSTM, and Transformer methods are presented in Figure 5B-E and 6B-E. Among these methods, RBF-HDMR (blue) generates vibration responses that closely approximate the measured data (black), as illustrated in Figures 5C and 6C. Transformer networks (blue) also provide improved overall waveform alignment, as shown in Figures 5E and 6E. Consequently, the responses produced by RBF-HDMR (blue), Transformer (blue), and the proposed AO-GRBF method (red) are closely aligned with the measured data, and visual inspection alone does not clearly indicate which method achieves the best identification performance.

To quantitatively assess performance, the Pearson correlation coefficients ($$ p_h $$, $$ p_v $$) and mean squared errors ($$ m_h $$, $$ m_v $$) were calculated based on the actual vibration responses, as shown in Table 6. Compared with the GRBF approach, AO-GRBF achieves improvements of 78.41% and 18.55% in $$ p_h $$ and $$ p_v $$, respectively, while substantially reducing $$ m_h $$ and $$ m_v $$. These improvements can be attributed to its adaptive identification capability, which overcomes the limitations of fixed kernel parameters in traditional RBF implementations.

Relative to RBF-HDMR, which adopts a dimensional decomposition strategy, the proposed method achieves additional improvements of 3.10% in $$ p_h $$ and 7.46% in $$ p_v $$, along with reductions in mean squared error (MSE). In comparison with deep learning approaches, the performance differences are more evident. When evaluated against LSTM networks, the AO-GRBF model increases $$ p_h $$ and $$ p_v $$ by 14.24% and 27.14%, respectively, while reducing both error metrics. This outcome is mainly attributed to the tendency of LSTM models to overfit when trained on limited data. By contrast, Transformer networks achieve a relatively high vertical correlation ($$ p_v = 0.88790 $$), owing to their attention mechanisms. Nevertheless, the AO-GRBF model maintains improvements of 7.74% for $$ p_h $$ and 2.85% for $$ p_v $$, along with further reductions in MSE. Overall, these gains are attributed to the integration of NSGA-II-based hyperparameter tuning and an error-aware adaptive sampling strategy, which together improve global generalization and local approximation accuracy.

4.3. Computational efficiency analysis

To evaluate real-time feasibility for practical suspension parameter identification, a computational efficiency analysis was conducted. Experiments were performed on a workstation equipped with an Intel Core i9-14900KF CPU, 64GB RAM, and an NVIDIA GeForce RTX 4060 Ti. A high-fidelity mechanistic train model was constructed on the SIMPACK 2021x/SIMPACK Post 2021x platform. Mechanism-based identification with NSGA-II was realized by coupling MATLAB R2021a with SIMPACK 2021x, whereas surrogate construction (GRBF, RBFNN, RBF-HDMR, LSTM, Transformer, AO-GRBF) and surrogate-based identification with NSGA-II were implemented in MATLAB R2021a.

For a fair comparison, the total runtime is defined separately for the two categories of methods. For the high-fidelity mechanistic approach, total runtime denotes the elapsed time of a complete identification run, in which NSGA-II, implemented in MATLAB R2021a, directly executes the SIMPACK-based train model. For surrogate-based approaches, total runtime equals the sum of surrogate construction time in MATLAB R2021a (including training and tuning) and the subsequent NSGA-II search performed on the surrogate.

Parameter identification with the proposed AO-GRBF method is completed in 1.25 min [Table 7], yielding a speedup of approximately 3, 010$$ \times $$ relative to the high-fidelity simulation time of 3, 762 min. The improvement arises because surrogate evaluations are executed in milliseconds during the search, whereas each run of the SIMPACK-based mechanistic model requires about 1 min, and a complete identification typically involves well over 1, 000 model runs.

Although classical surrogates (GRBF, RBFNN) complete the process more quickly due to their simpler structures, their predictive accuracy is lower [Table 6]. Deep learning-based surrogates (LSTM, Transformer) require more than 4 min, which is about 3.3 times the runtime of the proposed method. Overall, the AO-GRBF approach achieves a favorable balance between computational efficiency and estimation accuracy, enabling time-critical suspension parameter identification and supporting rapid detection of stiffness and damping changes for proactive safety management.

4.4. Method stability analysis

To evaluate the robustness of the proposed AO-GRBF framework under repeated surrogate retraining, a stability analysis was conducted across 20 independent trials. In each trial, the training and validation datasets $$ \mathcal{D}_{\mathrm{train}} $$ and $$ \mathcal{D}_{\mathrm{val}} $$ were regenerated within the feasible parameter domain defined in Table 2. The surrogate model was retrained to incorporate error-guided adaptive sampling and NSGA-II optimization. Subsequently, the optimal suspension parameter vector was identified by minimizing the deviation between surrogate predictions and measured vibration responses.

The distributions of the identified suspension parameters are shown in Figure 7. All parameters exhibit compact clustering, and the coefficients of variation are below 9%, indicating consistent identification performance across trials. In particular, $$ x_1 $$, $$ x_2 $$, $$ x_3 $$, and $$ x_4 $$ cluster tightly around their baseline values, confirming stable identification results. In comparison, $$ x_5 $$ displays the widest yet limited spread, whereas $$ x_6 $$ shows moderate dispersion. Taken together, these findings demonstrate that the optimization process converges reliably to stable parameter combinations, even when the training datasets are regenerated independently.

Figure 7. Distributions of the identified suspension parameters $$ x_1 $$ – $$ x_6 $$ across 20 independent trials (units: N/m). Each violin plot illustrates the range, density, and clustering of parameter values, indicating consistent convergence across different surrogate retraining instances.

To quantify stability, performance deviations were computed relative to the baseline results in Section 4.2. For each trial, the identified parameters were evaluated using the high-fidelity SIMPACK model to obtain the correlation coefficients $$ p_v $$ and $$ p_h $$. Deviations were defined as $$ \Delta p_v = p_v - p_v^{\ast} $$ and $$ \Delta p_h = p_h - p_h^{\ast} $$, where the asterisks denote baseline values. The relative parameter deviation magnitude (relative $$ \|\Delta x\|_2 $$) was also calculated as the normalized Euclidean distance from the baseline parameter vector.

The relationship between performance accuracy and the magnitude of parameter deviation is illustrated in Figure 8. For vertical vibration correlation, $$ \Delta p_v $$ is bounded within $$ \pm 0.025 $$ across all deviation levels, and most values are concentrated near zero. This indicates that vertical performance is relatively insensitive to moderate parameter variations, consistent with the smooth behavior of vertical vehicle dynamics discussed in Section 3.1.3. By contrast, the lateral correlation metric $$ \Delta p_h $$ exhibits greater variability. Although it remains within $$ \pm 0.025 $$, no systematic performance decline is observed as parameter deviations increase. These results suggest that, although lateral dynamics are more sensitive to parameter changes, the identification framework maintains acceptable accuracy across the explored parameter space.

Figure 8. Performance deviations versus relative parameter deviation magnitude. (A) Vertical correlation deviation $$ \Delta p_v $$ versus relative $$ \|\Delta x\|_2 $$, showing limited sensitivity to parameter perturbations; (B) Lateral correlation deviation $$ \Delta p_h $$ versus relative $$ \|\Delta x\|_2 $$, exhibiting higher variability but without systematic degradation trends.

The trial-wise performance variations are depicted in Figure 9. Both $$ \Delta p_v $$ and $$ \Delta p_h $$ fluctuate around zero, with standard deviations of approximately 0.010 and 0.012, respectively. The interquartile ranges for both metrics are below 0.015, confirming that most trials achieve performance within 1.5% of the baseline. Importantly, larger parameter deviations (indicated by marker size) do not consistently correspond to greater performance losses. This shows that the method delivers robust performance across different parameter identification outcomes. The horizontal distributions around zero further validate the stability of the AO-GRBF framework under repeated surrogate retraining.

Figure 9. Trial-wise performance deviations across 20 independent trials. (A) Vertical correlation deviation $$ \Delta p_v $$ by trial number; (B) Lateral correlation deviation $$ \Delta p_h $$ by trial number. Marker sizes are proportional to relative parameter deviation magnitude $$ \|\Delta x\|_2 $$, showing no systematic relationship between parameter drift and performance loss.

To more comprehensively assess model robustness, an error space analysis is conducted to quantify absolute accuracy relative to the identification objective. Figure 10 reports trial-wise deviations of the mean squared acceleration error, $$ \Delta m_v $$ and $$ \Delta m_h $$, across 20 independent surrogate retraining runs conducted under the same protocol. Additionally, the correlation stability is confirmed in Figure 9 to be accompanied by bounded variations in absolute error. Over the 20 trials, $$ \Delta m_v $$ lies within $$ [-9.5\times 10^{-5}, \, 1.20\times 10^{-4}] $$ and $$ \Delta m_h $$ within $$ [-9.7\times 10^{-5}, \, 1.66\times 10^{-4}] $$, with interquartile ranges of approximately $$ 6.0\times 10^{-5} $$ and $$ 1.2\times 10^{-4} $$, respectively. Furthermore, marker sizes are used to encode the relative parameter deviation magnitude $$ \|\Delta x\|_2 $$, and no monotonic association is observed with $$ \Delta m_v $$ or $$ \Delta m_h $$; the points are scattered around zero without systematic bias. Taken together with Figure 9, these bounded and non-monotonic patterns indicate that the AO-GRBF framework maintains stable and robust accuracy under repeated surrogate retraining and moderate parameter drift.

Figure 10. Trial-wise deviations of error-type metrics relative to the baseline. (A) Deviation in the mean squared error of vertical acceleration ($$ \Delta m_v $$) by trial number; (B) Deviation in the mean squared error of lateral acceleration ($$ \Delta m_h $$) by trial number. Marker sizes are proportional to the relative parameter deviation magnitude $$ \|\Delta x\|_2 $$.

The stability analysis indicates that consistent identification performance is achieved across independent trials, with deviations remaining within acceptable engineering tolerances. This robustness is attributed to the combined use of error-guided adaptive sampling and NSGA-II optimization, which ensures reliable convergence despite dataset variability typical of practical applications.

Authors’ contributions

Experimental analysis and manuscript writing: Jiang, S.

Overall study design, conceptual framework development, and methodological guidance: Liu, J.

Technical support: Yang, S. X.

Availability of data and materials

Restrictions apply to the availability of these data. Field measurements were obtained from a heavy-haul railway line using proprietary inspection equipment under a collaboration agreement with a commercial railway operator (third party); therefore, public release is not permitted. Data may be made available from the corresponding author upon reasonable request and with prior written permission from the data owner, subject to confidentiality and data use agreements.

AI and AI-assisted tools statement

Not applicable.

Financial support and sponsorship

This work was supported in part by the National Key R & D Program of China (Grant No. 2021YFF0501101), the National Natural Science Foundation of China (Grant No. 52272347), the Natural Science Foundation of Hunan Province (Grant No. 2024JJ7132), and the Scientific Research Project of Hunan Provincial Department of Education (Grant No. 22A0391).

Conflicts of interest

Yang, S. X. is the Editor-in-Chief of the journal Intelligence & Robotics. He is not involved in any steps of editorial processing, notably including reviewers’ selection, manuscript handling and decision making. The other authors declare that there are no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2026.