2.1 Modelling the inverse problem of structural dynamics

In the structural design and verification of modern aircraft, FEMs are extensively employed for structural dynamic analysis, modal characteristic prediction, flutter boundary assessment and aeroelastic optimisation [Reference Byun and Jun17–Reference Wright and Cooper19]. The structural dynamics problem is typically formulated as a second-order linear system [Reference Noh, Hwang and Valluzzi20]:

(1) \begin{align} M\ddot u(t) + C\dot u(t) + Ku(t) = f(t)\end{align}

where $M,C,K \in {R^{n \times n}}$ denote the system mass, damping and stiffness matrices, respectively, while $u(t)$ and $f(t)$ represent the displacement response and external excitation vectors. In modal analysis, neglecting damping and external forcing, the natural frequencies ${\omega _i}$ and mode shapes ${\varphi _i}$ are obtained by solving the generalised eigenvalue problem:

(2) \begin{align} \big(K - \omega _i^2M\big){\varphi _i} = 0\end{align}

Mathematically, model updating constitutes a measurement-based parameter identification inverse problem. To achieve a balanced update across all modes, the optimisation is formulated as a multi-objective problem rather than a single scalar minimisation. The non-dominated sorting genetic algorithm II (NSGA-II) is employed to simultaneously minimise two distinct objective functions: the mean relative error across all considered modes ( ${f_1}$ ) and the maximum single-mode relative error ( ${f_2}$ ).

