\title{Comparison of surrogate strategies for damage-aware optimization of buckling-delayed shear-link dampers with adaptive finite element validation}
\author[1]{J. Irazábal}
\author[1,2]{J. Ramírez}
\author[1,2]{J. M. González}
\author[1]{L. Lázaro}
\author[1,2]{F. Rastellini}
\author[2,3]{G. Bozzo}
\author[4]{L. Bozzo}
\authormark{IRAZÁBAL \textsc{et al.}}
\titlemark{COMPARISON OF SURROGATE STRATEGIES FOR DAMAGE-AWARE OPTIMIZATION OF SHEAR-LINK DAMPERS}
\address[1]{\orgname{Centre Internacional de Metodes Numerics en Enginyeria (CIMNE)}, \orgaddress{\city{Barcelona}, \country{Spain}}}
\address[2]{\orgname{Universitat Politecnica de Catalunya (UPC)}, \orgaddress{\city{Barcelona}, \country{Spain}}}
\abstract[Abstract]{Buckling-delayed shear-link dampers are used in seismic-resistant structures as passive devices that concentrate energy dissipation while limiting damage to the primary system. Their geometric optimization requires balancing high energy dissipation with strict control of local damage. Finite element models can accurately reproduce the nonlinear cyclic response of these devices and provide internal quantities such as damage indicators and local distortion, but their computational cost prevents their direct use within iterative optimization loops. This work proposes an adaptive surrogate-assisted optimization framework for buckling-delayed shear-link dampers. First, experimentally calibrated nonlinear finite element models are used to generate reference datasets for dampers with different geometric and mechanical configurations. Supervised learning models are initially evaluated, with support vector regression and Gaussian process regression consistently providing high predictive accuracy. This suggests that kernel-based, distance-dependent approximations are well suited to the problem, motivating the introduction of radial basis function surrogates as a computationally efficient alternative. The surrogate predictions are coupled with a differential evolution algorithm through a damage-aware objective function that controls local damage and uses dissipated energy as a performance criterion. Optimized geometries are finally re-evaluated with finite element simulations. When the surrogate error exceeds the adopted tolerances, the new simulation result is added to the dataset and the surrogate models are retrained. The proposed framework provides an efficient damage-aware optimization of seismic energy dissipation devices.}
\keywords{Buckling-delayed shear link, seismic energy dissipation, FEM validation, surrogate modelling, machine learning, radial basis functions, Differential Evolution.}
\jnlcitation{\cname{%
\author{Irazábal J.},
\author{Ramírez J.},
\author{González J. M.},
\author{Lázaro L.},
\author{Rastellini F.},
\author{Bozzo G.}, and
\author{Bozzo L.}}.
\ctitle{Damage-aware surrogate optimization of buckling-delayed shear-link dampers with adaptive finite element validation.}\cjournal{\it Earthquake Engineering \& Structural Dynamics.}\cvol{2026;00(00):1--18}.}
\maketitle
\renewcommand\thefootnote{\fnsymbol{footnote}}
\setcounter{footnote}{1}
\section{Introduction}\label{sec:introduction}
Shear-link beam (SLB) dampers are passive energy dissipation devices widely used in seismic-resistant structures, designed to undergo stable inelastic deformations while limiting damage to the primary system. Their configuration concentrates inelastic demand in replaceable components, enabling high energy dissipation, ductile response and stable hysteretic behavior under severe cyclic loading \cite{Malley1984,Okazaki2007}. Within this family, buckling-delayed shear-link (BDSL) dampers incorporate a mechanical configuration that promotes shear-dominated behaviour while delaying local and global buckling.
The global hysteretic response of BDSL dampers can be characterized through experimental testing. However, laboratory campaigns present important limitations when addressing design optimization problems: internal state variables such as local plastic strains, stress triaxiality or damage indicators cannot be directly measured and the high cost and logistical complexity of experimental programs restrict the number of geometric configurations that can be explored. These limitations become even more critical when the design objective is not only to increase global energy dissipation and maximum displacement supported by the device, but also to control where and how damage develops within the device, making systematic optimization impractical.
In this context, finite element method (FEM) simulations provide a powerful framework for generating high-fidelity datasets across a wide range of geometric configurations, making them an ideal foundation for surrogate-based optimization strategies. Advanced nonlinear FEM models can accurately reproduce experimental cyclic behaviour, including plasticity, geometric nonlinearity, contact interactions, local instability and damage evolution, while providing access to both global response quantities and local indicators governing failure mechanisms. However, their high computational cost makes their direct use within optimization loops inefficient. This limitation motivates the use of surrogate models trained on FEM-generated data, which can approximate the structural response with reduced computational effort. The use of surrogate models enables the possibility of evaluating a large number of design configurations, facilitating systematic optimization.
Over the past decades, FEM has been widely used to study BDSL dampers and other passive seismic energy dissipation devices, providing detailed insight into nonlinear cyclic response, stiffness degradation and local inelastic mechanisms \cite{Deng2014a,Deng2015}. It has also supported the optimization and parametric analysis of these devices by enabling systematic exploration of geometric configurations and performance criteria under prescribed loading \cite{Deng2014,Deng2015a}. In this context, FEM-based studies have been extensively applied to characterize the mechanical response of metallic dampers. Motamedi et al. \cite{Motamedi2018} investigated accordion metallic dampers through combined experimental and numerical analyses, assessing the influence of key geometric variables on stiffness, strength and energy dissipation. Ghamari et al. \cite{Ghamari2021} studied I-shaped shear links in concentrically braced frames, while Xiong et al. \cite{Xiong2024} examined replaceable steel shear links with different short-length ratios, highlighting the strong influence of geometry on cyclic performance and failure modes. Simplified analytical and semi-empirical models have also been proposed to reduce computational cost \cite{Deng2014b} and more recent simulation-based studies have further explored the role of geometric and material variables in damper performance \cite{Kim2022}.
Geometric optimization has also been extensively explored. Zhang et al. \cite{Zhang2017} proposed a Kriging-assisted framework to maximize hysteretic energy in coupling beam dampers. Farzampour et al. \cite{Farzampour2019} optimized butterfly-shaped shear links by maximizing the ratio between dissipated energy and plastic strain, while Khatibinia et al. \cite{Khatibinia2019,Khatibinia2021} developed efficient strategies for U-shaped dampers using FEM and surrogate models. Shi et al. \cite{Shi2019} introduced a non-parametric shape optimization framework for shear panel dampers and Saleh et al. \cite{Saleh2024,Saleh2026} extended this line through topology optimization of shear-link configurations. More recent contributions include the hybrid cellular automata approach by Mendoza-Cuy et al. \cite{MendozaCuy2025} and the statistical optimization framework by Rios et al. \cite{Rios2025}.
Data-driven approaches have mainly focused on response or property prediction. Chan et al. \cite{Chan2015} used nonlinear autoregressive exogenous (NARX) models to reproduce hysteretic behaviour. Bae et al. \cite{Bae2020} developed models for low-cycle fatigue estimation and Almasabha et al. \cite{Almasabha2022} predicted the shear strength of short steel links using ML. Elgammal et al. \cite{Elgammal2024} modelled hysteretic restoring forces using data-driven approaches, while Hu et al. \cite{Hu2023} proposed explainable ML models for the probabilistic prediction of buckling stress. Physics-informed approaches have also been explored, such as the PINN framework proposed by Hu et al. \cite{Hu2022}.
All these works demonstrate the increasing interest in applying FEM-based and data-driven approaches, as well as in combining both, to analyse, understand and optimize seismic energy dissipation devices. However, most of these studies focus either on the prediction of the hysteretic response or on maximizing energy dissipation, leaving a critical aspect insufficiently explored: the need to control local damage while maintaining adequate dissipative capacity. In practice, excessive local damage may compromise structural integrity, reduce durability and lead to premature failure, even when global energy dissipation is improved.
The present work addresses this gap through a damage-aware surrogate-assisted optimization methodology in which the objective is not only to maximize distortion or energy dissipation, but also to balance dissipative performance with damage indicators derived from FEM simulations. The proposed approach combines: (i) experimentally calibrated nonlinear FEM models used as numerical ground truth; (ii) supervised surrogate models trained to predict local damage and distortion indicators; (iii) a Differential Evolution (DE) optimizer; and (iv) an adaptive FEM validation and retraining loop.
Figure \ref{fig:MethodologyFlowChart} summarizes the proposed workflow. The different stages of the methodology, together with the surrogate modelling, optimization strategy, validation procedure and corresponding results, are described in the following sections.
The BDSL dampers analysed in this work, with one representative configuration shown in Figure \ref{fig:Device}, are designed to concentrate energy dissipation in localized reduced-thickness zones, hereafter referred to as windows, while preserving the structural integrity of the surrounding frame. The optimization variables correspond to the window thicknesses, whereas the frame dimensions are kept fixed. The dissipative element is connected to a load-transfer system through a mechanism that allows imposed in-plane displacement while preventing axial force transmission, thereby promoting a shear-dominated response. Under cyclic loading, plastic deformation is intended to concentrate in the windows, whereas the frame provides load transfer and stability.
\caption{Representative BDSL damper configuration considered in the optimization.}
\label{fig:Device}
\end{figure}
This separation of functions leads to a non-trivial design problem. Thin windows may enhance ductility and dissipative activation, but they may also promote excessive damage localization. Conversely, thicker windows may increase strength while transferring inelastic demand to the frame. Since severe frame damage may compromise the structural integrity of the device, frame damage must be penalized more strongly than window damage. At the same time, the dissipative demand should be distributed as uniformly as possible among the windows, avoiding configurations in which a single window absorbs most of the deformation while the remaining windows stay underused. Consequently, the design problem cannot be reduced to maximizing force or total dissipated energy alone, but must also control where damage develops and how the windows participate in the dissipative process.
The design variables considered in this work are the window thicknesses
where $W$ denotes the number of windows. The width and height identifiers of the device are represented by $B$ and $H$, respectively. Five geometry families are considered, as shown in Figure \ref{fig:GeometryFamilies}: H30\_B29, H30\_B34, H45\_B29, H45\_B34 and H60\_B34. Devices with $H=30$ cm have two windows, those with $H=45$ cm have three windows and those with $H=60$ cm have five windows.
\caption{BDSL families considered for optimization in the current study. From left to right: H30\_B29, H30\_B34, H45\_B29, H45\_B34 and H60\_B34.}
\label{fig:GeometryFamilies}
\end{figure}
The admissible thickness ranges are defined according to the geometry family and manufacturing constraints. The main characteristics of each family are summarized in Table~\ref{tab:families}. In all cases, the frame thickness is kept constant at 30 mm.
\begin{table}[htbp]
\centering
\caption{Geometry families considered for optimization in the current study and admissible window thickness ranges.}
\label{tab:families}
\begin{tabular}{lllllll}
\toprule
Family & Height $H$& Width $B$& Windows & Frame thickness & Design variables & Thickness bounds \\
\midrule
H30\_B29 & 30 cm & 29 cm & 2 & 30 mm &$t_{w,1},t_{w,2}$& 10--20 mm \\
H30\_B34 & 30 cm & 34 cm & 2 & 30 mm &$t_{w,1},t_{w,2}$& 10--20 mm \\
H45\_B29 & 45 cm & 29 cm & 3 & 30 mm &$t_{w,1},t_{w,2},t_{w,3}$& 5--14 mm \\
H45\_B34 & 45 cm & 34 cm & 3 & 30 mm &$t_{w,1},t_{w,2},t_{w,3}$& 5--14 mm \\
H60\_B34 & 60 cm & 34 cm & 5 & 30 mm &$t_{w,1},\ldots,t_{w,5}$& 5--12 mm \\
\bottomrule
\end{tabular}
\end{table}
\section{Validation of the FEM numerical model}\label{sec:fem}
The surrogate models developed in this work were trained using data generated from three-dimensional FEM simulations. The numerical model is based on a previously calibrated and validated representation of the BDSL device, described in detail in Ramirez et al. \cite{RamirezMachado2025}. An example of the numerical setup is shown in Figure \ref{fig:FEMsetup}. The simulations were carried out using the COMPACK code, an explicit dynamic FEM solver based on an incremental formulation \cite{Martinez2011}. The model accounts for large displacements, material and geometric nonlinearities, contact interactions and the boundary conditions associated with the experimental configuration.
The dissipative steel component is modelled as ASTM A36 steel, whose cyclic plastic behaviour is represented by the Yoshida--Uemori model \cite{Yoshida2002,Jia2014}. This constitutive law allows the model to reproduce cyclic hardening, softening and Bauschinger-type effects under large plastic deformation. The steel component is discretized using linear eight-node hexahedral solid elements, providing a structured three-dimensional mesh suitable for extracting local stress, strain and damage-related fields.
The imposed displacement is applied through an actuator-like connector that transfers horizontal displacement while preventing axial load transmission, thereby reproducing the kinematic condition required to promote a shear-dominated response. Three experimental tests were performed following the loading protocols defined in the ANSI/AISC 341-16 \cite{American2002} and ASCE 7-22 \cite{American2017} qualification standards for seismic energy dissipation devices. Those protocols involve cyclic loading with progressively increasing amplitudes. Additional contact and confinement conditions are included to represent the experimental anti-buckling system and to control local and global instability.
\caption{FEM validation model of the BDSL device: mesh discretization, main components, boundary conditions and local/global buckling control.}
\label{fig:FEMsetup}
\end{figure}
The model was calibrated and validated against cyclic experimental tests performed on representative BDSL specimens. The calibration involved the material parameters, assembled geometry, contact definitions, support flexibility and boundary conditions. The validated model accurately reproduces the main global experimental responses. Figure \ref{fig:FEM_validation_comparison} shows the comparison between experimental and numerical results, confirming the suitability of the FEM model as a numerical reference for configurations beyond those experimentally tested.
\caption{Experimental--numerical validation of the BDSL model: hysteretic response (left) and cumulative dissipated energy (right).}
\label{fig:FEM_validation_comparison}
\end{figure}
Once validated, the FEM model is used to generate the datasets required for surrogate training and optimization. The optimization strategy relies on local damage indicators in the dissipative windows and surrounding frame, together with local distortion measures associated with the activation of the dissipative mechanism. Since these quantities are difficult to measure experimentally, the use of a high-fidelity FEM model provides valuable access to the internal state variables and local fields governing damage evolution and energy dissipation.
The damage indicator adopted in this work is the Triaxial Failure Damage Map (TFDMap) \cite{Rastellini2016}. This stress-triaxiality-based indicator evaluates the proximity of each material point to ductile failure by comparing its stress triaxiality and accumulated equivalent plastic strain with a reference failure envelope \cite{Rice1969,Bao2004,Wierzbicki2005,Bai2008}. In this study, the TFDMap is used as a post-processing damage-screening indicator, not as a constitutive fracture criterion. Its purpose is therefore not to explicitly predict crack initiation, but to compare geometrical configurations and ensure that optimized designs remain within acceptable damage levels.
\section{Dataset generation and surrogate modelling}\label{sec:surrogates}
\subsection{Design of experiments}\label{subsec:doe}
The FEM campaign is planned to cover the admissible parameter domain of each device family while ensuring a homogeneous exploration of the multidimensional space. The design variables are the window thicknesses $t_{w,i}$, whose combinations are generated through a Design of Experiments (DoE) strategy based on Latin Hypercube Sampling (LHS) optimized with the maximin criterion \cite{Joseph2008}. This approach provides a near-random yet space-filling distribution of samples, reducing clustering and improving the representation of the admissible domain.
For each geometry family, the optimization cycle starts from an initial set of FEM simulations selected to provide reasonable coverage of the design space while keeping the computational cost as low as possible. The number of initial samples is defined according to the dimensionality of each family: 8 samples for two-window devices, 16 samples for three-window devices and 64 samples for five-window devices. In all cases, the number of samples is selected as a power of two, which facilitates the potential use of Progressive Latin Hypercube Sampling (PLHS) \cite{Sheikholeslami2017} in future iterations. This would allow the DoE to be expanded while preserving its space-filling properties and avoiding the need to repeat previously computed simulations.
To improve surrogate robustness near the admissible limits, the sampling domain is extended slightly beyond the actual optimization ranges reported in Table~\ref{tab:families}. This reduces the risk of extrapolation when evaluating candidate solutions close to the true design limits. The extended thickness ranges and the number of initial simulations for each family used in the DoE are summarized in Table~\ref{tab:families_doe}.
\begin{table}[htbp]
\centering
\caption{Geometry families, window thickness ranges and number of simulations considered during the DoE generation.}
H30\_B29 & 30 cm & 29 cm & 2 & 30 mm &$t_{w,1},t_{w,2}$& 8--22 mm & 8 \\
H30\_B34 & 30 cm & 34 cm & 2 & 30 mm &$t_{w,1},t_{w,2}$& 8--22 mm & 8 \\
H45\_B29 & 45 cm & 29 cm & 3 & 30 mm &$t_{w,1},t_{w,2},t_{w,3}$& 4--16 mm & 16 \\
H45\_B34 & 45 cm & 34 cm & 3 & 30 mm &$t_{w,1},t_{w,2},t_{w,3}$& 4--16 mm & 16 \\
H60\_B34 & 60 cm & 34 cm & 5 & 30 mm &$t_{w,1},\ldots,t_{w,5}$& 4--14 mm & 64 \\
\bottomrule
\end{tabular}
\end{table}
For every sampled configuration, a FEM simulation is performed under a displacement-controlled cyclic loading protocol. Since the admissible deformation demand depends on the size of the device, different loading patterns are adopted according to the device height. As shown in Figure~\ref{fig:LoadPatterns}, the maximum displacement amplitude increase with the device height, consistently with the expected performance range of each geometry family.
\caption{Displacement-controlled cyclic loading patterns adopted for the different device heights considered in the FEM campaign.}
\label{fig:LoadPatterns}
\end{figure}
The resulting dataset stores the input variables and the structural response quantities extracted at the final time of each simulation. The target outputs include damage indicators in the windows and in the frame, together with local distortion measures associated with dissipative activation. For compactness in the optimization formulation, the aggregated TFDMap indicator in window $i$ is denoted by $\TFD_i$, whereas the corresponding frame indicator is denoted by $\TFD_f$. Both quantities are computed as the average TFDMap value of the 12 nodes with the highest values within each region. This allows to capture the most critical damage levels while reducing sensitivity to isolated numerical peaks.
The maximum local shear distortion in each window is denoted by $\varepsilon_{xy,i}$ and is used as an indicator of the energy dissipation capacity. When included in the objective function, each window contribution is weighted according to its effective geometric volume, so that the optimization accounts for both local strain intensity and the material volume involved in the dissipation process.
\subsection{Supervised ML surrogate models}\label{subsec:ml_models}
This work compares several supervised surrogate models for predicting FEM-derived damage and distortion indicators. The considered models cover three families: tree-based methods, including Random Forest (RF) \cite{Breiman2001}, Gradient Boosting Regression (GBR) \cite{Friedman2001} and XGBoost \cite{Chen2016}; kernel-based methods, including Support Vector Regression (SVR) \cite{Drucker1996} and Gaussian Process Regression (GPR) \cite{Williams1995}; and neural-network models, represented by the Multilayer Perceptron (MLP) \cite{Rosenblatt1958,Rumelhart1986}. RF relies on bootstrap aggregation of decision trees, GBR and XGBoost are sequential boosting approaches, SVR and GPR exploit kernel functions to model nonlinear relationships and MLP approximates nonlinear input--output mappings through interconnected layers.
For every geometry family, a separate regression model is trained for each target output. The input vector contains the window thicknesses of the corresponding device, while the outputs are the local distortion indicators in each window, $\varepsilon_{xy,i}$, the window damage indicators, $\TFD_i$, and the frame damage indicator, $\TFD_f$. Therefore, for a geometry family with $W$ windows, $2W+1$ surrogate models are trained: $W$ models for window distortion, $W$ models for window damage and one model for frame damage. For each output, all candidate algorithms are evaluated independently and the selected model is retained for the optimization stage.
Hyperparameters are optimized using Bayesian optimization \cite{Snoek2012}, with 40 evaluations per model and RMSE as the refit criterion. The first 10 evaluations are randomly sampled to explore the search space, while the remaining 30 are guided by the Bayesian surrogate model. This strategy provides a more efficient alternative to exhaustive grid search, particularly considering the number of geometry families, target outputs and candidate algorithms analysed. It also allows a broader exploration of continuous hyperparameter ranges. Table~\ref{tab:cv_hyperparameter_settings} summarizes the hyperparameter search spaces considered.
The cross-validation strategy is adapted to the dataset size. Leave-One-Out validation is used for $N\leq20$, repeated four-fold cross-validation with five repetitions for $21\leq N\leq80$, and shuffled five-fold cross-validation for larger datasets. For small datasets, the search spaces of tree-based models are additionally restricted to reduce overfitting.
\begin{table*}[ht!]
\centering
\caption{Hyperparameter search spaces used for Bayesian optimization of the supervised surrogate models.}
\label{tab:cv_hyperparameter_settings}
\scriptsize
\setlength{\tabcolsep}{5pt}
\renewcommand{\arraystretch}{1.10}
\begin{tabular}{llll}
\toprule
Model & Preprocessing / kernel & Hyperparameter & Search space \\
\parbox{0.95\textwidth}{\footnotesize\textit{Note:} For tree-based models, the search spaces are automatically restricted for small datasets to reduce overfitting.}
\end{table*}
Model selection is performed in two stages. First, for each candidate algorithm, Bayesian hyperparameter optimization is carried out using cross-validated Root Mean Squared Error (RMSE) as the refit criterion. The best hyperparameter configuration is therefore the one with the lowest mean RMSE. Then, the best configurations from all candidate algorithms are compared. The model with the lowest mean RMSE defines a competitive threshold and all models with an RMSE within 5\% of this value are retained. Among these competitive models, the final selection is based on the lowest relative RMSE dispersion, computed as the standard deviation of the fold-wise RMSE divided by the mean RMSE. If two models have the same dispersion, the one with the lower mean RMSE is preferred. This procedure, summarized in Figure \ref{fig:BayesianSearchCV}, favours surrogates that are both accurate and stable.
\caption{Workflow of the supervised surrogate training and selection strategy. For each output variable, the cross-validation strategy is adapted to the dataset size, Bayesian optimization is used to tune each candidate model and the final surrogate is selected according to RMSE accuracy and fold-wise RMSE dispersion.}
\label{fig:BayesianSearchCV}
\end{figure}
Preliminary executions of the proposed workflow show that SVR and GPR always provide the highest, or second-highest, predictive accuracy for the considered datasets. These results suggest that kernel-based models, and in particular distance-based similarity measures, are well suited to approximate the FEM response surfaces involved in this problem. This observation motivates the assessment of Radial Basis Function (RBF) interpolation \cite{Gutmann2001} as a simpler and computationally efficient surrogate alternative. Although surrogate evaluation is negligible compared with FEM simulations, the training and hyperparameter optimization of complex supervised models can still become relevant when several outputs, geometry families and adaptive iterations are considered. RBF interpolation provides a non-parametric and fast-to-train alternative that can capture nonlinear response surfaces, making it attractive for low-dimensional and moderately sampled design spaces.
The RBF surrogate is implemented using the interpolator available in the \textit{SciPy} library \cite{Virtanen2025}. For each output variable, the model is trained using the same input features employed by the supervised ML surrogates. In the current implementation, a multiquadric radial basis function is adopted with zero smoothing and automatic shape-parameter selection.
where $\mathbf{x}_j$ denotes the FEM-sampled geometries, $\lambda_j$ represents the interpolation weights and $\phi$ corresponds to the selected radial basis function.
For each output variable, a final RBF surrogate is trained using all available FEM samples in the current iteration and stored for subsequent use in the optimization process. Its predictive performance is assessed through Leave-One-Out validation: each FEM sample is iteratively removed from the dataset, a temporary RBF interpolant is trained with the remaining samples and the excluded sample is predicted. The resulting out-of-sample predictions are then used to compute RMSE, MAE and $R^2$ in order to assess the interpolation accuracy.
The proposed methodology seeks to minimize local damage in both the dissipative windows and the surrounding frame, while promoting a balanced contribution of all windows to the energy dissipation process. The geometric optimization is carried out using DE \cite{Storn1997}, a population-based global optimizer that does not require gradient information and is therefore suitable for nonlinear and non-convex surrogate response surfaces. In the current implementation, DE is run with a maximum of 500 iterations, a population size factor of 25 and a convergence tolerance of $10^{-6}$. Once an optimal candidate is obtained, an adaptive FEM validation loop is applied to verify the predicted geometry before acceptance.
For each candidate geometry $\mathbf{x}$, the trained surrogate models predict the window distortions $\hat{\varepsilon}_{xy,i}$, the window damage indicators $\hat{\mathcal{D}}_i$, and the frame damage indicator $\hat{\mathcal{D}}_f$. Damage is therefore controlled in all regions of the device, but with different mechanical relevance: frame damage is penalized more severely because it may compromise the structural integrity of the damper, whereas the window penalties are formulated to promote comparable damage levels among windows and avoid concentrating the dissipative demand in a single region. The dissipative contribution is estimated from $\hat{\varepsilon}_{xy,i}^2$, the window thickness, and the corresponding area factor (Table~\ref{tab:volume_factors}), since the energy dissipated by each window depends not only on the distortion level but also on the amount of material involved. This term is several orders of magnitude smaller than the damage penalties and is intentionally left unscaled. As a result, damage control remains the dominant criterion, while the dissipative term acts as a tie-breaker among geometries with similar damage performance, favouring those with higher distortion and, consequently, greater energy dissipation capacity.
\begin{table}[htbp]
\centering
\caption{Window area factors used to weight the distortion contribution in the objective function.}
where $W$ is the number of windows, $t_{w,i}$ is the thickness of window $i$, $A_i$ is the corresponding area factor, $\mathcal{D}_w^{\star}$ is the target damage level for the windows and $\mathcal{D}_f^{\max}$ is the maximum admissible frame damage threshold. The first term is negative because the optimizer minimizes $J$; therefore, larger energy dissipation contributions reduce the objective value and are favoured once the damage criteria are satisfied.
This formulation penalizes values above the target quadratically, while also discouraging excessively underused windows through a linear distance to the target. As a result, the optimizer tends to balance the damage levels among windows rather than forcing all windows to remain far below the admissible level.
The cubic penalty reflects the greater severity of frame damage compared with localized damage in a window. Since a damage value $\mathcal{D}=100$ represents complete failure, $\mathcal{D}_f^{\max}=90$ is adopted in the current implementation. The target window threshold is also set to $\mathcal{D}_w^{\star}=90$.
The surrogate-optimized geometry is not accepted directly. Instead, once an optimal geometry is identified by the surrogate-assisted optimizer, it is re-evaluated with FEM to verify that the surrogate remains accurate in the region of the design space where the optimum lies. This validation step checks whether the surrogate has remained reliable in the region of the design space selected by the optimizer. The candidate geometry is accepted only if: (i) the prediction error of all damage and distortion variables remains below the prescribed tolerance, equal to 5\%; (ii) the absolute error of the objective function remains within the admissible limit, set to 10; and (iii) the optimized window thicknesses remain stable between consecutive optimization iterations, with variations smaller than the 5\% of the full design range. If all criteria are satisfied, the FEM-validated geometry is accepted as the optimized design. If at least one criterion is not satisfied, the new FEM result is added to the dataset, the surrogate models are retrained and the DE optimization is repeated. This loop, that reduces the risk of accepting a geometry that is optimal only because of surrogate extrapolation error, is summarized in Figure~\ref{fig:OptimizationFlowChart}.
\caption{Surrogate-assisted optimization and FEM validation retraining loop.}
\label{fig:OptimizationFlowChart}
\end{figure}
\section{Numerical results and discussion}\label{sec:results}
The supervised-learning comparison shows a clear hierarchy among the candidate surrogate models. A total of 100 output-specific training problems were considered, corresponding to the $2W+1$ target variables required for each geometry family and adaptive iteration. The H30\_B29 family required two optimization iterations, whereas the remaining families required three. Across all geometry families, iterations and output variables, SVR was the most frequently selected model and also the model that most often achieved the lowest cross-validated RMSE. As summarized in Figure~\ref{fig:surrogate_selection_summary_barplot}, SVR was selected in 71 cases, followed by GPR, GBR, XGBoost and MLP, while Random Forest was not selected in any case. This dominance was particularly clear for the two-window families, where only one output was assigned to a GPR model and for the damage-related outputs, for which SVR was selected in 47 out of 57 cases. For more complex devices, those with three and five windows, and for distortion-related outputs, the model selection became more heterogeneous, although SVR still provided the best overall performance. From the computational point of view, SVR also offered the lowest median training times among the supervised models, whereas MLP required substantially longer training times without providing gain in accuracy.
\caption{Summary of the supervised surrogate selection over all geometry families, adaptive iterations and target outputs. The training time corresponds to the median time required for one Bayesian hyperparameter search for a single output variable and, for visualization purposes, is in logarithmic scale.}
\label{fig:surrogate_selection_summary_barplot}
\end{figure}
These results indicate that kernel-based models are particularly well suited to the present surrogate task. SVR provides the best compromise between accuracy and computational cost, while GPR is the second most competitive supervised strategy, especially in some higher-dimensional cases. Tree-based models, although robust, are less frequently selected and MLP models are not competitive in terms of computational efficiency for the dataset sizes considered here. As mentioned in previous sections, this behaviour motivated the additional evaluation of RBF interpolation as a simpler surrogate alternative. In contrast to the supervised models, RBF models were trained in less than one second per output, making them especially attractive for repeated surrogate updates within the adaptive optimization loop.
The FEM validation of the optimized geometries is summarized in Table~\ref{tab:final_surrogate_comparison}. For each geometry family and surrogate strategy, the table reports the final accepted adaptive iteration, the optimized window thicknesses $\mathbf{t}_w^{\star}$, the surrogate-predicted objective value ($J_{\mathrm{surr}}$), the corresponding FEM-recomputed objective value ($J_{\mathrm{FEM}}$) and the associated validation errors ($|e_J|$ and $e_{\max}$). The optimization process required between two and three adaptive iterations depending on the geometry family and surrogate type, with most cases converging after three iterations. No systematic difference in the number of iterations was observed between RBF and supervised ML surrogates. The maximum variable error, $e_{\max}$, is defined as the largest relative error among all quantities entering the objective function, namely ${\Exy}_i$, $\TFD_i$ and $\TFD_f$, whereas the objective-function error, $|e_J|$, is reported in absolute value.
\begin{table*}[ht!]
\centering
\caption{Final FEM validation of the optimized geometries obtained with supervised ML and RBF surrogates.}
\label{tab:final_surrogate_comparison}
\scriptsize
\begin{tabular}{llcccccc}
\toprule
Family & Surrogate & Iterations &$\mathbf{t}_w^{\star}$ [mm] &$J_{\mathrm{surr}}$&$J_{\mathrm{FEM}}$&$|e_J|$&$e_{\max}$ [\%] \\
\midrule
H30\_B29 & RBF & 3
&$[12.56,\;14.79]$
& 0.00 & 0.41 & 0.41 & 0.29 \\
H30\_B29 & Supervised ML & 2
&$[12.53,\;14.75]$
& 0.00 & 0.08 & 0.08 & 3.14 \\
\midrule
H30\_B34 & RBF & 3
&$[14.77,\;18.95]$
& 688.53 & 691.06 & 2.53 & 0.10 \\
H30\_B34 & Supervised ML & 3
&$[15.20,\;20.00]$
& 796.29 & 802.64 & 6.35 & 0.64 \\
\midrule
H45\_B29 & RBF & 3
&$[5.81,\;7.88,\;8.98]$
& 0.00 & 0.93 & 0.93 & 1.70 \\
H45\_B29 & Supervised ML & 3
&$[5.72,\;7.87,\;8.96]$
& 0.00 & 7.47 & 7.47 & 2.85 \\
\midrule
H45\_B34 & RBF & 2
&$[6.81,\;9.02,\;9.65]$
& 0.00 & 1.39 & 1.39 & 4.21 \\
H45\_B34 & Supervised ML & 3
&$[6.84,\;9.05,\;9.73]$
& 0.00 & 1.35 & 1.35 & 2.91 \\
\midrule
H60\_B34 & RBF & 3
&$[5.77,\;7.38,\;8.37,\;6.18,\;5.00]$
& 17.74 & 19.26 & 1.51 & 2.66 \\
H60\_B34 & Supervised ML & 3
&$[5.74,\;7.46,\;8.50,\;6.37,\;5.00]$
& 16.90 & 26.87 & 9.97 & 3.81 \\
\bottomrule
\end{tabular}
\end{table*}
The final validation results show that both surrogate strategies provide FEM-consistent optimized geometries after only a few adaptive iterations. RBF surrogates generally lead to lower objective-function discrepancies, with an average $|e_J|$ of 1.36, compared with 5.04 for the supervised ML surrogates. The average maximum variable error is also lower for RBF, with 1.79\% compared with 2.67\% for supervised ML. The RBF surrogate provides particularly accurate predictions for the two-window devices, with $e_{\max}$ below 0.3\%, while remaining below 4.3\% for the three- and five-window families. The supervised ML surrogates also satisfy all validation criteria, although larger discrepancies are observed in some cases, especially for the H60\_B34 family, where the objective-function error reaches 9.97.
The adaptive validation loop is essential to reach these levels of agreement. In several cases, the first surrogate-optimized candidate did not satisfy the prescribed error tolerances, particularly for the three- and five-window devices. After incorporating the additional FEM results and retraining the surrogates, the prediction errors decreased significantly. For instance, the maximum variable error of the RBF surrogate decreased from 13.70\% to 1.70\% in the H45\_B29 family and from 11.44\% to 2.66\% in the H60\_B34 family. A similar behaviour was observed for the supervised ML surrogates, whose final candidates also satisfied all acceptance criteria. This confirms that the adaptive loop reduces the risk of accepting geometries that appear optimal only because of surrogate prediction errors in sparsely sampled regions of the design space.
The optimized geometries obtained with RBF and supervised ML surrogates are very similar in most geometry families. In general, the differences in the optimized window thicknesses are below 0.1 mm and always below 0.2 mm, except for the H30\_B34 family. In this case, the differences reach 0.43 mm and 1.05 mm for the first and second windows, respectively. This family also presents the largest objective-function values, indicating that the prescribed damage targets, $\TFD_i$ and $\TFD_f$, cannot be fully maintained below the desired threshold within the admissible thickness range. Consequently, the optimizer is forced to find the best compromise between damage control and dissipative activation, which may amplify the differences between surrogate predictions.
Figure~\ref{fig:optimized_window_thickness_evolution} shows the evolution of the window thicknesses during the adaptive optimization process. In most cases, the largest adjustment occurs between the first and second iterations, while the changes between the second and third iterations are considerably smaller, indicating convergence of the proposed geometry. According to the acceptance criteria, the maximum variation in any window thickness between the last two iterations must remain below 5\% of the corresponding design range. This condition is satisfied in all validated cases, with a maximum observed variation of 4.44\%. Therefore, the adaptive process converges not only in terms of surrogate prediction accuracy, but also in terms of the optimized geometry proposed by the optimizer.
\caption{Evolution of the optimized window thicknesses during the adaptive optimization process.}
\label{fig:optimized_window_thickness_evolution}
\end{figure}
From a methodological point of view, these results highlight the trade-off between surrogate complexity, accuracy and computational efficiency. Supervised models, particularly SVR and GPR, provide high predictive accuracy and robustness, but require hyperparameter optimization and cross-validation for every output variable and adaptive iteration. RBF interpolation, in contrast, has a much lower training cost and provides very competitive final predictions once the relevant regions of the design domain have been adaptively sampled. Therefore, the comparison does not identify a universally superior surrogate strategy. Instead, it suggests that RBF interpolation is especially suitable for low- to moderate-dimensional design spaces with well-distributed FEM samples and relatively smooth input--output relationships, as occurs for the response variables analysed in this work. Supervised ML surrogates remain valuable when greater robustness is required, or when the response surface is expected to involve stronger nonlinear interactions, local irregularities or higher-dimensional dependencies.
The objective-function values should also be interpreted carefully. The objective-function error measures the consistency between surrogate predictions and FEM validation, not whether the final objective value is necessarily close to zero. Some geometry families, such as H30\_B34, retain non-negligible penalty contributions because the prescribed damage targets cannot be fully achieved within the admissible design bounds. Nevertheless, the close agreement between surrogate and FEM objective values indicates that the accepted designs are not artifacts of surrogate extrapolation, but FEM-consistent optimized candidates within the explored design space.
This behaviour is illustrated in Figure~\ref{fig:rbf_surface_evolution}, which shows the evolution of the RBF objective surface during the adaptive optimization process for the two-window families. The response surfaces remain relatively smooth, supporting the suitability of RBF interpolation for this problem. The evolution between adaptive iterations also shows how the surrogate surface is progressively corrected as new FEM information is incorporated near the optimized region.
\caption{Evolution of the RBF objective surface during the adaptive optimization process for the two-window families. Left: H30\_B29. Right: H30\_B34.}
\label{fig:rbf_surface_evolution}
\end{figure}
\section{Conclusions and future work}\label{sec:conclusions}
This work presented an adaptive surrogate-assisted optimization framework for buckling-delayed shear-link dampers subjected to cyclic seismic loading. The proposed methodology addresses the main limitation of direct FEM-based optimization, the high computational cost associated with repeatedly evaluating nonlinear cyclic simulations. By training surrogate models on FEM-generated datasets and validating the optimized candidates through additional FEM analyses, the optimizer can efficiently explore the design domain while remaining consistent with the mechanical response captured by the calibrated numerical model.
One of the key features of the proposed optimization framework, compared to approaches focused primarily on energy maximization, is the formulation of an objective function that takes damage into account. The aim is to prioritise local damage control in both the dissipative windows and the surrounding frame, whilst promoting a balanced contribution from all windows to the energy dissipation process. Thus, damage to the windows is permitted and expected, as they are intended to act as dissipative regions, provided that it remains controlled and reasonably distributed; conversely, damage to the frame is penalised more severely because it can compromise the structural integrity of the damper. Therefore, rather than merely maximizing the energy dissipated, the proposed formulation favours geometries that concentrate dissipative activation in the windows, prevent excessive damage localization in a single region and protect the frame from critical damage.
A comprehensive comparison of surrogate strategies was performed in terms of predictive accuracy and computational efficiency. The supervised ML results showed that SVR and GPR were the most competitive models, with SVR being the most frequently selected across the analysed outputs. RBF interpolation proved to be even a better alternative, due to its high efficiency: in the final adaptive iterations, it achieved validation errors comparable to, and in most cases lower than, those of the supervised ML surrogates, while requiring substantially lower training effort and no hyperparameter search. This performance is attributed to the characteristics of the present problem: low- to moderate-dimensional design spaces, well-distributed FEM samples and nonlinear but relatively smooth relationships between window thicknesses and response indicators.
The proposed adaptive validation loop proved to be necessary and effective. Several initially optimized candidates did not satisfy the prescribed error tolerances. After incorporating the new FEM results into the training dataset and retraining the surrogates, the prediction errors decreased and all final optimized geometries satisfied the acceptance criteria after only two or three iterations. Therefore, the final designs are not accepted solely on the basis of surrogate predictions, but are explicitly verified through FEM in the region of the design space where the optimum is located.
The proposed methodology also has some limitations that should be acknowledged. First, its reliability depends on the quality of the calibrated FEM model used to generate the training data and validate the optimized designs. Second, the TFDMap is used here as a post-processing damage indicator rather than as a constitutive fracture model; therefore, the optimized configurations should be interpreted in terms of relative damage control and proximity to critical states, not as direct predictions of crack initiation. Third, only the window thicknesses are considered as design variables. Although this leads to a controlled and interpretable optimization problem, it does not exploit the full geometric flexibility of BDSL dampers. Finally, the optimized geometries should ultimately be validated experimentally before being used to establish general design recommendations.
Future work should extend the design space by including additional geometric and mechanical variables, such as window height, window spacing, frame thickness or global device proportions. This extension would increase the dimensionality and complexity of the surrogate task. In those cases, the performance of RBF interpolation should therefore be reassessed. While RBF models performed very well in the present study, their efficiency and accuracy may decrease as the input space becomes larger or the response surfaces develop stronger local nonlinearities. In such cases, supervised ML models or hybrid surrogate strategies may become more advantageous.
Another relevant future direction is the incorporation of interpretability analyses, such as SHapley Additive exPlanations, to quantify the influence of each geometric variable on window damage, frame damage and dissipative activation. Although such analysis lies outside the main scope of the present study, it could provide valuable insight into the design drivers governing the behaviour of BDSL dampers and support more transparent engineering decision-making. Overall, the proposed methodology establishes a scalable basis for FEM-consistent, damage-aware optimization of seismic energy dissipation devices, while leaving room for broader design variables, richer surrogate strategies and experimental validation of the optimized configurations.
%\backmatter
\bmsection*{Author contributions}
Conceptualization: J. Irazábal, J. Ramírez, J. M. González, G. Bozzo, L. Bozzo. Methodology: J. Irazábal, J. Ramírez, J. M. González. Software and surrogate optimization: J. Irazábal. FEM modelling and validation: J. Ramírez, J. M. González, L. Lázaro, F. Rastellini. Writing--original draft: J. Irazábal. Writing--review and editing: all authors.
\bmsection*{Acknowledgments}
The authors acknowledge the financial support of Project ACE100/23/000022, ``Edificacions resilients equipades amb dissipadores Shear Link'', funded by the Government of Catalonia through ACCIO and with the support of the Catalan Office for Climate Change, with the participation of Luis Bozzo Estructuras y Proyectos S.L. and the Centre Internacional de Metodes Numerics en Enginyeria (CIMNE).
\bmsection*{Financial disclosure}
None reported.
\bmsection*{Conflict of interest}
The authors declare no potential conflict of interests.
\bibliography{../wileyNJD-AMA}
\bmsection*{Supporting information}
Additional supporting information may include the FEM database, trained surrogate models, optimization scripts and final FEM-validation simulations.
