# Neural Architecture Search with Bayesian Optimization: Designing LSTM Networks to Predict Damage Evolution in Composite Laminates

Orestis Friderikos, Arturo Mendoza, Emmanuel Baranger, John Mirisidis, Dimitrios Sagris, Constantine David, Zvonimir Tomičević

► To cite this version:

Orestis Friderikos, Arturo Mendoza, Emmanuel Baranger, John Mirisidis, Dimitrios Sagris, et al. Neural Architecture Search with Bayesian Optimization: Designing LSTM Networks to Predict Damage Evolution in Composite Laminates. 2026. hal-05496179

HAL Id: hal-05496179

https://hal.science/hal-05496179v1

Preprint submitted on 5 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Neural Architecture Search with Bayesian Optimization: Designing LSTM Networks to Predict Damage Evolution in Composite Laminates

Orestis Friderikos<fnref n="2" />, Arturo Mendoza<fnref n="3" />,<fnref n="1" />, Emmanuel Baranger<fnref n="1" />, John Mirisidis<fnref n="2" />, Dimitrios Sagris<fnref n="2" />, Constantine David<fnref n="2" />, Zvonimir Tomičević<fnref n="4" />

Received: date / Revised version: date

**Abstract** This work takes a step toward Automated Machine Learning (AutoML) by moving beyond the traditional black-box approach often seen in engineering applications of neural networks. The primary objective is to advance the understanding of Neural Architecture Search (NAS), with a specific emphasis on using Bayesian Optimization (BO) with conditional constraints to fine-tune neural architectures and hyperparameters. To the best of our knowledge, NAS frameworks remain underexplored in engineering contexts, where neural architectures are typically developed through ad hoc, manually crafted designs. In this context, Long Short-Term Memory (LSTM) recurrent neural architectures for predicting highly nonlinear delamination growth patterns in composite laminates are automatically designed using NAS. The sequential data, such as nonlinear degradation with increasing load, are generated through finite element analysis (FEA). The automated selection process produces neural network models specifically tailored to capture the complex, history-dependent nature of delamination growth. The study investigates both shallow and progressively deeper architectures through extensive systematic experiments, statistical analysis and dimensionality reduction techniques, offering insights into the inner workings of BO in optimizing neural architectures. NAS results reveal that the objective function in BO can be high irregular with no clearly defined local minima, and that deeper neural networks are not universally superior to shallower ones. Additionally, random search is used as a natural baseline to judge BO performance. The findings underscore the effectiveness of LSTM in modeling sequence-to-sequence relationships in composite damage evolution, and exhibit strong potential as surrogate models, delivering near real-time predictions, a capability challenging to achieve with traditional FEA.

## Contents

<table><tr><td>1 Introduction</td><td>2</td></tr><tr><td>2 LSTM Recurrent Neural Networks</td><td>5</td></tr><tr><td>3 Neural Architecture Search (NAS)</td><td>7</td></tr><tr><td>4 Bayesian Optimization with Gaussian Process Priors</td><td>7</td></tr><tr><td>4.1 Gaussian Process Regression</td><td>8</td></tr></table>

<footnote n="1">¹ Université Paris-Saclay, ENS Paris-Saclay, CentraleSupélec, CNRS, LMPS – Laboratoire de Mécanique Paris-Saclay, 4 avenue des sciences, 91190 Gif-sur-Yvette, France</footnote>

<footnote n="2">² Laboratory of Manufacturing Technology & Production Systems, Mechanical Engineering Department, International Hellenic University, 62124 Serres Campus, Greece</footnote>

<footnote n="3">³ Safran Tech, Rue des Jeunes Bois, 78772 Magny les Hameaux, France</footnote>

<footnote n="4">⁴ University of Zagreb, Faculty of Mechanical Engineering and Naval Architecture 10002 Zagreb, Croatia</footnote>
<table><tr><td>4.2 Kernel functions</td><td>8</td></tr><tr><td>4.3 Acquisition Functions</td><td>8</td></tr><tr><td>5 Datasets</td><td>9</td></tr><tr><td>6 Neural Architecture Search Experiments</td><td>10</td></tr><tr><td>6.1 Benchmark example</td><td>13</td></tr><tr><td>6.2 Single layer LSTM with single fully connected layer networks</td><td>14</td></tr><tr><td>6.3 Fixed neural architectures with stacked LSTM and fully connected layers</td><td>17</td></tr><tr><td>6.4 Varying deep neural architectures with stacked LSTM and fully connected layers</td><td>19</td></tr><tr><td>6.4.1 Modeling Costs</td><td>23</td></tr><tr><td>6.4.2 Principal Component Analysis of neural architectures optimized via Bayesian search: Depth and performance trends</td><td>23</td></tr><tr><td>7 Random Search Experiments</td><td>26</td></tr><tr><td>7.1 Principal Component Analysis of neural architectures optimized via random search: Depth and performance trends</td><td>28</td></tr><tr><td>8 LSTM-based metamodeling for damage evolution</td><td>30</td></tr><tr><td>9 Conclusions</td><td>34</td></tr><tr><td>10 Acknowledgments</td><td>35</td></tr><tr><td>11 Appendix</td><td>35</td></tr><tr><td>A Overview of the method</td><td>35</td></tr><tr><td>B Mechanical modeling of cohesive zone delamination</td><td>36</td></tr><tr><td>B.1 Ply drop parametric model</td><td>36</td></tr><tr><td>B.2 Ply drop FEM model</td><td>36</td></tr><tr><td>B.3 Damage accumulation due to delamination in the cohesive zones</td><td>37</td></tr></table>

# 1 Introduction

Despite the continuing growth of computing power, high non-linear FEM analysis seems to struggle to keep pace with computing advances. In the case of design optimization, the computational effort of progressively analyzing more complex models making concept exploration time-consuming and often infeasible. However, using metamodels instead of actual FEM analysis provides excellent benefits. A metamodel is an approximate model used to capture complex and non-linear input-output interactions exhibited by another more complex model. In nonlinear FEM for damage of laminate-type structures under loading excitation, a fast prediction of damage evolution patterns is of significant interest. An effective metamodel could establish an efficient and accurate mapping from the sequence of loadings to the nonlinear structural response.

The application of deep neural networks (DNNs) in the design and damage modeling of composite laminates is an emerging and rapidly advancing field. DNNs excel in tackling challenging tasks such as speech [1] and visual object recognition [2]. DNNs offer significant potential for enhancing the prediction of complex material behaviors, structural responses, and failure mechanisms in composites, which are often highly nonlinear. Unfortunately, despite their flexibility and robustness, DNNs are limited by their reliance on fixed-dimensional vector inputs and outputs. This restriction poses a significant challenge because many problems are naturally represented by sequences of variable lengths, which are not known in advance [3]. On the other hand, LSTM are a type of recurrent neural networks (RNN) specifically designed to handle sequences of data. They have an architecture that allows them to remember information over long periods and manage the flow of information through memory cells and gating mechanisms. It is noted that there exists a connection between the modeling of damage mechanics, which relies on the introduction of not observable internal variables, and the LSTM framework, in which the memory variables are implicitly related to these internal variables. LSTM can be viewed as a form of integration scheme, where multiple LSTM layers can be thought of as performing multiple levels of integration, capturing increasingly complex patterns in the time series data. Both LSTM networks for time series and numerical integration schemes share a fundamental concept of accumulating information over time. While integration schemes do this through mathematical techniques to compute integrals, LSTM achieve it through their architectural design that processes sequential data and captures temporal dependencies. Contrary to DNNs, LSTM support input data and output data
with varying sequence lengths. Therefore, LSTM represent a suitable approach for capturing temporal dependencies in sequence prediction tasks such as modeling damage progression in composite laminates, owing to their inherent capability to accommodate sequence data of variable length and dimensionality. From an engineering standpoint, introduction of LSTM for the given metamodeling task is challenging due to high accuracy requirements in industrial design.

To address the limitations of manual architecture and hyperparameter tuning in deep neural networks, early efforts introduced automated brute-force techniques such as grid search and random search. Although grid search is simple to implement, it becomes computationally infeasible in high-dimensional spaces due to the *curse of dimensionality*, where the number of combinations increases exponentially with the number of parameters and their respective ranges [4]. On the other hand, granting random search the same computational budget, it finds models that are as good or better by effectively searching a larger, less promising configuration space [5].

In recent years, Automated Machine Learning (AutoML) which is one of the fastest growing and most prominent areas in artificial intelligence today, has emerged as a transformative paradigm, aiming to reduce the manual effort and expertise required for building effective machine learning models. At the heart of AutoML lies Neural Architecture Search (NAS).

NAS has emerged as a powerful framework for automating the design of neural architectures and hyperparameter tuning [6], [7]. Unlike traditional approaches that depend heavily on expert intuition and extensive manual experimentation which are often leading to time-consuming and suboptimal designs, NAS employs sophisticated search algorithms to explore the architectural space systematically. NAS aims to advance the field of neural network design by leveraging computational techniques to discover superior architectures, thereby pushing the boundaries of what is achievable with machine learning models.

Despite its promise, the adoption of NAS is not without challenges. One of the primary concerns is the computational cost, as traditional NAS methods may require training and evaluating thousands of candidate architectures. Although recent advancements such as weight-sharing supernetworks, early stopping, and performance prediction models have substantially reduced the search cost, achieving a balance between computational efficiency and architecture quality remains an open problem. The reliance on large-scale infrastructure for even moderately complex tasks restricts accessibility for researchers and practitioners without high-performance computing resources.

An additional critical challenge lies in integrating NAS into real-world workflows where domain constraints, heterogeneous data types (e.g., static and sequential inputs), and limited data availability pose practical limitations. The application of NAS in specialized contexts such as recurrent neural networks for time-dependent systems or physics-informed modeling requires further methodological adaptation and validation. Adapting NAS to these scenarios necessitates novel search spaces, flexible encoding schemes, and constraint-aware optimization strategies, all of which increase the complexity of the pipeline.

In pursuit of more efficient and effective search strategies, various NAS methods have been proposed, including reinforcement learning-based approaches [8], [9], evolutionary algorithms [10], [11], network pruning strategies [12], convolutional neural fabrics [13], and Bayesian optimization [14]. These techniques have been successfully applied to automate neural network design, producing models that achieve state-of-the-art results in both computer vision and natural language processing tasks [6], [7], [15], [16], [17]. Their success underscores the strength of NAS in effectively navigating the vast architectural search space and uncovering novel, high-performing configurations.

Adaptive (sequential) model-based optimization methods and particularly Bayesian Optimization (BO), are perhaps more promising because they offer principled approaches to weighting the importance of each dimension. The field of NAS has seen significant advancements with the integration of BO. Below is a review of foundational and cutting-edge papers that have shaped this area.

Snoek et al. [14] introduced BO as a powerful tool for optimizing hyperparameters in machine learning models. The authors demonstrated how Gaussian Processes (GPs) and acquisition functions could be used to efficiently explore high-dimensional search spaces. White et al. [18] presents a BO framework that incorporates a neural predictor as the surrogate model. The study introduces a novel
path-based encoding scheme for neural architectures, enabling the surrogate model to effectively capture the complex relationships between different architectures and their performances. A Gaussian Process BO method for combinatorial search spaces involving ordinal and categorical variables is proposed by Oh et al. [19]. This method quantifies ‘smoothness’ of functions on combinatorial search spaces by using a combinatorial graph, where vertices represent variable assignments and edges are constructed via the graph Cartesian product. Klein et al. [20] propose a BO framework for hyperparameter tuning on large datasets. A BO procedure is constructed which models loss and training time as a function of dataset size and automatically trades off high information gain about the global optimum against computational cost. Perrone et al. [21] introduce a method to automatically design the BO search space by relying on evaluations of previous black-box functions. This approach departs from the common practice of defining a set of arbitrary search ranges a priori by considering search space geometries that are learned from historical data. A novel surrogate model for BO which combines independent GPs with a linear model that encodes a tree-based dependency structure and can transfer information between overlapping decision sequences is proposed by Jenatton et al. [22]. Additionally, a specialized two-step acquisition function is designed that explores the search space more effectively. A latent space BO NAS (LSBO-NAS) algorithm is developed by Rao et al. [23] where the neural architectures are represented as sequences. A variational auto-encoder (VAE) is trained to convert the discrete search space of NAS into a continuous latent space by learning the continuous representation of neural architectures. Nguyen et al. [24] addressed the challenge of optimizing functions defined over mixed-input spaces (categorical and continuous variables) in BO. To optimize such functions, the authors propose a method that formulates the problem as a multi-armed bandit problem, wherein each category corresponds to an arm with its reward distribution centered around the optimum of the objective function in continuous variables. Wang et al. [25] propose an ensemble BO (EBO) to operate on an ensemble of additive GP models, each of which possesses a randomized strategy to divide and conquer. Falkner et al. [26] propose a hyperparameter optimization method which combines the benefits of both BO and bandit-based methods which consistently outperforms both BO and Hyperband on a wide range of problem types.

In the realm of NAS, the development of specialized kernels to measure the similarity between neural architectures has been a subject of significant research. The selection of appropriate kernels is crucial for effectively modeling the complex, high-dimensional, and often discrete search spaces inherent to neural architectures. Various approaches have been introduced for defining and utilizing such kernels.

Kandasamy et. al [27], introduce a BO approach for NAS that employs an optimal transport distance-based kernel to measure similarities between neural architectures. The proposed method effectively captures the structural differences between architectures, facilitating a more informed and efficient search process. Ru et. al [28], introduce a BO framework for NAS that employs the Weisfeiler-Lehman (WL) graph kernel within a GP surrogate model. The approach effectively captures the topological structures of neural architectures, enabling data-efficient optimization and providing interpretability by identifying influential network motifs. The derivation of neural architectures inspired by kernels over combinatorial structures, such as sequences and graphs have been explored by Lei et al. [29]. The authors introduce a class of neural models that incorporate these kernel functions, enabling the design of architectures that are well-suited for specific structural data. Nguyen et al. [30] presents the use of optimal transport-based kernels to measure similarities between neural architectures. The proposed kernels are designed to handle both sequential and parallel structures, providing a flexible and theoretically grounded method for comparing different neural network designs.