(3) \begin{align} \left\{ \begin{array}{l}{f_1}(x) = \dfrac{1}{n}\sum\limits_{i = 1}^n {\left| {\dfrac{{\lambda _i^c(x) - \lambda _i^e}}{{\lambda _i^e}}} \right|} \\[10pt]{f_2}(x) = {\max _{i \in [1,n]}}\left| {\dfrac{{\lambda _i^c(x) - \lambda _i^e}}{{\lambda _i^e}}} \right|\end{array} \right.\end{align}

where $\Omega $ denotes the feasible design domain, $\omega _i^{FE}$ and $\omega _i^{Exp}$ represent the i-th computed and experimental modal frequencies, respectively, and ${w_i}$ is the weighting coefficient reflecting measurement confidence. Since K and M typically exhibit complex nonlinear dependencies on x, this optimisation problem is inherently non-convex and multi-modal, rendering direct solution methods susceptible to local optima entrapment. By decoupling these objectives, the algorithm prevents the ‘smoothing effect’ where significant local errors in individual modes are masked by a low overall average error.

In the context of multi-objective optimisation for model updating, the ‘Pareto front’ represents a set of optimal compromise solutions where no single objective (e.g. frequency error of a specific mode) can be improved without degrading another. Unlike traditional single-objective optimisation that seeks one global minimum, identifying the Pareto front allows engineers to balance the trade-offs between global mean error and local single-mode deviations, thereby providing a more robust physical insight into the structure’s behaviour.

2.2 Construction of Bayesian regularisation neural network surrogate model

Given the prohibitive computational cost of directly invoking the finite element solver during iterative optimisation, this study employs an MLP to construct a surrogate model approximating the nonlinear mapping from parameter space to response space:

(4) \begin{align} {\textbf{Y}} = {\mathcal{F}_{NN}}({\textbf{X}},{\textbf{W}},{\textbf{b}}NN)\end{align}

where ${\textbf{X}}$ denotes the input design variables, ${\textbf{Y}}$ represents the predicted modal frequency response and ${\textbf{W}}$ and ${\textbf{b}}NN$ are the network weight matrices and bias vectors, respectively. The MLP neural network architecture is shown in Fig. 2.

Figure 2.

Schematic of MLP neural network architecture.

Although standard backpropagation algorithms perform effectively with abundant data, FEMU scenarios typically constrain training set size due to high sample acquisition costs. To mitigate overfitting under limited training data and enhance generalisation capability, the Bayesian regularisation (BR) algorithm is employed for network training.

Unlike conventional approaches that solely minimise mean square error ${E_D}$ , the BR algorithm incorporates a weight decay term ${E_W}$ , reformulating the objective function as:

(5) \begin{align} F({\textbf{W}}) = \beta {E_D} + \alpha {E_W} = \beta \sum\limits_{k = 1}^N {{{({y_k} - {t_k})}^2}} + \alpha \sum\limits_{j = 1}^M {w_j^2} \end{align}

where $\alpha $ and $\beta $ are hyperparameters governing the trade-off between prior probability and likelihood function. Within the Bayesian framework, network weights are treated as random variables following Gaussian distributions. Through posterior probability density maximisation, the algorithm adaptively updates $\alpha $ and $\beta $ during training, automatically penalising excessive weight magnitudes. This mechanism effectively reduces model complexity, enhancing robustness on unseen data and making it particularly suitable for high-dimensional nonlinear mapping challenges.

2.3 Two-stage collaborative optimisation strategy based on trust region resampling

Addressing the precision bottleneck frequently encountered by conventional single-stage global surrogate models during late-stage convergence, this paper designs a serial collaborative optimisation mechanism comprising global exploration and local refinement. The strategy constructs a locally refined response surface to approximate the FE-evaluated response in the neighbourhood of the candidate optimum, thereby reducing prediction errors arising from the smoothing effect of global models in small local regions.

During the initial updating stage, the design space is extensive and exhibits strong nonlinearity. The primary objective is to rapidly identify the Pareto front containing the global optimum rather than pinpointing a single extremum. Within the initial broad design space ${D_{global}} \in [{x_{lb}},{x_{ub}}]$ , Latin hypercube sampling (LHS) generates the initial sample set S₁. Following finite element analysis (FEA) response evaluation, the global Bayesian neural network model ${M_{global}}$ is trained on S₁. Subsequently, employing ${M_{global}}$ as the objective function evaluator, NSGA-II performs global search through non-dominated sorting and crowding distance mechanisms, yielding the first-stage Pareto non-dominated solution set P₁. To balance full-spectrum modal trends with individual mode accuracy – a critical engineering requirement – a ‘knee point’ strategy selects a compromise solution $x_{opt}^{(1)}$ exhibiting balanced multi-objective performance from P₁. This solution serves both as a candidate optimum and, more importantly, as the ‘anchor point’ guiding second-stage local search.

Although the first stage successfully identifies the approximate optimal solution location, training sample sparsity causes the global surrogate model ${M_{global}}$ to exhibit prediction errors near $x_{opt}^{(1)}$ that exceed the target modal-frequency agreement required for reliable updating. To address this limitation, a trust region resampling strategy is introduced: centred on the first-stage knee-point solution $x_{opt}^{(1)}$ , a contracted local hypercube defines the refined design space ${D_{local}}$ . This region’s boundaries are defined as:

(6) \begin{align} x_{lb}^{local} = x_{opt}^{(1)} \cdot (1 - \delta ),\quad x_{ub}^{local} = x_{opt}^{(1)} \cdot (1 + \delta )\end{align}

where $\delta $ denotes the contraction factor. In this study, $\delta = 0.20$ was selected after considering two competing requirements. The trust region must be sufficiently large to contain the true FE optimum if the first-stage knee point is affected by global-surrogate bias, but sufficiently small to increase the local sample density and reduce response nonlinearity. A ± 20% region was therefore used as a conservative local refinement range: it is half of the global geometric design range (±40%) and remains comparable to the allowable material-property range, while still providing adequate coverage for the second-stage FE re-evaluation.

Subsequently, small-scale LHS resampling is performed within the local trust region ${D_{local}}$ , with the FE solver invoked to compute FE-verified modal responses for these samples. A local neural-network surrogate model ${M_{local}}$ is then trained on this newly generated high-fidelity local dataset. Since the local response surface is smoother and the local sample density is substantially increased, the local model ${M_{local}}$ captures objective-function variations within the candidate optimum neighbourhood more accurately than the global surrogate. Employing the local model ${M_{local}}$ as the fitness evaluator, NSGA-II executes a second optimisation cycle. The resulting second-stage Pareto front P₂ is then interpreted together with direct FE re-analysis to avoid confusing surrogate-predicted optima with FE-evaluated updating performance. The final solution extracted from this front constitutes the updated physical parameters.

3.1 Description of the benchmark structure

The GARTEUR SM-AG19 physical model comprises primarily rectangular cross-section beams fabricated from aviation-grade aluminum alloy. Its geometric topology encompasses the fuselage, wings, vertical stabiliser, horizontal stabiliser and winglets. Figure 3(a) shows the overall structure of the aircraft, including the relative positions of its various components; Fig. 3(b) shows the specific dimensions of each component of the aircraft. According to the information in Fig. 3(b), the wingspan of the aircraft is 2,000 mm, the fuselage length is 1,500 mm and the tailplane span is 400 mm. The most notable dynamic characteristics of this structure include high flexibility due to the high aspect ratio configuration and significant sensitivity of the joint stiffness at multiple component connection points.

Figure 3.

Geometric dimensions of the GARTEUR SM-AG19 aircraft: (a) diagram of the aircraft components; and (b) schematic diagram of aircraft component dimensions.

The aircraft finite element model (FEM) was developed using HyperMesh software. Given the structure’s slender beam characteristics, fuselage and wing primary structures were discretised using CBEAM beam elements. To accurately represent physical connection relationships between components (e.g. wing-fuselage and tail-fuselage junctions), RBE2 multi-point constraint rigid elements modeled joint coupling.

To improve the physical representation of structural joints, the RBE2 rigid elements at the wing-fuselage and tail-fuselage junctions were replaced by CBUSH spring elements. The translational degrees of freedom (K1–K3) and torsional degree of freedom (K4) were assigned large stiffness values to preserve kinematic constraints, while the rotational stiffness components (K5, K6) were varied parametrically from 10⁴ to 10¹⁷ N·mm/rad. The test showed that modal frequencies converge when K5 and K6 exceed approximately 10⁹ N·mm/rad, with further increases producing frequency variations of less than 0.1%. Accordingly, K5 and K6 were fixed at 10¹⁰ N·mm/rad in the final model. This result indicates that under the current structural configuration, connection rotational stiffness is not the dominant source of model discrepancy, and the primary updating targets remain the beam cross-sectional parameters and component-level material properties. The finite element mesh diagram of GARTEUR SM-AG19 is shown in Fig. 4. The final FEM consists of 312 elements. Specifically, the fuselage and lifting surfaces are discretised into 308 CBEAM elements, while the connections are modeled using 2 RBE2 rigid elements and 2 CBUSH elements.

Figure 4.

Finite element mesh of the GARTEUR SM-AG19 aircraft.

The material properties assumed for the initial model are those of an isotropic aluminum alloy. Regarding the boundary conditions, the numerical modal analysis was performed under free-free conditions to match the setup of the GVT, ensuring that the rigid-body modes (first six modes) were separated from the elastic modes used in the updating process. Although the initial FEM was constructed from design drawings, inevitable discrepancies arise between eigenvalue analysis results and GVT data due to manufacturing tolerances, assembly clearances and beam element simplification errors, particularly within high-order and densely spaced modal frequency bands.

3.2 Parameter selection and network configuration

Based on structural mechanics principles and sensitivity analysis, variations in cross-sectional geometric parameters (width b and height h) directly alter the section moment of inertia and area, thereby modifying the structural stiffness matrix K and mass matrix M. Recognising that model errors primarily originate from machining tolerances, beam idealisation and equivalent stiffness deviations at connection interfaces, this study selected 16 design variables x for updating, encompassing rectangular cross-sectional dimensions of the fuselage, wings, winglets and tail, as well as material properties (Young’s modulus E and material density rho). The shear modulus G is no longer treated as an independent variable but is derived using a fixed Poisson’s ratio to maintain isotropic material consistency. To avoid repetitive listing, the variables are grouped by component in Table 1.

Table 1.

Definition of design variables for model updating

For the 16-dimensional input space and first 10 modal frequency outputs, the MLP neural network described in Section 2.2 was instantiated. Following hyperparameter optimisation, the final network topology was configured as follows: the input layer comprises 16 neurons corresponding to design variables X1–X16 in Table 1; the hidden layer adopts a two-layer architecture with 64 neurons per layer. The network weights W and biases b are initialised using the Xavier initialisation strategy to ensure variance consistency across layers and prevent gradient vanishing during early training epochs. The hidden layers employ the hyperbolic tangent (Tanh) activation function [Reference LeCun, Bottou, Orr, Müller, Montavon, Orr and Müller21] to capture strong nonlinear characteristics within the design space, while the output layer employs linear activation (Purelin) [Reference Goodfellow, Bengio and Courville22] to encompass the broad frequency response range.

3.3 Training and validation of two-stage surrogate models

To balance computational efficiency and local convergence accuracy, this study adopted a ‘Global Coarse Search – Local Refinement’ sequential modelling strategy. The specific sampling and training processes were executed in the following two stages.

3.3.1 Construction of global surrogate model

To capture the global mapping relationship between structural parameters and modal frequencies within the broad initial design space, a global surrogate model was first constructed. The variation range of each geometric design variable was set to ±40% of its initial nominal value, while the variation ranges of each elastic modulus and density were set to ±15%, defining the initial design space ${D_{global}}$ . LHS generated 1,500 spatially uniform input samples within this domain. The OptiStruct FE solver performed forward modal analysis on these 1,500 samples, extracting the first 10 natural frequencies as response labels. The dataset was randomly partitioned into training and test sets at an 85%:15% ratio. The global MLP model was trained using the Bayesian regularisation algorithm. Figure 5 presents regression analysis results for each modal frequency on the test set, and Fig. 6 presents the root mean square error (RMSE) of each vibration mode. The coefficient of determination (R²) exceeded 0.99 for all modes, confirming the global model’s adequacy for supporting NSGA-II preliminary Pareto-front identification and ROI localisation, but not for final FE-verified parameter selection.

Figure 5.

Scatter plot comparing predicted values and true values of modal frequencies.

Figure 6.

Global surrogate-model RMSE for each modal frequency.

To explicitly demonstrate the necessity of the second-stage local optimisation, Fig. 7 compares the surrogate-predicted Pareto solutions with their direct FE-evaluated responses. The blue points represent the optimal solutions predicted by the global MLP surrogate model. When these predicted optimal design variables are substituted back into the high-fidelity OptiStruct FE solver, the true structural responses (red points) exhibit a visible drift. Therefore, the first-stage Pareto front is used only to identify the ROI and define the trust-region centre, rather than being interpreted as the final FE-verified Pareto front.

Figure 7.

Pareto front comparison.

This discrepancy illustrates the residual epistemic uncertainty of the global surrogate model. Because the global MLP is trained across a broad design space (±40%), it inherently suffers from a smoothing effect and lacks the resolution to capture small, highly nonlinear gradient variations in the immediate vicinity of the FE-evaluated optimum. Consequently, while the global surrogate successfully narrows down the ROI, relying on it alone yields FE-reanalysed responses that do not fully satisfy the desired modal-frequency agreement. This mismatch between surrogate-predicted optimality and FE-evaluated performance makes the second-stage trust-region resampling strictly necessary. By confining the design space to a ± 20% neighbourhood around the knee point and generating new high-fidelity samples directly from the FE solver, the local model substantially reduces this prediction bias.

3.3.2 Trust region resampling and local modelling

Despite the global model’s satisfactory test-set accuracy, its prediction residual within the small neighbourhood of the candidate optimum may still exceed the accuracy required for reliable modal-frequency updating. Therefore, local high-fidelity modelling was implemented based on the first-stage knee-point solution.

Figure 8.

Scatter plot comparing predicted values and true values of modal frequencies.

Figure 9.

Local surrogate-model RMSE for each modal frequency after trust-region resampling.

Centred on the first-stage knee-point solution $x_{opt}^{(1)}$ , a local trust region ${D_{local}}$ with ±20% variable variation was established. This region substantially contracts the design space and dramatically reduces response nonlinearity while retaining sufficient coverage around the candidate optimum. LHS resampling within the local region ${D_{local}}$ generated 200 new samples. Critically, this step bypassed global-surrogate interpolation and directly invoked the FE solver for FE-verified calculations. This high-fidelity re-evaluation strategy reduces epistemic uncertainty inherited from the first stage. A local MLP model was trained on these 200 high-fidelity local samples. As shown in Fig. 8, the predicted modal frequencies of the local surrogate model agree well with the corresponding FE-calculated values, indicating that the locally refined model captures the parameter-frequency relationship within the trust region with high accuracy. Owing to the smoother local response surface and substantially increased sample density, the local surrogate model’s RMSE on the validation set decreased by an order of magnitude relative to the global model, providing a more accurate local response approximation for second-stage refined optimisation. The specific RMSE for each vibration mode is shown in Fig. 9.