## 2 Materials and Methods
In addition, there is a broad spectrum of weight and bias optimizers, among which stochastic gradient descent (SGD), SGD with momentum (e.g., NADAM--Nesterov adaptive momentum), and limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithms are one of the most common and well-known in the literature [32]. This complexity may seem as an obstacle; however, it provides a strong, flexible tool to adjust the network to achieve its highest efficiency. The main disadvantage, in this case, is the time-consuming process of finding the optimal network architecture and its set of properties. Once the artificial neural network is trained on a given set of data, the estimation of the objective function using unknown vectors is carried out using the optimized weights and biases from the training procedure. The implementation of ANN and FEA in problems similar to the one presented in this paper was described in [33], where the critical loads for the structures were investigated. The authors created the ANN metamodel using the python programming language (PyTorch) framework, which enables efficient parallelization, assisting in the reduction in model training time significantly. The second metamodel is called the Kriging, after [PERSON], who developed this interpolation metamodel for geostatistics as his master thesis [34]. Although it was designed for mining purposes, such as gold exploration using sampling in a two-dimensional space (map), the method may be adapted for spaces with more dimensions. Kriging models are widely used as surrogate models for costly FEM-based simulations, even when severe nonlinearities are present [35; 36]. There are several types of Kriging models, including simple, ordinary, universal, and blind [37], with the first being used in this paper for its generality and ease of implementation (sci-kit learn Python framework). The general form of the Kriging is the following:  
\[Y(x)=f(x)+Z(x), \tag{7}\]  
where:  
\(Y(x)\)--predicted value;  
\(f(x)\)--trend function;  
\(Z(x)\)--realization of the normally distributed Gaussian random process. Most commonly, the Gaussian correlation function is incorporated in the elements of the covariance matrix (defining the \(Z(x)\)):  
\[\psi^{(i)}=e^{-\sum_{j=1}^{k}\theta_{j}\left|x_{j}^{(i)}-x_{j}\right|^{2}}, \tag{8}\]  
where \(\theta_{j}\) is theta parameter, scaling the width of the basis function, and k is the number of design variables (dimensions). Once the training data are shown to the Kriging model, the hyperparameters of the basis functions are optimized by maximizing the ln-likelihood function. The predictions on the unseen data are calculated using both the training data set and the new, unseen data. Based on the previous studies, the authors decided to focus on a combination of two algorithms, both being a part of the evolutionary optimization group, the genetic algorithm in its multi-island version (MIGA), and the evolution strategy (EVOL). This mixture of algorithms is not unique, as other tailor-made optimization schemes utilize a similar strategy [38]. Gradient methods were not included in this study, as they are not as efficient in optimizing systems, the response of which can be seen as quantified, i.e., the objective function is not continuous in the whole parameter domain, but it has constant values for wide ranges of parameter values, making the gradient calculation incorrect, as presented in Figure 2:0  
The multi-island version of the genetic algorithm is an extension of the work carried out by [PERSON] in 1975. Compared to the basic version of the algorithm, the MIGA consists of several autonomous populations, referred to as islands [39]. Genetic operations, such as selection, cross-over, and mutation are performed (with a certain probability) only among individuals from a single island. Once a fixed number of generations is achieved (migration interval), migration occurs with a certain percentage amount of individuals from each island (migration rate).