The above studies highlight several practical guidelines for achieving effective and computationally efficient search. A well-defined and structured search space is essential to ensure that the surrogate model can generalize across candidate architectures, while the choice of surrogate should match the dimensionality and type of the variables involved. GP remain effective for continuous, low-dimensional spaces, whereas neural predictors or tree based dependency structures are better suited to discrete or higher-dimensional cases. Balanced acquisition functions such as Expected Improvement (EI) or Upper Confidence Bound (UCB) are commonly recommended to maintain an appropriate exploration-exploitation trade-off. To manage the high computational cost of NAS, low-fidelity evaluations such as
limiting the number of training epochs or employing early stopping are widely adopted and have been shown to preserve relative performance rankings under constrained budgets. Following these principles, our experimental design employs a constrained yet expressive search space, an appropriate surrogate-acquisition pair, and controlled training budgets to align with established best practices in BO-based NAS.

This research investigates the integration of NAS to systematically design and optimize neural architectures that are well-suited for capturing the complex temporal dependencies associated with damage progression, which involves complex dynamics and heterogeneous input features. The questions that inspired this paper was the following: What constitutes an effective neural architecture for accurately modeling complex engineering systems? To what extent can manually designed architectures be assumed to be optimal for such tasks? How deep/wide should a neural network be? To answer these questions, the study aims to provide deeper insight into the application of NAS frameworks for automated model design. To the best of the authors' knowledge, the use of NAS has not been systematically explored in the context of neural network-based modeling for composite laminates, thereby representing a novel contribution to the field.

By examining both the opportunities and the limitations of NAS, this work contributes to the advancement of more efficient, interpretable, and scalable machine learning pipelines for engineering applications. The investigation is grounded in a rigorous experimental framework, employing BO as the core search strategy. The NAS process is further enhanced through comprehensive statistical analyses aimed at elucidating the performance characteristics and decision patterns of BO-driven architectural exploration. Through this approach, the research seeks to advance the understanding of how automated neural architecture design can be effectively applied to solve challenging engineering problems.

The remainder of the paper is organized as follows:

Section 2 offers a general introduction to LSTM networks. NAS and BO with GP priors are introduced in Section 3 and Section 4, respectively. Information regarding the input and FEM-generated output dataset used to train the neural network is provided in Section 5. The NAS experiments using BO are elucidated in Section 6. Benchmark NAS experiments using random search are illustrated in Section 7. The LSTM-based metamodel for damage evolution for various combinations of input variables is shown in Section 8. Lastly, Section 9 concludes by summarizing the main outcomes, and highlighting some important perspectives.

## 2 LSTM Recurrent Neural Networks

Recurrent neural networks (RNN) with Long Short-Term Memory have emerged as effective and scalable models for sequence-to-sequence modeling tasks [31], [32]. Given a sequence of inputs $x_t$, LSTM computes a sequence of outputs $y_t$ by iterating the *hidden state* $h_t$, where the subscript $t$ represents time. The core concept behind LSTM is the structure of *cells* and the regulation of information flow within these cells. At time step $t$, the LSTM *block* uses the current state of the network ($c_{t-1}$, $h_{t-1}$) and the next time step of the sequence $x_t$ to compute the output and update the *cell state* $c_t$.

The LSTM cell consists of a connection mechanism between four interacting units, including three controlling *gates*, e.g., an input gate $i_t$, a forget gate $f_t$, an output gate $o_t$, and an internal cell state $c_t$ (Fig. 1). The relationship among these defined variables can be described by the equations

$$i_t = \sigma(W^{xi}x_t + R^{hi}h_{t-1} + b_i) \qquad (1)$$

$$f_t = \sigma(W^{xf}x_t + R^{hf}h_{t-1} + b_f) \qquad (2)$$

$$o_t = \sigma(W^{xo}x_t + R^{ho}h_{t-1} + b_o) \qquad (3)$$

$$\bar{c}_t = \tanh(W^{xc}x_t + R^{hc}h_{t-1} + b_c), \quad (\text{intermediate cell candidate}) \qquad (4)$$

where the matrices $W^{x*}$, $R^{h*}$ are the input weights and recurrent weights, respectively, while $b_c$ represents the corresponding bias vector within different gates $* \in \{i, f, o, c\}$. These learnable parameters
are randomly initialized and will be iteratively updated during neural network training. Function $\sigma$ denotes the logistic sigmoid function, which outputs values between 0 and 1.

The cell state (long-term memory) at time step $t$ is given by

$$c_t = f_t \odot c_{t-1} + i_t \odot \bar{c}_t \quad (5)$$

The hidden state (short-term memory) at time step $t$ is given by

$$h_t = o_t \odot \tanh(c_t) \quad (6)$$

where the symbol $\odot$ denotes the Hadamard product (element-wise product of vectors).

The input gate $i_t$ controls the updating of the cell state by deciding which new information from the current input ($x_t$) and the previous hidden state ($h_{t-1}$) should be added. The forget gate $f_t$ decides what information from the previous cell state ($c_{t-1}$) should be discarded. The output gate $o_t$ determines what part of the cell state should be outputted as the hidden state. This hidden state will then be used for the next time step and as part of the final output. The cell state update $\bar{c}_t$ generates new candidate values that could be added to the cell state, based on the current input and the previous hidden state. The cell state $c_t$ acts as a memory line that runs through the entire sequence. The hidden state $h_t$ is the short-term memory of the LSTM. By combining these mechanisms, LSTMs can effectively remember and utilize long-term dependencies in sequential data. Summarizing, a conventional LSTM is constructed by defining a state transition function $f_h$ and an output function $f_o$ as

$$h_t = f_h(x_t, h_{t-1}) \quad (7)$$

$$y_t = f_o(h_t, x_t) \quad (8)$$

The complex state transition function $f_h$ renders the LSTM network powerful in sequence modeling, while the output function $f_o$, such us a fully connected layer (multilayer perceptron) proves beneficial for mapping the temporal features to the corresponding output space.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Architecture of a typical LSTM cell which consists of cell input $x_t$, cell output $y_t$, cell state $c_t$, hidden state $h_t$, and four gate variables $i_t$, $f_t$, $o_t$, $\bar{c}_t$. The bias terms (which are added to four gate variables) have been omitted for clarity.</figcaption></figure>
## 3 Neural Architecture Search (NAS)

NAS focuses on automating the design of neural architectures and hyperparameter tuning. It aims to advance the field of neural network design by leveraging computational techniques to discover superior architectures, thereby pushing the boundaries of what is achievable with machine learning models. NAS methods can be categorized according to three dimensions: *search space*, *search strategy*, and *performance estimation strategy* [7] (Fig. 2). The *search space* $\mathcal{X}$, i.e., the set of feasible solutions, is a key component of NAS. It influences both the search duration and the quality of the solution. It is thus important to constrain the search space to achieve a good compromise between computation time and effectiveness. The *search strategy* details how to explore the search space (which is often exponentially large or even unbounded). It encompasses the classical exploration-exploitation trade-off since, on the one hand, it is desirable to find well-performing architectures quickly, while on the other hand, premature convergence to a region of suboptimal architectures should be avoided. *Performance estimation strategy* defines how to estimate/predict the performance of a given neural architecture in the design space. The simplest option is to perform a standard training and validation of the architecture on data, but this is unfortunately computationally expensive and limits the number of architectures that can be explored. Much recent research therefore focuses on developing methods that reduce the cost of these performance estimations.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Abstract illustration of Neural Architecture Search methods.</figcaption></figure>

## 4 Bayesian Optimization with Gaussian Process Priors

BO emerged as an intelligent search strategy to explore the search space $\mathcal{X}$, serving as an alternative to manual, grid, or random search. In BO for NAS, the goal is to find the neural network topology and/or set of hyperparameters in $\mathcal{X}$ that lead to the best validation accuracy. Formally, BO is a class of machine-learning-based optimization methods focused on solving the problem

$$ \mathbf{x}^* = \underset{\mathbf{x} \in \mathcal{X}}{\operatorname{argmax}} f(\mathbf{x}), \quad f: \mathcal{X} \subseteq \mathbb{R}^D \to \mathbb{R} \qquad (9) $$

where the index set $\mathcal{X}$ is the set of possible inputs, which could be $\mathbb{R}^D$. BO is a sequential search algorithm well-suited for optimizing hyperparameters of classification and regression models [33], [14]. It is used to optimize expensive black-box functions that are nondifferentiable, discrete or continuous, and time-consuming to evaluate. BO has recently became the standard strategy for Automated Machine Learning (AutoML) and NAS [34]. Moreover, experimental results show that BO outperforms other global optimization algorithms [35].

BO proceeds by constructing a *probabilistic surrogate model* that combines all available information from previous evaluations and does not rely solely on local gradient and Hessian approximations. Then, it employs an *acquisition function* that combines information about the expectation and uncertainty of a function value to make decisions where to sample next (Algorithm 1). BO is able to achieve a balance between exploitation and exploration [36]. When using BO, the surrogate model is often chosen to be a Gaussian Process (GP) prior as a probability model due to its flexibility, ability to provide uncertainty
estimates, and solid theoretical foundations [37]. Gaussian Processes are particularly well-suited for problems where the objective function is smooth and relatively low-dimensional. However, the choice of the surrogate model in BO is arbitrary and can be tailored to fit the specific characteristics of the problem at hand. While Gaussian Processes are a common choice, other models such as Random Forests, Bayesian Neural Networks, or other machine learning models can also serve as effective surrogates in BO.

## 4.1 Gaussian Process Regression

A GP is a collection of random variables, any finite number of which have a joint Gaussian distribution [37]. If a real process {$f(\mathbf{x}), \mathbf{x} \in \mathbb{R}^D$} is a GP, then given any finite collection of $n$ observations {$\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_n$}, the joint distribution of the random variables $f(\mathbf{x}_1), f(\mathbf{x}_2), \dots, f(\mathbf{x}_n)$ is assumed to be jointly Gaussian. GP is completely specified by a mean function $m(\mathbf{x})$ and a covariance (kernel) function $k(\mathbf{x}, \mathbf{x}')$ of a real process $f(\mathbf{x})$ as

$$m(\mathbf{x}) = \mathbb{E}[f(\mathbf{x})], \quad \mathcal{X} \to \mathbb{R} \qquad (10)$$

$$k(\mathbf{x}, \mathbf{x}') = \mathbb{E}[(f(\mathbf{x}) - m(\mathbf{x}))(f(\mathbf{x}') - m(\mathbf{x}'))], \quad \mathcal{X} \times \mathcal{X} \to \mathbb{R} \qquad (11)$$

Thus we can write the GP as

$$f(\mathbf{x}) \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}')) \qquad (12)$$

## 4.2 Kernel functions

GP uses a positive-definite kernel function $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$. The kernel function expresses the covariance between two latent variables $f(\mathbf{x}_i)$ and $f(\mathbf{x}_j)$, where $\mathbf{x}_i, \mathbf{x}_j \in \mathbb{R}^D$. In other words, it can be considered as a similarity metric which determines how the response at one point $\mathbf{x}_i$ is affected by the response at other point $\mathbf{x}_j$. The covariance function $k(\mathbf{x}_i, \mathbf{x}_j)$ can be defined by various kernel functions a priori, and the efficiency of BO will rely on this expert choice. Matlab function bayesopt uses the popular automatic relevance determination (ARD) Matérn 5/2 kernel [38] with a separate length scale $\sigma_m$ for each predictor $m = 1, 2, \dots, D$. The characteristic length scales briefly define how far apart the input values $\mathbf{x}_i$ can be for the response values to become uncorrelated, e.g., the distance before the response can change significantly. Additionally, the kernel can be parameterized in terms of the GP kernel hyperparameters in vector $\theta$ and it is possible to express the covariance function as

$$k(\mathbf{x}_i, \mathbf{x}_j | \theta) = \sigma_f^2 \left(1 + \sqrt{5}r + \frac{5}{3}r^2\right) \exp\left(-\sqrt{5}r\right) \qquad (13)$$

where $\sigma_f$ is the signal standard deviation and

$$r = \sqrt{\sum_{m=1}^{d} \frac{(x_{im} - x_{jm})^2}{\sigma_m^2}} \qquad (14)$$

## 4.3 Acquisition Functions

For a given GP approximation of $f(\mathbf{x})$, the acquisition function $\alpha(\mathbf{x}) : \mathcal{X} \to \mathbb{R}$, determines which point in $\mathcal{X}$ should be evaluated next. In this work the Expected Improvement (EI) [39] is chosen

$$\alpha(\mathbf{x}) = \mathbb{E}[\max(f(\mathbf{x}^*) - f(\mathbf{x}), 0)] \qquad (15)$$
where $\mathbf{x}^* = \operatorname*{argmin}_{\mathbf{x}} f(\mathbf{x})$ is the location of the lowest observed posterior mean and $f^*(\mathbf{x}) = f(\mathbf{x}^*)$ is the lowest value of the posterior mean. Other common acquisition functions include lower or upper confidence bounds (UCB) and probability of improvement. EI is better at exploring the search space and finding global minima, while others, such as UCB, are better at exploiting the current best point and finding local minima. Ultimately, the goal of BO is to identify effective hyperparameter settings as quickly as possible. While acquisition procedures such as EI optimize progress in the next function evaluation, in practice, wall-clock time is the more relevant measure. To improve performance in terms of wall-clock time, 'expected-improvement-per-second-plus' is utilized in bayesopt function. Thus, bayesopt can obtain better improvement by using time-weighting in its acquisition function which prefers to acquire points that are not only likely to be good, but that are also likely to be evaluated quickly.

<caption>**Algorithm 1:** Bayesian Optimization for Neural Architecture Search (BO-NAS)</caption>

**Input:** Architecture search space $\mathcal{X}$; expensive evaluator $\mathcal{E}(\mathbf{x}, \zeta)$ returning performance $f(\mathbf{x})$ at fidelity $\zeta$ (e.g., epochs); surrogate prior (GP) with hyperparameters $\theta$; acquisition function $\alpha(\cdot | \mathcal{D})$; total BO iterations $N$.
**Output:** Best architecture $\mathbf{x}^*$ and its observed performance $f^*(\mathbf{x})$.
1 Set dataset $\mathcal{D}_0 \leftarrow \emptyset$;
2 Specify surrogate prior $f(\mathbf{x}) \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}')); \theta$;
3 Set $n \leftarrow 1$;
4 **while** $n \le N$ **do**
5 Compute acquisition function $\alpha(\mathbf{x} | \mathcal{D}_{n-1})$ using the current surrogate posterior. Select next architecture $\mathbf{x}_n$
$$ \mathbf{x}_n = \arg\max_{\mathbf{x} \in \mathcal{X}} \alpha(\mathbf{x} | \mathcal{D}_{n-1}, \theta) $$
6 Evaluate performance $f_n(\mathbf{x}) \leftarrow \mathcal{E}(\mathbf{x}_n, \zeta)$;
7 Augment data: $\mathcal{D}_n \leftarrow \mathcal{D}_{n-1} \cup \{(\mathbf{x}_n, f_n(\mathbf{x}))\}$;
8 Update surrogate posterior $f_n(\mathbf{x}) | \mathcal{D}_{n-1}, \mathbf{x}_n \sim \mathcal{GP}(\mu(\mathbf{x}_n | \mathcal{D}_{n-1}), \sigma(\mathbf{x}_n | \mathcal{D}_{n-1}); \theta)$;
9 $n \leftarrow n + 1$;
10 Let $(\mathbf{x}^*, f^*(\mathbf{x})) = \arg\max_{(\mathbf{x}, f(\mathbf{x})) \in \mathcal{D}_N} f(\mathbf{x})$; // best observed architecture
11 **return** $\mathbf{x}^*, f^*(\mathbf{x})$

