## 1 Introduction  
In inverse engineering applications, causal factors that govern the response of a currently analyzed system must be determined. For example, one of the most common inverse problems in mechanics is to identify material properties (e.g., Young's modulus, viscoelastic or hyperelastic constants, etc.) and/or stiffness properties (e.g., tension/bending/shear terms in shells, fiber orientations and ply thicknesses in composite laminates, etc.) corresponding to the given displacements {\(u\)(\(x\),\(y\),\(z\)}, \(v\)(\(x\),\(y\),\(z\)), \(w\)(\(x\),\(y\),\(z\))}. The most popular approaches adopted in material identification are the finite element model updating method (FEMU) [1; 2; 3; 4; 5] and the virtual fields method (VFM) [3; 4; 5; 6]. FEMU has a higher computational cost, while VFM can be easily implemented. Measured displacements are compared at \(N_{\text{CNT}}\) selected control points with the corresponding displacement values computed by an FE model simulating the experiment. Measured and computed displacements match if actual materials/structural properties are inputted to the FE model: the difference is the error functional \(\Omega\), which depends on the \(NDV\) properties to be identified. The inverse problem is hence stated as an optimization problem:\[\left\{ \begin{array}{l}Min\left[\Omega(X_{1},X_{2},\ldots,X_{\text{NDV}})=\sqrt{ \frac{1}{N_{\text{CNT}}}}\sum\limits_{k=i}^{\text{NCNT}}\left(\frac{g_{\text{ FDA}i}-\overline{\delta^{i}}}{\delta^{i}}\right)^{2}\right]\\ \text{G}_{\text{P}}(X_{1},X_{2},\ldots,X_{\text{NDV}})\geq 0\\ \chi_{i}^{\text{L}}\leq X_{i}\leq X_{i}^{\text{U}}\end{array} \right. \tag{1}\]  
In Equation (1), the design vector \(\mathbf{X}_{\text{PROP}}(X_{1},\,X_{2},\,\ldots,\,X_{\text{NDV}})\) includes the \(NDV\) unknown properties ranging between lower and upper bounds \"L\" and \"U\", respectively; \(\delta_{\text{FEM}}\) and \(\overline{\delta^{i}}\), respectively, denote computed and target measured displacements for the \(i\)th control point; and constraint functions \(G_{\text{P}}(X_{\text{PROP}})\) in the unknown properties guarantee the existence of a solution for the inverse problem. The formulation of error functional \(\Omega\) depends on the type of identification problem; for example, in indentation/nanoindentation problems, reaction forces at the indenter/body interface also may be considered, while in dynamic problems, natural frequencies and modal shapes can be utilized as target quantities; strains and stresses also may be involved in the definition of \(\Omega\). [PERSON] and his collaborators [7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21] carried out pioneering work on mechanical identification of materials and structures, considering a variety of cases such as, for example, elasto-plasticity (also accounting for material anisotropy, inhomogeneity and functional grading), indentation, residual stresses, multiaxial fatigue, damage and fracture, thin foil mechanics, etc. The inverse problem defined by Equation (1) may be highly nonlinear and nonconvex if there are many unknown material/structural parameters to be determined and each new evaluation of error functional \(\Omega\) requires a new FE analysis. Since the \(\Omega\) functional may be sensitive only to some of the unknown parameters, sensitivity analyses should be carried out to select experimental data most related to these parameters. Such a strategy makes the inverse problem well posed but increases the computational cost of the identification process (see, for example, [10; 15; 18; 22; 23] for detailed discussions on this issue). However, model reduction techniques may be used to speed up computations [24]. The non-smooth optimization problem (1) can be solved using metaheuristic optimization algorithms (MOAs) that randomly generate trial designs according to some principles inspired by physics, evolution, biology, social sciences, human activities, etc.