## 5 Datasets

This study builds upon the database generated in previous works [40,41] and specifically focuses on the automated design of LSTM architectures using NAS. A concise overview of the overall framework is presented in Appendix A. The parametric ply-drop laminate FEM model (see Appendix B) incorporates *static* variables related to geometry ($t_1, t_2, t_3$) and material orientations ($\theta_1, \theta_2, \theta_3$). The *dynamic* variables $u_1, u_2, u_3$ are the imposed displacements and rotation to the laminate structure [40,41]. On the other hand, the output space $\mathcal{Y}$ contains the response vectors $\mathbf{y}^{(s)} \in \mathbb{R}, s = 1, \dots, S$. Different nonlinear output responses of the structure degradation (spaces $\mathcal{Y}$) have been considered in previous works [40, 41]. The method proposed herein addresses the evolution of stiffness loss associated with rotational loading during the LPA steps. The LSTM maps the input vector sequence $\mathbf{x}^{(s)}, s = 1, \dots, S$ (ensemble) to the stiffness loss evolution sequence $\mathbf{y}^{(s)} \in \mathcal{Y}$ (sequence-to-sequence input-output regression).

The effective set $D_{all}$ is composed of 49588 sequences (FEM simulations). This set is initially split into a training-validation set $D_{train / val}$ and a test set $D_{test}$ using randomization with a 75-25% split. Then, the training set $D_{train}$ and the validation set $D_{val}$ are defined through randomization with a 90-10% split of the $D_{train / val}$ set, respectively. Consequently, $N_{train} = 33471, N_{val} = 3720$ and $N_{test} = 12397$ represent the number of sequences of each set, respectively. Furthermore, the design variables are normalized to be bounded within the range $[-1, +1]$ using min-max normalization scaling across each variable.

**Remark 1:** The LSTM developed does not incorporate any physical laws or constraints into its design or training process, e.g., zero initial damage or positive increments of damage. It purely relies
on data-driven learning. This is in opposition to the current trend in machine learning where models are designed or constrained using physical principles or domain knowledge [42]. This includes ensuring properties like convexity, conservation laws, or other physical constraints are respected by the model.

**Remark 2:** The LSTM utilizes the functional approach where the history of the material's behavior is directly integrated into a functional form which leads to the prediction of the system's response without relying on internal variables. This method is commonly used in viscoelastic materials, where the history of stress and strain influences current behavior [43]. The objective is to define the best functional representation of a model that traditionally uses internal variables. By focusing on the functional representation of the history within the LSTM, the research aims to leverage the strengths of LSTMs in handling sequential data to capture and predict complex material behaviors.

# 6 Neural Architecture Search Experiments

**Search Space.** There is no existing prior knowledge about the neural architectures well-suited to reduce the size of the search space and simplify the search for the problem under consideration. To address this, a common practice of defining a set of arbitrary search ranges a priori is employed. In this work, the search space $\mathcal{X}$ is a hyper-rectangle $\{\mathbf{x} \in \mathbb{R}^D, \alpha^i \le x^i \le \beta^i\}$.

*Network Topology.* To handle both *static* and *dynamic* inputs, the *static* (non-time series) variables are repeated across time steps and concatenated with the *dynamic* variables at each time step (Fig. 3). Hyperparameters for the neural architectures are the number of stacked LSTM layers $\alpha^{LSTM}$, the width of each LSTM layer $A_i^{LSTM}$ (number of units in hidden layer $i$), the number of additional fully connected dense layers $\alpha^{FC}$ for consolidation after the LSTM layers, and their corresponding number of units $A_j^{FC}$ (in layer $j$).

*Hyperparameters.* Training hyperparameters govern the mini batch size $\xi$, learning rate $\gamma$, and $L_2$-regularization [44] for deeper models.

A network-level search space with fixed connection patterns is used, e.g., the entire architecture is explored and optimized, rather than focusing on fine-grained components like individual layers or operations. Layer width scaling is not realized for the neural network topology. Moreover, only scalar hyperparameters are searched. All categorical hyperparameters remain fixed for all experiments (Table 1).

Different network topologies might require different settings of hyperparameters to achieve optimal performance. Recently, it has been reported that using NAS in a two step process, e.g., tuning separately neural network topology and hyperparameters in a post-hoc step can lead to suboptimal solutions [45]. Likewise, it has been demonstrated that the common practice of using very few epochs during the main NAS and much larger numbers of epochs during a post-processing step is inefficient due to little correlation in the relative rankings for these two training regimes. On the other hand, joint NAS is better suited for exploration and discovering novel architectures and can lead to models that generalize better across a wider range of hyperparameter choices. For these reasons, an end-to-end approach is chosen, i.e., jointly tuning neural architectures and hyperparameters, even though it can be more complex to manage and computationally expensive as it involves a larger search space $\mathcal{X}$ (Table 1).

In this work, neural network depth is gradually increased by introducing additional layers, each with an arbitrary number of units (i.e., arbitrary width). Following Pascanu et al. [46], two complementary strategies are adopted to extend a RNN into a deep RNN: (i) deep transition functions, realized by stacking multiple LSTM hidden layers on top of each other to increase the expressive power of the recurrent dynamics [47], and (ii) deep output functions, realized by incorporating multiple fully connected layers between the final hidden state and the output layer, thereby enriching the hidden-to-output mapping. These two approaches are illustrated in Fig. 3.

The search space $\mathcal{X}$ encompasses two main configurations: (i) fixing the neural network depth while exploring optimization hyperparameters, and (ii) jointly searching over neural network depth and optimization hyperparameters. The design space explored in the NAS experiments is summarized in Table 2. Architectures were grouped into four batches: B.1, B.2, B.3, and B.4. The neural network depth increases progressively across these batches, from the shallow models in B.1 to deeper models in B.4. In
<table><caption>Table 1: NAS configuration space for architecture and hyperparameter optimization (black color indicates fixed parameters whereas blue color indicates parameters used in BO)</caption><thead><tr><th>Architecture/Hyperparameters</th><th>Type/Range of param./values</th></tr></thead><tbody><tr><td>Stacked LSTM layers (h<sub>t</sub> → ... → z<sub>t</sub>)</td><td>α<sup>LSTM</sup> ∈ [1, 3], [1, 4]</td></tr><tr><td>Width of LSTM layers</td><td>A<sub>i</sub><sup>LSTM</sup> ∈ [20, 240]</td></tr><tr><td>Intermediate fully connected layers (hidden-to-output h<sub>t</sub> → y<sub>t</sub>)</td><td>α<sup>FC</sup> ∈ [1, 2], [1, 6]</td></tr><tr><td>Number of units at fully connected layers</td><td>A<sub>j</sub><sup>FC</sup> ∈ [20, 150]</td></tr><tr><td>Learning rate</td><td>γ ∈ [9e<sup>-5</sup>, 5e<sup>-3</sup>]</td></tr><tr><td>Minibatch</td><td>ξ ∈ [32, 256]</td></tr><tr><td>L<sub>2</sub>-Regularization</td><td>λ ∈ [1e<sup>-9</sup>, 1e<sup>-3</sup>]</td></tr><tr><td>Epochs</td><td>{250, 500, 750, 1200}</td></tr><tr><td>Activation function</td><td>LeakyReLU</td></tr><tr><td>Solver</td><td>Adam optimizer</td></tr><tr><td>Shuffle data</td><td>Every epoch</td></tr><tr><td>Validation frequency</td><td>Every 3 epochs</td></tr><tr><td>Learning rate scheduler</td><td>Piecewise</td></tr><tr><td>Regression layer</td><td>Loss function: RMSE</td></tr></tbody></table>

batch B.1, simple baseline architectures consisting of a single LSTM layer followed by a single fully connected layer were investigated. Batches B.2 and B.3 introduced moderately deeper architectures with fixed depths, allowing the impact of increased complexity to be assessed systematically. In contrast to the previous cases, in B.4, both the number of LSTM layers and the number of intermediate fully connected (dense) layers were treated as tunable hyperparameters in BO. The last column of Table 2 lists the specific hyperparameters optimized via BO for each batch.

<table><caption>Table 2: Depth of neural networks in NAS experiments (black color indicates fixed parameters whereas blue color indicates parameters used in BO)</caption><thead><tr><th>Batch of experiments<br>B.#</th><th>Stacked LSTM layers</th><th>Fully connected layers</th><th>BO trials<br>T.#</th><th>Epochs</th><th>BO hyper-parameters</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1</td><td>5</td><td>{250, 500, 750}</td><td>4</td></tr><tr><td>2</td><td>2</td><td>2</td><td>5</td><td>{250, 500, 750}</td><td>7</td></tr><tr><td>3</td><td>3</td><td>2</td><td>5</td><td>{250, 500, 750}</td><td>8</td></tr><tr><td>4</td><td>[1, 4]</td><td>[1, 6]</td><td>1</td><td>1200</td><td>15</td></tr></tbody></table>

**Remark 3:** We refer to the conceptual and functional similarities between LSTM and numerical integration schemes used in modeling of physical systems. Since LSTM are particularly good at capturing dependencies over long sequences, this ability is somewhat analogous to integration schemes in numerical analysis, which accumulate values over time to provide a continuous output. When multiple LSTM layers are stacked, each layer processes the sequence data, adding another level of integration or accumulation of information. This ‘integration of integration’ means that the LSTM layers can capture increasingly complex and hierarchical dependencies in the data. By using multiple layers, the network effectively builds a sophisticated functional representation of the system’s history. Thus, LSTM naturally align with the functional approach by continuously refining the representation of the system’s behavior as it processes sequential data.

**Remark 4:** Before proposing a deeper LSTM, it is essential to obtain quantitative evidence that this is required, and not a much simpler model would be better suited for the specific task at hand. It has been theoretically shown that a single layer LSTM with saturated linear activation is Turing-
complete [48], which means that the network has the capability to perform any computation that a Turing machine can, given enough resources such as time and memory.

**Search Strategy.** BO is used as a search strategy over a range of potential network topologies and hyperparameters ([33], [14]). The most important BO variables are summarized in Table 3. BO is performed through the Matlab function *baysopt*, which internally uses the Gaussian Process Regression (GPR) function *fitrgp* to model the data.

**Performance Estimation Strategy.** Each model is trained from scratch which inevitably incurs a significant computational cost. The Bayesian optimizer returns the RMSE on the validation set as the objective function to be minimized.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Architecture design in the NAS design space.</figcaption></figure>

<table><caption>Table 3: Bayesian optimization variables</caption><tr><td>Optimizer type (surrogate)</td><td>Gaussian Process Regression</td></tr><tr><td>Kernel function</td><td>Automatic relevance determination (ARD) Matérn 5/2 kernel</td></tr><tr><td>Acquisition function</td><td>Expected-improvement-per-second-plus</td></tr><tr><td>Max objective function evaluations</td><td>M = 120</td></tr></table>

For ease of neural network comparison, the training $D_{train}$, validation $D_{val}$ and test $D_{test}$ sets are the same for all NAS experiments. The Adam optimizer is used with a mini-batch size $\xi \in [32, 256]$. Learning rate is searched in the range $\gamma \in [9e^{-5}, 5e^{-3}]$ using a logarithmic transform. The training $D_{train}$ set is shuffled before each epoch to maximize feature learning, e.g., reduce sequence bias, prevent overfitting, improve generalization, enhance optimization by avoiding local minima, and ensure that the training is less sensitive to the order of the data. The rectified linear unit (ReLU) is used as an activation function. In approximation theory, it is shown that neural networks with such unbounded activation functions have the universal approximation property [49].
It has to be remarked that re-running BO to the same dataset does not necessarily lead to the same result as BO is stochastic method. Thus, optimization is repeated several times (see trials used in each experiment in Table 2). NAS experiments are evaluated using a variable number of epochs with values {250, 500, 750, 1200}. Each trial was given a fixed computational budget of M = 120 function evaluations.

The network training has been performed in Matlab using Deep learning Toolbox, which provides a well-documented framework for machine learning applications. It offers flexible dataflow architecture enabling high-performance training of various types of neural networks on a variety of platforms.

## 6.1 Benchmark example

Before proceeding to the NAS experiments described in the previous section, an illustrative example is presented using a low dimensional design space consisting only two hyperparameters. A fixed topology is used for the search space in terms of layers. Specifically, an architecture with a single LSTM layer and a single fully connected layer is considered, with two hyperparameters selected for BO: the number of LSTM units and the learning rate. Low-fidelity estimates are employed for the performance estimation strategy, where early stopping (e.g., 50 epochs) is combined with BO (e.g., 60 objective function evaluations). This approach reduces computational cost while still providing reasonable approximations of model performance during the optimization process. For the size of the recurrent layer, the search is conducted over $A_1^{LSTM} \in [20, 700]$ units. Additionally, the fully connected layer is fixed at $A_1^{FC} = 60$ units, followed by a Leaky ReLU activation function. The learning rate is searched within the range $\gamma \in [5e^{-5}, 5e^{-2}]$ using a logarithmic transform. The batch size is fixed at $\xi = 50$.

The performance of BO versus the number of function evaluations is shown in Fig. 4. This convergence plot illustrates the improvement during optimization, providing insight into when the optimizer is still making progress and when further improvements are minimal or nonexistent. The vertical axis represents the prediction accuracy (RMSE) of (i) the estimated minimum objective value (posterior mean) and (ii) the minimum observed objective value on the validation set. The horizontal axis represents the number of function evaluations performed by the BO algorithm. From Fig. 4, it can be seen that the smallest value is observed at evaluation point 20, after which no further improvement is observed.

Additionally, Fig. 5 and Fig. 6 show the objective function model after 30 and 60 function evaluations, respectively. The surface represents the posterior mean of the GP model for the objective function over the 2D hyperparameter space. The scattered blue dots represent the points where the objective function has been evaluated. In the first 30 function evaluations, the optimizer tends to exploit regions close to the observed minimum rather than exploring the design space broadly (evidenced by the clustering of points near the minimum observed value and the sparse distribution of points in other regions). After 60 function evaluations, the objective function becomes nearly flat, indicating that the neural network performance landscape is highly irregular, with no clearly defined local minima. One can observe, that the estimated minimum objective value after 27 evaluations has reached a plateau (Fig. 4). This behavior is a common challenge in NAS, where the objective function landscape is often complex and lacks clear structure.

The sequential evaluation points of the BO, considering the sets $D_{train}$, $D_{val}$, $D_{test}$ is shown in Fig. 7. The horizontal axis represents the iterations of BO algorithm, and the vertical axis represents the prediction accuracy in terms of RMSE for the train, validation and test sets at each function evaluation. The curve corresponding to the validation set $D_{val}$ shows the history of objective function values (RMSE) over all iterations, where each point on the plot represents the objective function value at a specific iteration. The fluctuations observed in this plot is a natural outcome in BO which indicates that the method explores different regions in the design space.
<figure><img src="image_5.png" /><figcaption>Fig. 4: Minimum objective value versus functions evaluations. The horizontal axis represents the it-
erations of BO algorithm, while the vertical axis shows the prediction accuracy on the validation set
(RMSE).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5: Objective function model after 30 function evaluations.</figcaption></figure>

6.2 Single layer LSTM with single fully connected layer networks

This model represents the simplest baseline neural architecture, and the aim is to provide details on
the training hyperparameters (in terms of units, learning rate, etc) and performance. In this batch of
experiments (B.1), BO is used to jointly tune 4 hyperparameters of each neural network (Table 2).
Table 1 summarizes the hyperparameter selection, and the corresponding training options with the
relative bounds (the blue color indicates parameters for BO). Given the exponential growth of the search
space in terms of the number and range of hyperparameters, the following constraints are imposed:
<figure><img src="image_7.png" /><figcaption>Fig. 6: Objective function model after 60 function evaluations.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7: The sequential evaluation points of BO considering the *D*<sub>train</sub>, *D*<sub>val</sub>, *D*<sub>test</sub> versus the number of function evaluations; the horizontal axis represents the iterations of BO algorithm; the vertical axis represents the prediction accuracy (RMSE) on the train (red), validation (green) and test (blue) set for each function evaluation.</figcaption></figure>

the model architectures consist of a single LSTM layer $\alpha^{LSTM} = 1$ with a single fully connected layer $\alpha^{FC} = 1$. For the size of recurrent layer, search is conducted over $A_1^{LSTM} \in [20, 240]$ units. Additionally, $A_1^{FC} \in [20, 150]$ units for the fully connected network are allowed, followed by a leaky ReLU activation function. We experiment by adjusting different $L_2$-regularization strengths, represented by different values of $\lambda$, to identify the optimal setting for the specific task. However, it is found that $L_2$-regularization did not align well with the model structure, as it degenerates performance by causing underfitting. This was evidenced by persistently high training and validation losses. Based on these experiments, it was concluded that the model does not benefit from $L_2$-regularization, and therefore, it was not used.

After BO, the neural architecture that performs best on the validation set (minimum RMSE) is used for inference on new, unseen data. Hence, the five BO trials are ranked, and from the top two
BO trials, the top five models are selected to evaluate their performance on the test set $D_{test}$. The ranking of the top five feasible points observed during the first and second best-performing BO trial T.1 (750 epochs) and T.2 (500 epochs) in batch B.1, along with their corresponding hyperparameters, is presented in Table 4 and Table 5, respectively. As shown in Table 4, the model with the minimum observed objective function value (RMSE = 0.0068556) resulted in the largest prediction error on the test set $D_{test}$ (RSME =0.0087131) among the top five feasible points.

<table><caption>Table 4: Ranking of the top five feasible points observed during the best-performing BO trial T.1 in experiment batch B.1.</caption><thead><tr><th>Architecture A.#</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>N. of units at layer A<sub>1</sub><sup>LSTM</sup></td><td>194</td><td>222</td><td>81</td><td>239</td><td>78</td></tr><tr><td>N. of units at layer A<sub>1</sub><sup>FC</sup></td><td>32</td><td>97</td><td>103</td><td>83</td><td>50</td></tr><tr><td>Learning rate γ</td><td>9.2279e-04</td><td>1.9205e-04</td><td>0.0019</td><td>2.7112e-04</td><td>7.0418e-04</td></tr><tr><td>Minibatch ξ</td><td>102</td><td>137</td><td>106</td><td>96</td><td>103</td></tr><tr><td>D<sub>val</sub> (min objective RMSE)</td><td><strong>0.0068556</strong></td><td>0.0069445</td><td>0.0071253</td><td>0.0072028</td><td>0.0072117</td></tr><tr><td>D<sub>train</sub> (RMSE)</td><td>0.0036077</td><td>0.0047615</td><td>0.0035527</td><td><strong>0.0031552</strong></td><td>0.0046375</td></tr><tr><td>D<sub>test</sub> (RMSE)</td><td>0.0087131</td><td>0.0082648</td><td><strong>0.0076794</strong></td><td>0.0081216</td><td>0.0079892</td></tr></tbody></table>

<table><caption>Table 5: Ranking of the top five feasible points observed in the second-best BO trial T.2 in experiment batch B.1</caption><thead><tr><th>Architecture #</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>N. of units at layer A<sub>1</sub><sup>LSTM</sup></td><td>185</td><td>125</td><td>183</td><td>160</td><td>118</td></tr><tr><td>N. of units at layer A<sub>1</sub><sup>FC</sup></td><td>20</td><td>88</td><td>26</td><td>56</td><td>144</td></tr><tr><td>Learning rate γ</td><td>5.7286e-04</td><td>0.0049</td><td>0.0025</td><td>3.9587e-04</td><td>3.6026e-04</td></tr><tr><td>Minibatch ξ</td><td>189</td><td>204</td><td>190</td><td>102</td><td>73</td></tr><tr><td>D<sub>val</sub> (min objective RMSE)</td><td><strong>0.0070433</strong></td><td>0.0071908</td><td>0.0072919</td><td>0.0074762</td><td>0.0074916</td></tr><tr><td>D<sub>train</sub> (RMSE)</td><td>0.0041591</td><td><strong>0.0029043</strong></td><td>0.0042039</td><td>0.0044284</td><td>0.0051709</td></tr><tr><td>D<sub>test</sub> (RMSE)</td><td>0.0079463</td><td><strong>0.0076855</strong></td><td>0.0083283</td><td>0.0078231</td><td>0.0083223</td></tr></tbody></table>

The performance of BO versus the number of function evaluations for the BO trial T.2 is displayed in Fig. 8. Although there is a discrepancy between the estimated minimum objective value and the minimum observed objective value, the estimated minimum objective value exhibits a decreasing trend that aligns with the decrease in the minimum observed objective value until approximately 60 function evaluations. Beyond this point, the GP surrogate model reaches a plateau. This plateau indicates the highly irregular landscape of the neural network models, where further improvements become negligible.

The sequential evaluation points of the BO, considering the sets $D_{train}$, $D_{val}$, $D_{test}$ versus the number of function evaluations for the second best-performing trial T.2 is shown in Fig. 9. One can observe the convergence of the curves associated to the validation $D_{val}$ and test $D_{test}$ sets which confirms i) the generalization ability of the validation set, e.g., it is representative of the true distribution of the data, and ii) the BO performance is not biased by predictions on the validation set. Note that the fluctuations observed in these plot is a natural outcome in BO which indicates that the method explores different regions in the design space.
<figure><img src="image_9.png" /><figcaption>Fig. 8: Minimum objective value versus functions evaluations for the second best-performing BO trial T.2. The horizontal axis represents the iterations of BO algorithm, while the vertical axis shows the prediction accuracy on the validation set (RMSE).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: The sequential evaluation points of BO considering the $D_{train}$, $D_{val}$, $D_{test}$ versus the number of function evaluations for the second best-performing trial T.2; the horizontal axis represents the iterations of BO algorithm; the vertical axis represents the prediction accuracy (RMSE) on the train (red), validation (green) and test (blue) set for each function evaluation.</figcaption></figure>

## 6.3 Fixed neural architectures with stacked LSTM and fully connected layers

Fine-tuning deep neural networks is challenging due to the substantial computational effort required. In this series of experiments, the performance of predefined deeper neural architectures with fixed numbers of layers is assessed. In the second batch of experiments (B.2), neural architectures with two stacked LSTM layers ($\alpha^{LSTM} = 2$) and two dense fully connected layers ($\alpha^{FC} = 2$) are investigated (Table 2). For the third batch (B.3), the model architectures are extended to three stacked LSTM layers ($\alpha^{LSTM} = 3$) and two dense fully connected layers ($\alpha^{FC} = 2$) (Table 2). Regarding layer size, search is
performed over $A_i^{LSTM} = [20, 240]$ units for the recurrent layers and $A_j^{FC} = [20, 150]$ units for the fully connected layers. Each dense layer is followed by a leaky ReLU activation function. Given the neural architecture depth, weak $L_2$-regularization is included as a hyperparameter to prevent overfitting, with values in the range $\lambda \in [1e^{-14}, 1e^{-6}]$.

In batches B.2 and B.3, each consisting of five independent trials with a variable number of epochs, BO is applied to jointly tune 7 and 8 hyperparameters, respectively (refer to Table 2). Table 1 provides a summary of the search space for the neural architectures, along with the training configurations and relative parameter bounds (parameters for BO are highlighted in blue).

In each BO experiment, the neural architecture that achieves the lowest RMSE on the validation set $D_{val}$ is selected and used for inference on new, unseen data. To this end, the BO trials are ranked according to their validation performance. As expected, the best-performing trials are obtained when using the largest number of epochs (e.g., 750). From the top-performing BO trial, the neural architectures are ranked and the five best models are used to evaluate their generalization on the test set $D_{test}$. For the batches B.2 ($\alpha^{LSTM} = 2, \alpha^{FC} = 2$) and B.3 ($\alpha^{LSTM} = 3, \alpha^{FC} = 2$), the ranking of best observed feasible points, along with the corresponding hyperparameters from the BO process, is shown in Table 6 and Table 7, respectively. The best hyperparameters achieved in B.2 yield a validation error that is 6.23% lower than that obtained in B.3. Nevertheless, both experiments result in nearly identical errors on the test set. A similar conclusion can also be drawn from the rest models, indicating that increased depth does not confer performance benefits within the explored search space.

<table><caption>Table 6: Ranking of the top five feasible points observed during the best-performing BO trial T.1 in experiment batch B.2 (2 LSTM - 2 FCL)</caption><thead><tr><th>Architecture #</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>N. of units at layer A<sub>1</sub><sup>LSTM</sup></td><td>66</td><td>33</td><td>64</td><td>35</td><td>31</td></tr><tr><td>N. of units at layer A<sub>2</sub><sup>LSTM</sup></td><td>44</td><td>37</td><td>35</td><td>34</td><td>71</td></tr><tr><td>N. of units at layer A<sub>1</sub><sup>FC</sup></td><td>102</td><td>106</td><td>140</td><td>111</td><td>87</td></tr><tr><td>N. of units at layer A<sub>2</sub><sup>FC</sup></td><td>123</td><td>40</td><td>57</td><td>76</td><td>64</td></tr><tr><td>Learning rate &gamma;</td><td>0.00232</td><td>0.0025</td><td>3.0779e-04</td><td>0.0014</td><td>7.7797e-04</td></tr><tr><td>Minibatch &xi;</td><td>231</td><td>231</td><td>93</td><td>93</td><td>99</td></tr><tr><td>D<sub>val</sub> (min objective RMSE)</td><td><strong>0.0059782</strong></td><td>0.0067312</td><td>0.0067510</td><td>0.0067776</td><td>0.0068506</td></tr><tr><td>D<sub>train</sub> (RMSE)</td><td><strong>0.0035021</strong></td><td>0.0036389</td><td>0.0047327</td><td>0.0044772</td><td>0.0045040</td></tr><tr><td>D<sub>test</sub> (RMSE)</td><td><strong>0.0072065</strong></td><td>0.0078754</td><td>0.0076774</td><td>0.0078516</td><td>0.0077298</td></tr></tbody></table>

<table><caption>Table 7: Ranking of the top five feasible points observed during the best-performing BO trial T.1 in experiment batch B.3 (3 LSTM - 2 FCL).</caption><thead><tr><th>Architecture #</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>N. of units at layer A<sub>1</sub><sup>LSTM</sup></td><td>60</td><td>66</td><td>78</td><td>81</td><td>63</td></tr><tr><td>N. of units at layer A<sub>2</sub><sup>LSTM</sup></td><td>41</td><td>27</td><td>31</td><td>22</td><td>44</td></tr><tr><td>N. of units at layer A<sub>1</sub><sup>FC</sup></td><td>22</td><td>56</td><td>74</td><td>65</td><td>39</td></tr><tr><td>N. of units at layer A<sub>2</sub><sup>FC</sup></td><td>130</td><td>99</td><td>61</td><td>149</td><td>22</td></tr><tr><td>N. of units at layer A<sub>1</sub><sup>FC</sup></td><td>123</td><td>100</td><td>31</td><td>51</td><td>129</td></tr><tr><td>Learning rate &gamma;</td><td>0.0046370</td><td>4.6925e-04</td><td>4.9004e-04</td><td>4.8288e-04</td><td>4.7060e-04</td></tr><tr><td>Minibatch &xi;</td><td>209</td><td>144</td><td>184</td><td>204</td><td>164</td></tr><tr><td>D<sub>val</sub> (min objective RMSE)</td><td><strong>0.0063792</strong></td><td>0.0066088</td><td>0.0066679</td><td>0.0067048</td><td>0.0068114</td></tr><tr><td>D<sub>train</sub> (RMSE)</td><td><strong>0.0027985</strong></td><td>0.004154</td><td>0.0041963</td><td>0.0043583</td><td>0.0043682</td></tr><tr><td>D<sub>test</sub> (RMSE)</td><td><strong>0.0072469</strong></td><td>0.007554</td><td>0.0078424</td><td>0.0077534</td><td>0.0079657</td></tr></tbody></table>
The performance of BO versus number of function evaluations for the best-performing trials T.1 of the second B.2 and third B.3 batch of experiments is shown in Fig. 10 and Fig. 11, respectively. As can been seen in Fig. 10, during the initial BO evaluations, the estimated minimum objective function shows a clear decreasing trend, closely following the minimum observed objective value. However, beyond these initial function evaluations, and with the exception of a few points, the GP surrogate model performance becomes unstable, exhibiting significant fluctuations. These may occur due to sudden changes in the GP posterior mean hyper-surface which is a natural consequence of the GP's sensitivity to new observations, especially in highly irregular and noisy objective function landscapes. A different BO result is shown in Fig. 11 where the estimated minimum objective function value i) remarkably coincides with the minimum observed objective function value, and ii) shows a decreasing trend, indicating convergence of the BO algorithm. This suggests that the BO algorithm is effectively leveraging the probabilistic model to guide the search towards the optimal solution.

<figure><img src="image_11.png" /><figcaption>Fig. 10: Minimum objective value versus functions evaluations (T.1, B.2). The horizontal axis represents the iterations of BO algorithm, while the vertical axis shows the prediction accuracy on the validation set (RMSE); (2 stacked LSTM - 2 FCL).</figcaption></figure>

The sequential evaluation points of the BO, for the best-performing trials of B.2 and B.3, considering the sets $D_{train}$, $D_{val}$, $D_{test}$, is depicted in Fig. 12 and Fig. 13, respectively. The horizontal axis shows the number of BO iterations, while the vertical axis represents the prediction accuracy (RMSE) for each function evaluation. Notably, the convergence of the curves for the validation $D_{val}$ and test $D_{test}$ sets highlights two key points: (i) the strong generalization ability of the model, as the validation set reliably reflects performance on unseen data, and, (ii) the absence of bias in BO performance, as the optimization process is not overly reliant on the predictions from the validation set. Additionally, the fluctuations observed in this plot serve as an indicator of the exploratory nature of the BO within the design space. These variations reflect the method's systematic search for optimal solutions across the space of possible configurations.

## 6.4 Varying deep neural architectures with stacked LSTM and fully connected layers

Deeper architectures and larger search spaces offer the potential for discovering highly performant models. Unfortunately, as the dimensionality increases, the search space grows exponentially, making BO significantly more complex. In this section, the primary objective is to evaluate the BO performance
Min objective vs. Number of function evaluations

<figure><img src="image_12.png" /><figcaption>Fig. 11: Minimum objective value versus functions evaluations (T.1, B.3). The horizontal axis represents the iterations of BO algorithm, while the vertical axis shows the prediction accuracy on the validation set (RMSE) (3 stacked LSTM - 2 FCL).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12: The sequential evaluation points during the BO process for the best-performing trial T.1 of the second B.2 experiment batch, showing the prediction accuracy (RMSE) across the training $D_{train}$ (red), validation $D_{val}$ (green), and test $D_{test}$ (blue) sets (2 stacked LSTM - 2 FCL).</figcaption></figure>

within an extensive search space, specifically one that involves jointly tuning 15 hyperparameters (Table 2). The key difference in this experiment, is that the neural architectures are not predefined, and different number of layers can be considered during optimization. Specifically, for each candidate architecture within the search space $\mathcal{X}$, the number of stacked LSTM layers is introduced as a scalar hyperparameter, denoted by $\alpha^{LSTM} \in [1, 4]$. Similarly, the number of dense fully connected layers is treated as a scalar hyperparameter, represented by $\alpha^{FC} \in [1, 6]$. For the size of recurrent and dense layers units, the search is conducted over the intervals $A_t^{LSTM} = [20, 240]$ and $A_j^{FC} = [20, 150]$, respectively. Each dense layer is followed by a leaky ReLU activation function. Additionally, an $L_2$-regularization penalty is applied, whose strength was drawn from $1e^{-9}$ to $1e^{-3}$, to prevent overfitting in deeper models.
<figure><img src="image_14.png" /><figcaption>Fig. 13: The sequential evaluation points during the BO process for the best-performing trial T.1 of the third B.3 experiment batch, showing the prediction accuracy (RMSE) across the training $D_{train}$, validation $D_{val}$, and test $D_{test}$ sets versus the number of function evaluations (3 stacked LSTM - 2 FCL).</figcaption></figure>

els. Note that the search spaces defined in previous experiments (B.1, B.2 and B.3) are subspaces of the broader search space considered in this experiment (B.4). Table 1 summarizes the hyperparameter selection for the deep networks, along with the corresponding training options and the relative bounds of the parameters (BO parameters are indicated in blue).

During the optimization process, not ignoring the units of unused layers (e.g., if $A_i^{LSTM}, i > 2$ are not ignored when $\alpha^{LSTM} = 2$) can significantly degrade BO performance. The GP surrogate model wastes capacity modeling irrelevant dimensions, reducing its accuracy since BO assumes all hyperparameters influence the objective. Unused layer units introduce artificial noise, making it harder to identify true optima. To maintain BO efficiency, dynamic hyperparameter masking is implemented via conditional search spaces by using conditional constraint functions inside Matlab's *bayeopt* function. Thus, inactive variables are not modeled at all, and the GP surrogate focuses only on relevant dimensions. In this way, the optimization trace reflects the true hyperparameter space explored.

The top five best-performing models are selected and assessed for accuracy on the test set $D_{test}$. Table 8 shows the ranking of best observed feasible points, along with the corresponding hyperparameters. The minimum RMSE observed across all objective function evaluations is 0.0061214.
<table><caption>Table 8: Ranking of the top five feasible points observed during experiment batch B.4.</caption><thead><tr><th>Architecture #</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>N. of units at layer A<sub>1</sub><sup>STM</sup></td><td>72</td><td>50</td><td>58</td><td>24</td><td>63</td></tr><tr><td>N. of units at layer A<sub>2</sub><sup>STM</sup></td><td>139</td><td>219</td><td>60</td><td>49</td><td>112</td></tr><tr><td>N. of units at layer A<sub>3</sub><sup>STM</sup></td><td>116</td><td>202</td><td>225</td><td>180</td><td>174</td></tr><tr><td>N. of units at layer A<sub>4</sub><sup>STM</sup></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>N. of units at layer A<sub>1</sub><sup>FC</sup></td><td>98</td><td>123</td><td>73</td><td>133</td><td>115</td></tr><tr><td>N. of units at layer A<sub>2</sub><sup>FC</sup></td><td>103</td><td>95</td><td>100</td><td>79</td><td>74</td></tr><tr><td>N. of units at layer A<sub>3</sub><sup>FC</sup></td><td>133</td><td>95</td><td>135</td><td>148</td><td>146</td></tr><tr><td>N. of units at layer A<sub>4</sub><sup>FC</sup></td><td>84</td><td>93</td><td>104</td><td>69</td><td>91</td></tr><tr><td>N. of units at layer A<sub>5</sub><sup>FC</sup></td><td>-</td><td>87</td><td>-</td><td>141</td><td>102</td></tr><tr><td>N. of units at layer A<sub>6</sub><sup>FC</sup></td><td>-</td><td>133</td><td>-</td><td>129</td><td>-</td></tr><tr><td>Learning rate γ</td><td>2.923e-04</td><td>2.319e-04</td><td>3.392e-04</td><td>8.296e-04</td><td>4.689e-04</td></tr><tr><td>Minibatch ξ</td><td>139</td><td>207</td><td>232</td><td>184</td><td>231</td></tr><tr><td>L<sub>2</sub>-Regularization λ</td><td>2.718e-09</td><td>1.662e-09</td><td>2.636e-09</td><td>2.284e-09</td><td>1.293e-09</td></tr><tr><td>D<sub>val</sub> (min objective RMSE)</td><td><b>0.0061214</b></td><td>0.0063736</td><td>0.0065123</td><td>0.0065722</td><td>0.0065869</td></tr><tr><td>D<sub>train</sub> (RMSE)</td><td><b>0.0024909</b></td><td>0.0032285</td><td>0.0026325</td><td>0.0029533</td><td>0.0031808</td></tr><tr><td>D<sub>test</sub> (RMSE)</td><td><b>0.0072654</b></td><td>0.0072847</td><td>0.0078189</td><td>0.0073120</td><td>0.0075317</td></tr></tbody></table>

The BO performance in relation to the number of function evaluations is depicted in Fig. 14. The BO demonstrates a rapid initial improvement within the first 10 evaluations, characterized by a notable reduction in both the observed and estimated objective values. This behavior suggests that the optimizer quickly identified a promising region within the search space. Subsequently, a plateau phase is observed between evaluations 10 and 80, during which only minor improvements were achieved. Around the 80th evaluation, a new minimum objective value was discovered, followed by another plateau in performance. Several fluctuations in the estimated objective between evaluations 80 and 120 reflect the surrogate model's continued exploration of regions with predicted potential, although these did not yield improvements over the previously observed minimum.

<figure><img src="image_15.png" /><figcaption>Fig. 14: Minimum objective value versus function evaluations. The horizontal axis represents the iterations of BO algorithm, while the vertical axis shows the prediction accuracy on the validation set (RMSE).</figcaption></figure>
The sequential evaluation points of the BO considering the sets $D_{train}$, $D_{val}$, $D_{test}$ is shown in Fig. 15. The horizontal axis represents the iterations of BO algorithm, and the vertical axis represents the prediction accuracy (RMSE) on the train, validation and test sets for each function evaluation. The curve associated to the validation set $D_{val}$ shows the history of objective function values (e.g., RMSE) across all iterations, where each point on the plot represents the objective function value at a specific iteration. Again, one can observe the convergence of the curves associated with the validation $D_{val}$ and test $D_{test}$ sets, which confirms i) the generalization ability of the validation set, and ii) that BO performance is not biased by predictions on the validation set. Again, the variability in the values observed at the evaluation points reflects the method's active exploration of the design space in search of optimal solutions.

<figure><img src="image_16.png" /><figcaption>Fig. 15: The sequential evaluation points of the BO showing the prediction accuracy (RMSE) across the training $D_{train}$, validation $D_{val}$, and test $D_{test}$ sets.</figcaption></figure>

### 6.4.1 Modeling Costs

Execution times vary greatly across different regions of the hyperparameter space due to varying depth, regularization, learning rate, etc. Fig. 16 shows the total function evaluation time, which represents the time spent evaluating the objective function at each iteration. It can be observed that the evaluation times exhibit significant variation, ranging from approximately 1275 s to $3.5 \times 10^4$ s on a workstation equipped with an Intel® Core™ i9-14900K CPU (24 cores, 3.2 GHz) 96 GB RAM.

### 6.4.2 Principal Component Analysis of neural architectures optimized via Bayesian search: Depth and performance trends

Understanding the structure of the resulting search space and the characteristics of discovered architectures remains a key challenge. For this reason, Principal Component Analysis (PCA) [50] is employed which provides a means of reducing the dimensionality of hyperparameter space, enabling visual inspection and qualitative analysis. PCA is a parametric linear dimensionality reduction technique. It consists in finding the hyper-plane that maximizes the variance of the projected data. The principal components can be computed using a singular value decomposition

$$ \tilde{Z}^T \tilde{Z} = V \Sigma^2 V^T \qquad (16) $$
<figure><img src="image_17.png" /><figcaption>Fig. 16: Function evaluation run time versus the number of function evaluations.</figcaption></figure>

where matrix $\tilde{Z}$ is obtained after mean-centering $Z$, and $\Sigma$ and $V$ are the respective (sorted) eigenvalues and eigenvectors matrices. Then, the partial projection matrix $V^{(k)}$ can be obtained by truncating $V$ up to the $k$-eigenvector (column). As such, the projected data is obtained by $\hat{Z}^{(k)} = \tilde{Z}V^{(k)}$.

In the NAS specific case, the data comprises neural architectures characterized by its hyperparameters, e.g., depth (number of layers) and number of layer units. In this analysis, the neural network training hyperparameters, e.g., learning rate, batch size and $L_2$-regularization are not included. The data assembled into a matrix $Z \in \mathbb{R}^{m \times p}$, where $m$ denotes the number of architectures and $p$ the number of hyperparameters. Clearly, data of matrix $Z$ are contained in spaces of high dimension, which makes their analysis not straightforward. Thus, dimensionality reduction is employed so as to transform the matrix $Z$ into a low dimensional representation. We annotate the reduced PCA projection matrix $\hat{Z}^{(2)} \in \mathbb{R}^{m \times 2}$ to explore correlations between architecture complexity and learning efficacy. As such, the respective manifolds can be easily visualized.

The 2D projection is shown in Fig. 17 (a) following a color-scale that indicates the corresponding BO objective function RMSE on a scale from $\approx 6 \times 10^{-3}$ to $12 \times 10^{-3}$. Thus, architecture performance is encoded using a color map, with darker blue colors indicating better performance (lower validation error). It should be noted that the RMSE values are not used for the PCA, they are shown here only for visualization purposes.

The PCA projection reveals a structured and nonlinear trajectory of architectures in the reduced space. A distinct curved pattern is observed, suggesting that the underlying hyperparameter space is organized along a low-dimensional manifold, despite the linear nature of PCA. The depth of each network (total number of LSTM and fully connected dense layers) is indicated by numerical labels next to the corresponding point. Shallower networks (depths 2 to 4) are largely located in the left portion of the PCA space (PC1 < -2), while deeper networks (depths 7 to 10) cluster more tightly in the center-right region (PC1 from 0.5 to 2.5). This distribution suggests that the optimization process transitioned from shallow to deeper architectures, likely driven by performance feedback, as BO focused on promising configurations. Most of the best architectures (deep blue) are clustered in the center-right region, especially around PC1 $\approx [0.5, 2]$, corresponding to networks with depths between 7 and 9. Conversely, architectures with poor performance (orange/red points) are more dispersed and appear frequently among both shallow and overly deep models, such as certain configurations with depths 9 and 10. This indicates that while deeper networks are generally better, there is a sweet spot in the depth range where performance peaks. The increasing density of points in the high-performing PCA region
implies that BO effectively exploited areas of the hyperparameter space with promising results more heavily over time. The arc-like formation and clustering pattern indicate that BO may have converged toward an optimal subspace characterized by moderate depth and certain topological features.

<figure><img src="image_18.png" /></figure>

To gain deeper insights into the distribution and structure of neural architectures discovered via BO, K-means clustering is applied (with $k = 4$) to the architectures projected into the two-dimensional PCA space (Fig. 17) (b). Each point in the space represents a neural architecture encoded by its topological configuration. The depth of each network (total number of LSTM and fully connected dense layers) is indicated by numerical labels next to the corresponding point. The resulting clusters reveal distinct families of architectures, which align closely with network depth and likely correspond to different regions of the performance landscape. Cluster 1 (blue) located in the positive PC1 and PC2 region dominated by shallow to mid-depth networks, composed mostly of architectures with depths 6–7. Cluster 2 (orange) is tightly packed in the mid-left region and similarly dominated by shallow networks. Its concentrated nature may point to a local optimum with suboptimal performance. Cluster 3 (yellow), located in the bottom right quadrant, is characterized by deeper architectures (depths 8–10) and reflects a shift in the architecture search towards increased complexity. Cluster 4 (purple) contains moderately deep networks (depths 7–9) with high density, aligning with regions previously observed to yield superior performance. This segmentation confirms that the architecture space exhibits structured, non-random variability, which PCA can effectively reveal and K-means can exploit. The separation between clusters supports the hypothesis that neural architecture design space contains distinct subregions corresponding to different design motifs and performance regimes. The clustering has effectively segmented the space into: i) Simpler, possibly suboptimal architectures (clusters 1, 2), and ii) Advanced, diverse architecture families (clusters 3, 4). These clusters can be further leveraged for *pruning, meta-modeling, or adaptive sampling* in future iterations of the NAS pipeline.

The choice of a 2D visualization using PCA limits the insight into the complete manifolds. Indeed, as Fig. 18 shows, the percentage of the total variance explained by the first two principal components is 48.55% (35.44% by PC1 and 13.11% by PC2). It can be seen that the remaining principal components do contain some non-negligible variance of the data. Table 9 presents the contribution of hyperparameters to principal components PC1 and PC2. The number of LSTM layers and fully connected layers
dominate PC1, contributing 17.94% and 17.88%, respectively. For PC2, the number of units in fully connected layers 1 and 6 contribute most, with 29.35% and 21.93%, respectively.

<figure><img src="image_19.png" /><figcaption>Fig. 18: Percentage of the total variance explained by each principal component and cumulative explained variance.</figcaption></figure>

<table><caption>Table 9: Contribution of hyperparameters to principal components PC1 and PC2.</caption><thead><tr><th></th><th>α<sup>LSTM</sup></th><th>A<sub>1</sub><sup>LSTM</sup></th><th>A<sub>2</sub><sup>LSTM</sup></th><th>A<sub>3</sub><sup>LSTM</sup></th><th>A<sub>4</sub><sup>LSTM</sup></th><th>α<sup>FC</sup></th><th>A<sub>1</sub><sup>FC</sup></th><th>A<sub>2</sub><sup>FC</sup></th><th>A<sub>3</sub><sup>FC</sup></th><th>A<sub>4</sub><sup>FC</sup></th><th>A<sub>5</sub><sup>FC</sup></th><th>A<sub>6</sub><sup>FC</sup></th></tr></thead><tbody><tr><td>PC1 (%)</td><td><strong>17.94</strong></td><td>0.02</td><td>7.94</td><td>10.95</td><td>0.38</td><td><strong>17.88</strong></td><td>0.20</td><td>1.10</td><td>13.82</td><td>14.31</td><td>10.33</td><td>5.12</td></tr><tr><td>PC2 (%)</td><td>2.91</td><td>8.56</td><td>3.75</td><td>5.72</td><td>4.00</td><td>6.27</td><td><strong>29.35</strong></td><td>0.001</td><td>7.11</td><td>0.02</td><td>10.36</td><td><strong>21.93</strong></td></tr></tbody></table>

## 7 Random Search Experiments

This section investigates the effectiveness of random search in discovering performant architectures over the same 15-dimensional search space explored using BO. By analyzing the distribution, performance, and structural characteristics of randomly sampled architectures, we aim to establish foundational insights and reference baselines to more sophisticated NAS strategies such as BO.

In the context of NAS, random search operates by uniformly sampling points (architectures and hyperparameters) at random in each dimension from a predefined search space without bias toward previously evaluated configurations (Algorithm 2). Despite its simplicity, it has been shown to outperform more complex search strategies under certain conditions, particularly when the number of tunable dimensions is large and the performance landscape is irregular or noisy. This is supported by theoretical results such as those by Liam et al. [51], Yu et al. [52] and Bergstra et al. [5] demonstrating that, in high-dimensional hyperparameter spaces, random search can identify competitive solutions faster than grid-based approaches by exploring more diverse regions of the space.

When applied to neural architecture design, random search offers several advantages: it is easy to parallelize, does not require surrogate models or acquisition functions, and avoids the risk of overfitting to a biased exploration strategy. However, it also suffers from key limitations, including the lack of guided refinement and inefficiency in sample utilization, which can be critical in compute-constrained NAS settings. As such, random search is often used as a *natural baseline or bootstrapping mechanism* in hybrid NAS workflows, particularly when combined with more advanced methods such as BO, reinforcement learning, or evolutionary algorithms.
In the random search experiment, 120 neural architectures were uniformly drawn from the design space. Table 10 presents the five highest-ranked architectures discovered, ordered by their minimum validation error ($D_{val}$, RMSE). One can observe that the number of units per layer varying significantly (e.g., LSTM: 34-240 units; FC: 29-147 units). Both deep LSTM configurations (e.g., architectures #1 and #3) and shallower LSTM configurations (e.g., architectures #2 and #5) appear among the top performers, indicating multiple viable structural patterns. The learning rates range from 9.435 × 10⁻⁵ to 5.531 × 10⁻⁴. Minibatch sizes span from 65 to 180, while L₂-regularization coefficients vary over several orders of magnitude (10⁻¹⁰ to 10⁻⁷), demonstrating the effectiveness of diverse hyperparameter settings in achieving competitive results. Architecture #1 achieves the lowest validation error (RMSE = 0.0068867) however with the worst generalization performance among the five highest-ranked models (test RMSE = 0.0096811). Notably, architecture #5, while ranked fifth on validation, achieves the lowest test RMSE (0.007553), indicating superior generalization performance. All architectures exhibit low training RMSE (≤ 0.0043), suggesting effective optimization across the board. These results demonstrate that high-performing architectures can differ substantially in structure and hyperparameter configuration. The diversity of these top models underscores the capability of random search to uncover a broad set of competitive solutions within the neural architecture space. Interestingly, granting random search the same search space and computational budget, random search found relatively comparative models with BO. This can be attributed to the function $f(\mathbf{x})$ which has a low effective dimensionality; essentially, $f(\mathbf{x})$ is more sensitive to changes in some dimensions than others [53].

<caption>**Algorithm 2:** Random Search for NAS</caption>

**Input:** Search space $\mathcal{X}$; hyperparameters $\mathbf{x} = (x_1, \dots, x_p)$; parameter search space ranges $S_i, i = 1, \dots, p$; expensive evaluator $\mathcal{E}(\mathbf{x}, \zeta)$ returning performance $f(\mathbf{x})$ at fidelity $\zeta$ (e.g., epochs); total iterations $N$.
**Output:** Architectures $\mathbf{x}$ and their performance $\mathcal{E}$.
1 &nbsp;&nbsp;Set $n \leftarrow 1$;
2 &nbsp;&nbsp;**while** $n \le N$ **do**
3 &nbsp;&nbsp;&nbsp;&nbsp;Sample hyperparameters $x_i \leftarrow \text{Uniform}(S_i)$;
4 &nbsp;&nbsp;&nbsp;&nbsp;Create architecture structure and training hyperparameters:
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\mathbf{x}_n : \text{arch}.x_i \leftarrow x_i$$;
6 &nbsp;&nbsp;&nbsp;&nbsp;Evaluate architecture performance $\mathcal{E}(\mathbf{x}_n, \zeta)$;
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$n \leftarrow n + 1;$$
8 &nbsp;&nbsp;**return** $\mathbf{x}, \mathcal{E}(\mathbf{x}, \zeta)$
<table><caption>Table 10: Ranking of the top five neural architectures observed during random search.</caption><thead><tr><th>Architecture #</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>N. of units at layer A<sub>1</sub><sup>LSTM</sup></td><td>69</td><td>224</td><td>34</td><td>158</td><td>203</td></tr><tr><td>N. of units at layer A<sub>2</sub><sup>LSTM</sup></td><td>213</td><td>-</td><td>48</td><td>183</td><td>-</td></tr><tr><td>N. of units at layer A<sub>3</sub><sup>LSTM</sup></td><td>223</td><td>-</td><td>124</td><td>240</td><td>-</td></tr><tr><td>N. of units at layer A<sub>4</sub><sup>LSTM</sup></td><td>217</td><td>-</td><td>234</td><td>168</td><td>-</td></tr><tr><td>N. of units at layer A<sub>1</sub><sup>FC</sup></td><td>96</td><td>60</td><td>37</td><td>138</td><td>147</td></tr><tr><td>N. of units at layer A<sub>2</sub><sup>FC</sup></td><td>-</td><td>111</td><td>84</td><td>48</td><td>144</td></tr><tr><td>N. of units at layer A<sub>3</sub><sup>FC</sup></td><td>-</td><td>48</td><td>97</td><td>-</td><td>65</td></tr><tr><td>N. of units at layer A<sub>4</sub><sup>FC</sup></td><td>-</td><td>29</td><td>94</td><td>-</td><td>-</td></tr><tr><td>N. of units at layer A<sub>5</sub><sup>FC</sup></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>N. of units at layer A<sub>6</sub><sup>FC</sup></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Learning rate γ</td><td>4.540e-04</td><td>9.435e-05</td><td>1.512e-04</td><td>4.595e-04</td><td>5.5310e-04</td></tr><tr><td>Minibatch ξ</td><td>65</td><td>132</td><td>100</td><td>108</td><td>180</td></tr><tr><td>L<sub>2</sub>-Regularization λ</td><td>3.921e-09</td><td>2.879e-10</td><td>4.682e-09</td><td>5.536e-10</td><td>8.169e-07</td></tr><tr><td>D<sub>val</sub> (min objective RSRMSE)</td><td><strong>0.0068867</strong></td><td>0.0069133</td><td>0.0069454</td><td>0.0070759</td><td>0.0071001</td></tr><tr><td>D<sub>train</sub> (RMSE)</td><td>0.0012773</td><td>0.0043167</td><td>0.0041166</td><td><strong>0.00077318</strong></td><td>0.0044249</td></tr><tr><td>D<sub>test</sub> (RMSE)</td><td>0.0096811</td><td>0.0075697</td><td>0.0076971</td><td>0.008524</td><td><strong>0.007553</strong></td></tr></tbody></table>

## 7.1 Principal Component Analysis of neural architectures optimized via random search: Depth and performance trends

For the random search experiment, the PCA was conducted on data representing neural architectures characterized by their structural hyperparameters (e.g., depth, number of units per layer) and training-related hyperparameters including learning rate, batch size, and L₂-regularization. Fig. 19 (a) presents a PCA projection which aims to visualize the structural diversity of architectures evaluated during the search. Each point in the scatter plot represents a distinct architecture, projected into a 2D space formed by the first two principal components of the design space (the x-axis and y-axis represent the first and second principal components, respectively). The points are color-coded by performance (RMSE), while the accompanying numerical labels indicate network depth in terms of the number of layers. Importantly, well-performing architectures are distributed throughout the space, suggesting that competitive performance can be achieved by a variety of structural configurations rather than being confined to a single architectural pattern. Poorly performing architectures (in red to yellow) are more isolated and less frequent. The best-performing architectures exhibit a range of depths, with moderate-depth networks (4-7 layers) showing consistently good performance. Very deep networks (9-10 layers) appear at the right part of the space but do not show superior performance, suggesting that simply increasing depth does not guarantee better results in this search space.

The random search effectively explored a broad and diverse subset of the architecture space, as evidenced by the wide spread of points in the PCA projection. The existence of high-performing architectures in multiple disconnected regions implies that there may be multiple optima or architectural motifs yielding low RMSE. The clustering of blue-colored, low-RMSE architectures in specific PCA regions highlights the potential of random search to uncover performant models, while also emphasizing the importance of architectural structure beyond just depth.

Fig. 19 (b) presents the results of applying K-means clustering (with $k = 4$) to neural architectures that were evaluated during random search. The architectures are first projected into a two-dimensional space using PCA. The clustering was then performed in this reduced space to identify distinct structural groupings among architectures. Each point in the scatter plot corresponds to a unique architecture. Color and marker shape represent the assigned cluster label (from 1 to 4). Cluster centroids are marked by large black ‘x’ symbols. Each architecture is annotated with its depth (i.e., number of layers), providing insight into how network depth is distributed across clusters. The data naturally segregated into four distinct regions of the PCA space, suggesting that architectures can be grouped based on
latent structural characteristics. Clusters are well-separated in the PCA space, especially clusters 1 (blue, upper region), 2 (orange, right), and 4 (purple, left); Cluster 1 (blue) contains a mix of mid-depth architectures (e.g., depths 5-8) and is located in the upper region of the PCA space. Cluster 2 (orange) is dominated by deeper architectures (e.g., depth 8-10), occupying the right region. This suggests that deeper architectures may share common design features leading to this spatial grouping. Cluster 3 (yellow) primarily composed of mid-depth architectures (e.g., depths 5-6), positioned in the lower-central area. Cluster 4 (purple) mostly contains shallow architectures (e.g., depths 2-5), located in the left-central region.

The cluster centroids are centrally positioned within their respective groups, reinforcing that each cluster is reasonably compact and distinct. The relative distances between cluster centroids indicate inter-cluster dissimilarity. Notably, clusters 2 and 4 are the most distant from each other in PCA space, suggesting that architectures in these groups are structurally very different. In contrast, clusters 3 and 4 are positioned more closely, implying more similarity in their architectural characteristics. The clustering analysis reveals that architectural diversity explored by random search is non-uniformly distributed-some regions of the design space concentrate similar depth and structural patterns. The clear separation of clusters in the PCA space suggests that latent architectural features (not just depth) significantly contribute to network behavior and design similarity.

<figure><img src="image_20.png" /><figcaption>Fig. 19: Principal component projection of the matrix $Z \in \mathbb{R}^{m \times p}$.</figcaption></figure>

Moreover, Fig. 20 summarizes the variance explained by the principal components obtained in the PCA analysis. The scree plot (left) shows that the first two components account for approximately 40% of the total variance, with PC1 and PC2 contributing around 23% and 18%, respectively. Subsequent components explain progressively smaller proportions, with individual contributions falling below 3% beyond the tenth component. The cumulative explained variance plot (right) indicates that the first five components capture over 70% of the variance, while approximately ten components are required to exceed the 90% threshold. These results highlight that although two-dimensional projections provide a useful visualization of the architecture space, a higher number of components is necessary to preserve most of the variability in the dataset.
<figure><img src="image_21.png" /><figcaption>Fig. 20: Percentage of the total variance explained by each principal component and cumulative explained variance.</figcaption></figure>

## 8 LSTM-based metamodeling for damage evolution

This section presents several examples of the LSTM-based metamodel for predicting the damage evolution. To illustrate its capabilities, the optimal neural architecture A.1 from the best-performing trial T.1 of the second experiment batch B.2 (see Table 6) obtained in BO is used for comparison with FEM simulations.

The prediction performance of the optimal neural architecture on the training, validation and test sets is shown in Fig. 21, Fig. 22 and Fig. 23, respectively. To generate these scatter-type plots for each set, points from all sequences are gathered into a single set, and then sorted by their damage values. It can been observed that the neural network predictions on the training set closely match the actual values, demonstrating that it has effectively learned the underlying patterns in the data. Despite a small number of outliers, the results on the validation and test sets show that the model generalizes well from the training to test data points with high accuracy.

Proceeding further, the quality and comprehensiveness of the sampling data used to train, validate, and test the neural networks are examined. Specifically, this experiment aims to determine how accurately and thoroughly this data covers the entire range of possible scenarios or conditions (the design space D) that the neural network is expected to understand and work with. New points can expose structural dynamics or nonlinearities not captured during training, providing valuable feedback about the strengths and weaknesses of the neural network model. Thus, the neural network predicted stiffness loss evolution curves are compared against the reference FE results for eight new data points randomly sampled from the design space (see Fig. 24 and Fig. 25). The neural network achieves high overall accuracy and effectively learns diverse types of nonlinear stiffness loss evolution curves. However, it performs poorly on a small subset of points, such as points 4 and 5 (see Fig. 24 and Fig. 25), where the FEM reference values are overestimated by approximately 100% and 50% at the final loading state, respectively. One possible reason may be the insufficient representation in training data where certain regions of the design space D may not be well-covered. Thus, the model learns the dominant patterns but fails in sparsely sampled regions of the data space. Additionally, the absence of constraints allows initial stiffness loss values slightly above zero (see point 8 in Fig. 25). Interestingly, the predicted stiffness loss curves conform to monotonically increasing functions and is captured by the neural network without any explicit monotonicity constraints.

Additionally, the LSTM-based metamodel is illustrated by using several combinations of sampling points (design variables in $R^9$). First, the metamodel is constructed to represent damage evolution for random sets of material orientation angles and ply drop geometries under various loading directions. It is important to note that these sampled points are randomly drawn from the design space and may not
<figure><img src="image_22.png" /><figcaption>Fig. 21: Prediction performance on the train set $D_{\text{train}}$ of the optimal LSTM neural network (A.1, T.1, B.2).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 22: Prediction performance on the validation set $D_{\text{val}}$ of the optimal LSTM neural network (A.1, T.1, B.2).</figcaption></figure>

correspond to actual manufacturing scenarios. Figure 26 depicts the metamodel (surface plot) obtained using a triangulation-based cubic interpolation algorithm applied to the neural network’s predicted response data. The interpolation is performed using Matlab’s *griddata* function, ensuring that the interpolated surface passes through the data points with $C^2$ continuity. The plots illustrate the damage evolution as a function of $u_1$ and $u_3$, with constant $u_2 = 0$ (see Fig. 28 for the related variables). The LSTM-based metamodel effectively captures the damage evolution in laminates with random layup configurations. By analyzing a specific set of structural variables, one can clearly identify regions where the structure remains largely intact, detect transitional zones with abrupt jumps in damage evolution, and reveal systematic patterns that could potentially be transformed into rules or predictive models. The LSTM-based metamodel demonstrate the rich diversity of nonlinear damage evolution behaviors in the structure, showcasing the model’s strong capability for detecting and predicting damage progression.
<figure><img src="image_24.png" /><figcaption>Fig. 23: Prediction performance on the test set $D_{test}$ of the optimal LSTM neural network (A.1, T.1, B.2).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 24: Stiffness loss evolution curves. Comparison between the reference FE solutions and the neural network predictions for unseen points 1-4. Note that the vertical axes are not common to all plots (A.1, T.1, B.2).</figcaption></figure>

To further demonstrate the capabilities of the LSTM-based metamodel, it is insightful to visualize how ply drop material orientations influence the damage evolution under specific loading conditions and ply geometries. To achieve this, Latin hypercubes is used to generate a set of randomly distributed points within the design space, corresponding to different ply material orientations. Thus, Latin hypercube sampling with $d = 3$ dimensions and $SP = 400$ points is employed. Fig. 27 illustrates slice plots on different sectional planes within the material orientation design space defined by the angles $\theta_1, \theta_2, \theta_3$. It should be noted that for these plots, random loading conditions were applied, and the maximum
<figure><img src="image_26.png" /><figcaption>Fig. 25: stiffness loss evolution curves. Comparison between the reference FE solutions and the predicted damage function DSL₃ for unseen points 5-8. Note that the vertical axes are not common to all plots (A.1, T.1, B.2).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 26: Metamodel of damage evolution given as a function of $u_1$ and $u_3$ at constant $u_3 = 0$; $\theta_1 = 20.6^\circ, \theta_2 = 29.3^\circ, \theta_3 = 85.8^\circ$; (A.1, T.1, B.2).</figcaption></figure>

loading state was considered. In general, the LSTM-based metamodel assists in identifying the most influential material orientation angles that significantly affect the variation in the damage evolution function. Additionally, it highlights regions of interest within the design input space, providing valuable insights into the relationship between material orientation and damage behavior.
<figure><img src="image_28.png" /><figcaption>Fig. 27: Slices of the maximum damage value are shown with respect to the material orientation angles $\theta_1, \theta_2, \theta_3$. The corresponding ply thicknesses are $t_1 = 50.71, t_2 = 43.94, t_3 = 65.46$. For generating the graph, random loading conditions are applied, with the maximum loading state given as $u_1 = 0.50, u_2 = 0.10, u_3 = -0.165$; (A.1, T.1, B.2).</figcaption></figure>

## 9 Conclusions

To achieve optimal neural architecture design, NAS utilizing BO with conditional constraints is employed in a one-step process that explores the joint search space of neural architectures and hyperparameters. Experimental results demonstrated that the objective function landscape in BO for NAS is highly irregular and multi modal, with no clearly defined local minima. Unlike traditional optimization problems, this landscape presents significant challenges for the efficient identification of optimal architectures. The irregularity stems from the high-dimensional and non-convex nature of the search space. Moreover, it is revealed that various combinations of neural architectures and hyperparameters can result to small validation loss values of the objective function.

Given the computational constraints in the NAS experiments (i.e., number of epochs), the best-performing models in all cases were obtained at the maximum number of epochs (e.g., 750), where the training progress had nearly reached a plateau. Experiments showed that the relationship between network depth and performance is not strictly linear. As the neural network depth was progressively increased from shallow (B.1) to deeper architectures (B.4), the optimal performance was achieved in experiment B2 (see Table 2). The best-performing architecture consisted of two stacked LSTM layers followed by two fully connected layers, achieving an RMSE of 0.0059782 on the validation set. Moreover, all deeper models outperformed the baseline model (B.1). The results indicate that while increasing the depth of the network can enhance performance up to a certain point, the additional complexity introduced by deeper architectures does not always provide substantial benefits beyond an optimal depth. This suggests that there is a trade-off between model complexity and performance gains, and that excessively deep networks may lead to diminishing returns. The findings highlight the importance of carefully balancing depth and complexity in neural network design. Even though, BO provides a powerful framework for exploring this balance, it is crucial to consider computational constraints and the potential for diminishing returns with increased depth.

Additionally, random search is employed as a baseline method for neural architecture exploration to benchmark the performance of BO. While BO achieved superior overall results, random search still identified a diverse set of neural architectures with competitive performance.

Finally, it was demonstrated, that the optimal LSTM-based surrogate model can be used to generate complex mappings from discrete sequences of input variables to discrete nonlinear structural responses in near-real time.
The results motivate extending the study of deep recurrent neural networks to more complex delamination prediction scenarios, including the modeling of initial delamination patterns, multi-layer delamination, and the consideration of various boundary conditions. The proposed method can be also extended to other loading types, e.g., dynamic, fatigue, etc.

Future research in NAS should prioritize the development of advanced deep learning architectures that explicitly incorporate distinct pathways for processing static and dynamic variables. Such models are crucial for capturing the unique temporal dependencies and contextual significance associated with each variable type. Several architectural paradigms can be explored, including dual-stream LSTM networks, feature-fusion LSTM architectures, context-aware LSTM (where static variables are embedded as context vectors), attention-based LSTM with static-guided attention mechanisms, modular LSTM with adaptive routing strategies, graph-based LSTM enabling static-dynamic interactions, and Mixture of Experts (MoE) frameworks tailored for heterogeneous input processing. These architectures are designed to effectively integrate time-series and static features by dynamically modulating the influence of static variables at each time step, based on the evolving context within the sequence. This adaptive weighting mechanism enables the model to selectively emphasize different aspects of the static input, thereby enhancing the network's ability to capture complex interactions and improving its predictive performance in tasks involving heterogeneous data.

An obvious next step to enhance the effectiveness of NAS, is to expand the range of hyperparameters explored during the search process. A broader exploration of the hyperparameter space increases the likelihood of discovering optimal configurations and improves the overall quality of the search. This can be accomplished by developing more efficient sampling strategies, incorporating prior knowledge, or leveraging meta-learning techniques to mitigate the effects of high-dimensional search spaces and improve the scalability of BO for NAS. Another important area of focus is the fine-tuning of key hyperparameters in the GP kernel, such as those in the Matern 5/2 kernel. Properly adjusting these hyperparameters ensures that the GP accurately models the objective function, leading to more efficient and effective optimization. Exploring various acquisition functions, such as Upper Confidence Bound (UCB) and Lower Confidence Bound (LCB), can also provide valuable insights into their effectiveness for NAS. These functions balance exploration and exploitation differently, and their performance may vary depending on the specific problem and search space. Finally, defining a suitable distance metric to quantify the similarity or dissimilarity between neural architectures is essential. When combined with an appropriate kernel for the GP, this metric can significantly improve the modeling of the search space. Relevant work in this area includes studies by Kandasamy et al. [27] and Nguyen et al. [30], which provide a solid foundation for further advancements.

# 10 Acknowledgments

This work is supported by the Composite Innovation Openlab (COMPINNOV), a common laboratory between CETIM, LMPS and GeM Nantes. OF acknowledges the financial support from this entity.

# 11 Appendix

## A Overview of the method

A concise overview of the overall framework is presented as a step-by-step procedure comprising two main components: (i) the mechanical modeling of cohesive zone delamination in the laminate structure using the finite element method, and (ii) the automated design of LSTM architectures through NAS. The overarching concept encompasses the following implementation steps:

### Offline:

* Step 1. Definition of the design space $D$, e.g., loadings, mechanical properties and the reference parametric geometry of the investigated composite structure.

* Step 2. Design of Experiments, utilizing Latin hypercubes.
<figure><img src="image_29.png" /><figcaption>Fig. 28: Diagram and definitions of the parametric ply drop laminate model, including geometrical ($t_1, t_2, t_3$), material orientations ($\theta_1, \theta_2, \theta_3$) and loading ($u_1, u_2, u_3$) design variables.</figcaption></figure>

• Step 3. Damage analysis of the 2D parametric ply-drop model using nonlinear FEM by employing the Surface-Based Cohesive Behavior approach in Abaqus [54].

• Step 4. Establishment of a damage metric delineating the structural performance degradation due to delamination in the cohesive zones. This index is determined within the nonlinear finite element method utilizing Linear Perturbation Analysis (LPA) at several loading increments.

• Step 5. Derivation of stiffness loss evolution functions acquired from the LPA (sequential data) used in the LSTM.

• Step 6. Design of LSTM architectures and hyperparameter tuning using NAS.

**Online:**

• Step 7. Generation of a damage metamodel based on the LSTM predictions to establish a mapping from the sequence of input variables to the nonlinear structural response, e.g., $x^{(s)} \xrightarrow{\text{metamodel}} y^{(s)}$ in near real time (sequence to sequence regression).

## B Mechanical modeling of cohesive zone delamination

### B.1 Ply drop parametric model

The parametric model is depicted in Fig. 28 illustrating the geometric design variables for the ply drop thicknesses $t_i$ and fiber orientation angles $\theta_i$, where $i = 1, 2, 3$. Notably, the tapered angle and the length of the cohesive zones remain fixed. Additionally, the model features loading variables $u_1$ and $u_2$, which prescribe displacements, along with $u_3$, which dictates rotation. These loading conditions are applied to the right-most edge of both plies in the thin ply drop section.

### B.2 Ply drop FEM model

To enhance computational efficiency, a parametric 2D ply drop model using Surface-Based Cohesive Behavior (SBCB) [54] was implemented in a previous work [40]. SBCB offers a consistent methodology for modeling cohesive connections with negligibly small interface thickness, employing the traction-separation constitutive law. Cohesive zone modeling enables the prediction of damage initiation and propagation without requiring a pre-defined crack. Within the SBCB framework, cohesive behavior is realized as a surface interaction property and a pure master-slave formulation is imposed for surfaces exhibiting cohesive behavior. Given that delamination is the predominant failure mode in tapered laminates [55], intralaminar damage was neglected in the finite element analysis.

The reference 2D model consists of 2300 bilinear plane stress quadrilateral elements (of type CPS4) with a total number of 2563 nodes (Fig. 29 and Fig. 30). The element size in the direction of damage propagation is 0.1 mm. Note, that in Abaqus/Standard it is not possible to account for the different fiber orientations of the plies in the 2D case. Thus, a transformation of the material axes is required to
determine the full anisotropic tensor. Nevertheless, some three-dimensional effects, such as anticlastic (saddle-shaped) curvature or elastic coupling are omitted and besides, not all loading conditions can be applied [40]. Its worthy to emphasized that cohesive failure modeling of 3D tapered laminates is a straight-forward task.

<figure><img src="image_30.png" /><figcaption>Fig. 29: Geometry of a benchmark 2D equivalent ply drop model using SBCB.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 30: Mesh of a benchmark 2D equivalent ply drop model using SBCB.</figcaption></figure>

A displacement-controlled procedure is applied since it is more stable than load control for nonlinear analysis. The model is subjected to prescribed displacements imposed at the thin ply drop section while the kinematic constraints are applied to the thick section. Material properties of Celstran/Ticona/PA66 (provided by CETIM) are used in the model: $E_{11} = 38600$ MPa, $E_{22} = E_{33} = 8800$ MPa, $\nu_{12} = \nu_{13} = \nu_{23} = 0.3$, $G_{12} = G_{13} = G_{23} = 3500$ MPa. The elastic properties of the interface material are defined using an uncoupled traction-separation behavior with interface stiffness of $K_{nn} = K_{ss} = K_{tt} = 1 \times 10^6$ N/mm³. The relevant material data are as follows: $t_n^o = t_s^o = t_t^o = 40$ MPa represent the nominal cohesive strengths in the normal and two shear directions, respectively, $G_n^c = 2.2$ N/mm, $G_s^c = G_t^c = 3.0$ N/mm denote the critical fracture energies along the normal and the first and second shear directions, respectively. Additionally, $n = 1.57$ is the semi-empirical criterion exponent applied to delamination initiation and growth [56].

## B.3 Damage accumulation due to delamination in the cohesive zones

The objective is to define an interpretable damage index capable of effectively characterizing performance degradation in ply drop laminates. The proposed index should quantify the structural integrity loss associated with damage accumulation (delamination) within cohesive zones. To achieve this, the stiffness degradation is evaluated along a set of predefined reference directions under an applied loading state. Accordingly, a natural and physically meaningful criterion is formulated based on Directional Stiffness Loss (DSL), and evaluated within a Linear Perturbation Analysis (LPA) framework. The LPA framework decomposes a general nonlinear simulation into S incremental LPA steps. Each LPA step consists of three sequential stages: (i) a nonlinear loading step, (ii) a subsequent unloading step, and (iii) a series of three linear perturbation steps. In the s-th loading step, the structure transitions from an unloaded state $G^{(0)}$ to a deformed state $G^{(s)}$ by applying a fraction of the total load, defined as:

$$ \hat{u}_i^{(s)} = u_i \cdot \frac{s}{S} \qquad (17) $$

for $s \in [1, S]$. Following the loading phase, an unloading step is performed by applying an equal and opposite load, thereby returning the structure to its unloaded configuration $G^{(0)}$. Subsequently, three independent linear perturbation steps are carried out, each perturbing the structure from the state $G^{(0)}$ to distinct perturbed configurations $G_i^{(s)}$, for $i = 1, 2, 3$. Each perturbation corresponds to a prescribed load $\tilde{u}_i = 0.008$ (in either millimeters or radians), typically representing translational perturbations in the x- and y-directions and a rotational perturbation, respectively. After each perturbation, the structure is restored to the unloaded state $G^{(0)}$ in preparation for the subsequent LPA step, as illustrated in Fig. 31 and Fig. 32. Following each perturbation $\tilde{u}_i$, the corresponding structural stiffness K is evaluated based
<figure><img src="image_32.png" /><figcaption>Fig. 31: Illustration of delamination failure for a simulation that has reached the $G^{(8)}$ state (i.e., after the nonlinear loading). Displacement fields $u_1$ and $u_2$ are shown on the left and right sides respectively [41].</figcaption></figure>

on the secant modulus (defined as the slope of the force-displacement response). For the s-th LPA step, the stiffness in each perturbation direction is computed as:

$$K_i^{(s)} = \frac{R_i}{\tilde{u}_i} \qquad (18)$$

with the corresponding reaction force or moment $R_i$ computed as the sum of all nodal reactions along the left-most edge of all plies in the thick ply drop region. Additionally, a reference stiffness $K_i^{(0)}$ is calculated in the initial intact configuration of the structure, prior to the application of any nonlinear load and finally, the DSL is estimated as:

$$DSL_i^{(s)} = 1 - \frac{K_i^{(s)}}{K_i^{(0)}} \qquad (19)$$

Thus, the accumulation of damage is represented by the DSL functions, which are defined as a sequence of s data points corresponding to each LPA step. It is evident that, by employing these DSL functions as design criteria, a specific threshold of stiffness degradation (e.g., 10%) can be established as a design limit.

## References

1. George E. Dahl, Dong Yu, Li Deng, and Alex Acero. Context-dependent pre-trained deep neural networks for large-vocabulary speech recognition. *IEEE Transactions on Audio, Speech, and Language Processing*, 20:30–42, 2012.

2. Alex Krizhevsky, Ilya Sutskever, and Geoffrey E. Hinton. Imagenet classification with deep convolutional neural networks. *Commun. ACM*, 60(6):84–90, may 2017.

3. Ilya Sutskever, Oriol Vinyals, and Quoc V. Le. Sequence to sequence learning with neural networks, 2014.

4. R E Bellman. *Adaptive Control Processes: A Guided Tour*. 1961. First introduction of curse of dimensionality.

5. James Bergstra and Yoshua Bengio. Random search for hyper-parameter optimization. *Journal of Machine Learning Research*, 13(Feb):281–305, 2012.

6. Martin Wistuba, Ambrish Rawat, and Tejaswini Pedapati. A survey on neural architecture search, 2019.

7. Thomas Elsken, Jan Hendrik Metzen, and Frank Hutter. Neural architecture search: A survey, 2019.

8. Bowen Baker, Otkrist Gupta, Nikhil Naik, and Ramesh Raskar. Designing neural network architectures using reinforcement learning. *5th International Conference on Learning Representations, ICLR 2017 - Conference Track Proceedings*, 2019.

9. Barret Zoph and Quoc V. Le. Neural architecture search with reinforcement learning, 2017.
<figure><img src="image_33.png" /><figcaption>Fig. 32: Illustration of linear perturbation steps applied to a structure after the $G^{(8)}$ state has been reached. Different components ($u_i$) are shown for each case with a scale factor of 30 [41].</figcaption></figure>

10. J David Schaffer, Darrell Whitley, and Larry J Eshelman. Combinations of genetic algorithms and neural networks: A survey of the state of the art. *International Workshop on Combinations of Genetic Algorithms and Neural Networks*, pages 1–37, 1992.

11. H. KITANO. Designing Neural Networks using Genetic Algorithms with Graph Generation System. In *Complex Systems*, volume 4, pages 461–476, San Francisco, CA, USA, 1990. Morgan Kaufmann Publishers Inc.

12. Song Han, Huizi Mao, and William J. Dally. Deep compression: Compressing deep neural networks with pruning, trained quantization and Huffman coding. *4th International Conference on Learning Representations, ICLR 2016 - Conference Track Proceedings*, 2016.

13. Shreyas Saxena and Jakob Verbeek. Convolutional neural fabrics. In D D Lee, M Sugiyama, U V Luxburg, I Guyon, and R Garnett, editors, *Advances in Neural Information Processing Systems*, pages 4060–4068. Curran Associates, Inc., 2016.

14. Jasper Snoek, Hugo Larochelle, and Ryan P. Adams. Practical Bayesian optimization of machine learning algorithms. In *Advances in Neural Information Processing Systems*, volume 4, pages 2951–2959, 2012.

15. Hieu Pham, Melody Y. Guan, Barret Zoph, Quoc V. Le, and Jeff Dean. Efficient neural architecture search via parameter sharing, 2018.

16. Hanxiao Liu, Karen Simonyan, and Yiming Yang. Darts: Differentiable architecture search, 2019.

17. Haifeng Jin, Qingquan Song, and Xia Hu. Auto-keras: An efficient neural architecture search system, 2019.

18. Colin White, Willie Neiswanger, and Yash Savani. Bananas: Bayesian optimization with neural architectures for neural architecture search, 2020.

19. Changyong Oh, Jakub M. Tomczak, Efstratios Gavves, and Max Welling. Combinatorial bayesian optimization using the graph cartesian product, 2019.

20. Aaron Klein, Stefan Falkner, Simon Bartels, Philipp Hennig, and Frank Hutter. Fast bayesian optimization of machine learning hyperparameters on large datasets, 2017.

21. Valerio Perrone, Huibin Shen, Matthias Seeger, Cedric Archambeau, and Rodolphe Jenatton. Learning search spaces for bayesian optimization: Another view of hyperparameter transfer learning, 2019.
22. Rodolphe Jenatton, Cedric Archambeau, Javier Gonzalez, and Matthias Seeger. Bayesian optimization with tree-structured dependencies. In *Proceedings of the 34th International Conference on Machine Learning - Volume 70*, ICML'17, page 1655-1664. JMLR.org, 2017.

23. Xuan Rao, Songyi Xiao, Jiaxin Li, Qiuye Wu, Bo Zhao, and Derong Liu. Lsbo-nas: Latent space bayesian optimization for neural architecture search. In *2022 4th International Conference on Control and Robotics (ICCR)*, pages 22-27, 2022.

24. Dang Nguyen, Sunil Gupta, Santu Rana, Alistair Shilton, and Svetha Venkatesh. Bayesian optimization for categorical and category-specific continuous inputs, 2019.

25. Zi Wang, Clement Gehring, Pushmeet Kohli, and Stefanie Jegelka. Batched large-scale bayesian optimization in high-dimensional spaces, 2018.

26. Stefan Falkner, Aaron Klein, and Frank Hutter. Bohb: Robust and efficient hyperparameter optimization at scale, 2018.

27. Kirthivasan Kandasamy, Willie Neiswanger, Jeff Schneider, Barnabas Poczos, and Eric Xing. Neural architecture search with bayesian optimisation and optimal transport, 2019.

28. Binxin Ru, Xingchen Wan, Xiaowen Dong, and Michael Osborne. Interpretable neural architecture search via bayesian optimisation with weisfeiler-lehman kernels, 2021.

29. Tao Lei, Wengong Jin, Regina Barzilay, and Tommi Jaakkola. Deriving neural architectures from sequence and graph kernels. In Doina Precup and Yee Whye Teh, editors, *Proceedings of the 34th International Conference on Machine Learning*, volume 70 of *Proceedings of Machine Learning Research*, pages 2024-2033. PMLR, 06-11 Aug 2017.

30. Vu Nguyen, Tam Le, Makoto Yamada, and Michael A Osborne. Optimal transport kernels for sequential and parallel neural architecture search, 2021.

31. Klaus Greff, Rupesh K. Srivastava, Jan Koutník, Bas R. Steunebrink, and Jürgen Schmidhuber. Lstm: A search space odyssey. *IEEE Transactions on Neural Networks and Learning Systems*, 28(10):2222-2232, 2017.

32. Alex Graves. Generating sequences with recurrent neural networks, 2013.

33. J. Bergstra, D. Yamins, and D. D. Cox. Making a science of model search: Hyperparameter optimization in hundreds of dimensions for vision architectures. *30th International Conference on Machine Learning*, ICML 2013, (PART 1):115-123, 2013.

34. Frank Hutter, Lars Kotthoff, and Joaquin Vanschoren, editors. *Automatic Machine Learning: Methods, Systems, Challenges*. Springer, 2019.

35. Donald R. Jones. A taxonomy of global optimization methods based on response surfaces. *J. Glob. Optim.*, 21(4):345-383, 2001.

36. Peter I. Frazier. A tutorial on bayesian optimization, 2018.

37. Carl Edward Rasmussen and Christopher K. I. Williams. *Gaussian processes for machine learning*. Adaptive computation and machine learning. MIT Press, 2006.

38. Radford M Neal. *Bayesian learning for neural networks*, volume 118. Springer Science & Business Media, 2012.

39. Donald R. Jones, Matthias Schonlau, and William J. Welch. Efficient global optimization of expensive black-box functions. *Journal of Global optimization*, 13(4):455-492, 1998.

40. Orestis Friderikos, Emmanuel Baranger, and Damien Guillon. Regularized least squares for the building of a sizing criterion based on damage mechanics. *Composite Structures*, 234(October 2019):111653, 2020.

41. Arturo Mendoza, Orestis Friderikos, Roger Trullo, and Emmanuel Baranger. Design of artificial neural networks for damage estimation of composite laminates: Application to delamination failures in ply drops. *Composite Structures*, 304:116320, 2023.

42. M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational Physics*, 378:686-707, 2019.

43. C. Truesdell and W. Noll. *The Non-Linear Field Theories of Mechanics*, pages 1-541. Springer Berlin Heidelberg, Berlin, Heidelberg, 1965.

44. Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. Dropout: A simple way to prevent neural networks from overfitting. *Journal of Machine Learning Research*, 15(56):1929-1958, 2014.

45. Arber Zela, Aaron Klein, Stefan Falkner, and Frank Hutter. Towards Automated Deep Learning: Efficient Joint Neural Architecture and Hyperparameter Search. *arXiv preprint arXiv:1807.06906*, 2018.

46. Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho, and Yoshua Bengio. How to construct deep recurrent neural networks, 2013.
47. Alex Graves, Abdel-rahman Mohamed, and Geoffrey Hinton. Speech recognition with deep recurrent neural networks, 2013.

48. H.T. Siegelmann and E.D. Sontag. On the computational power of neural nets. *Journal of Computer and System Sciences*, 50(1):132–150, 1995.

49. Sho Sonoda and Noboru Murata. Neural network with unbounded activation functions is universal approximator. *Applied and Computational Harmonic Analysis*, 43(2):233–268, 2017.

50. Harold Hotelling. Analysis of a complex of statistical variables into principal components. *Journal of Educational Psychology*, 24:498–520, 1933.

51. Liam Li and Ameet Talwalkar. Random search and reproducibility for neural architecture search, 2019.

52. Kaicheng Yu, Christian Sciuto, Martin Jaggi, Claudiu Musat, and Mathieu Salzmann. Evaluating the search phase of neural architecture search, 2019.

53. Russel E. Caflisch, William J. Morokoff, and Art B. Owen. Valuation of mortgage-backed securities using brownian bridges to reduce effective dimension. *Journal of Computational Finance*, 1:27–46, 1997.

54. Dassault Systèmes. Abaqus Version 6.14.2. *User's Manual*, 2014.

55. K He, S.V Hoa, and R Ganesan. The study of tapered laminated composite structures: a review. *Composites Science and Technology*, 60(14):2643–2657, 2000.

56. M. L. Benzeggagh and M. Kenane. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. *Composites Science and Technology*, 56(4):439–449, 1996.