Multi-Scale Surface Roughness Optimization Through Genetic Algorithms

Artificial intelligence is changing perspectives of industries about manufacturing of components, introducing emerging techniques such as additive manufacturing technologies. These techniques can be exploited to manufacture not only precision mechanical components, but also interfaces. In this context, we investigate the use of artificial intelligence and in particular genetic algorithms to identify optimal multi-scale roughness features to design prototype surfaces achieving a target contact mechanics response. Exploiting an analogy with biology, the features of roughness at a given length scale are described through model profiles named chromosomes. In the present work, the mathematical description of chromosomes is firstly provided, then three genetic algorithms are proposed to superimpose and combine them in order to identify optimal roughness features. The three methods are compared, discussing the topological and spectral features of roughness obtained in each case.


INTRODUCTION
It is well-established in the literature that surface topology and texture are important for enhancing the tribological behavior of contacts. Therefore, optimization of topological features is considered as a research topic of paramount interest for industrial applications. For instance, optimization of the shape and position of periodic grooved textures has been investigated in Buscaglia et al. (2007) in relation to lubrication. Even more ambitious are the recent attempts to tailor roughness of engineering surfaces by controlling end-milling operations, as discussed in Zhang et al. (2007), or by cutting techniques as in Nouari et al. (2018). Artificial intelligence based on genetic algorithms (Zain et al., 2008) and artificial neural networks (Moghri et al., 2014) have been recently exploited to control milling operations and surface roughness manufacturing by material removal. As another strategy, additive manufacturing (Brettel et al., 2014) is opening new perspectives to produce surfaces with specified roughness (see e.g., Farina et al., 2016;Ko et al., 2019;Wüst et al., 2020).
In addition to tribology, the role of an interface between material constituents/phases is becoming progressively dominant over the one of bulk properties, thanks to the increasing trend in the miniaturization of components, and to the significant progress in the design of mechanical systems and materials, starting from the sub-micrometer scale. The interface has very different properties than the bulk, and is important in a variety of mechanical, physical, and chemical processes. Interfaces can have significant effects also on the properties of composite materials (Moroni et al., 2013;Guo et al., 2014Guo et al., , 2017Kwon et al., 2019). Consequently, industrial applications often require the design of surface textures able to achieve given target responses and/or to enable rapid manufacturing and morphology changes for in-line control of mechanical components, as in Bora et al. (2005) for Micro-Electro-Mechanical Systems (MEMS). Nature is also offering interesting perspectives for the identification of optimal surface topologies for selected applications, see e.g., the successful attempts to create artificial super-adhesives like Gecko's pads (Sherge and Gorb, 2001;Gao et al., 2005;Nosonovsky and Bushan, 2008).
In the present study, a theoretical multi-scale description of roughness is considered in relation to the Multivariate Weierstrass-Mandelbrot (MWM) function as a model for a fractal rough profile (Cinat et al., 2019). The aim of the research is to establish a mathematical and a computational framework for the robust identification of the profile model parameters that lead to identified topologies whose contact response is matching any requested real contact area-load and contact conductance-load curves that can be specified by the user. Although the authors are aware that not all the natural and manufactured surfaces do obey to the fractal scaling, as pointed out by Carr and Benzer (1991), Wen and Sinding-Larsen (1997), Borri and Paggi (2015), Xiaohan et al. (2017), and Gujrati et al. (2018), the choice here is motivated by the fact that the MWM function has been widely studied in relation to contact mechanics, and, for instance, we expect to find that optimization in the low frequency range is the major concern for controlling the thermal/electric contact conductance, while fine scale features are governing the real area of contact (see e.g., Ciavarella et al., 2000Ciavarella et al., , 2004Paggi and Barber, 2011). However, it has to be highlighted that, although the MWM model is herein adopted for validating identification predictions in relation to benchmark results established in the literature, the overall proposed computational procedure can be applied also to other profile or surface models as well, without specific restrictions.
Motivated by an analogy with biology, the expression surface roughness genome is herein used to denote the ensemble of parameters, the genes, associated with the set of elementary waves that describes the main features of roughness over multiple length scales. Although such a terminology imported from biology is not usually adopted in the contact mechanics community, it is introduced here for its common use in the research area of genetic algorithms. Within this framework, we optimize the parameters above to identify a suitable genome that produces a rough profile having a frictionless elastic normal contact mechanics response close to a requested target one. To achieve this goal, we propose various genetic algorithms, combining mechanical considerations and suitable optimization tools. The contact problem is solved in terms of the boundary element method (BEM) (Vakis et al., 2018;Paggi and Hills, 2020), based on the formulation validated in Bemporad and Paggi (2015). See also Bemporad and Paggi (2015) for a wide comparison of other BEM techniques that could also be effectively applied to solve the same contact problem.
Regarding the genetic algorithm research contribution, differently from the previous publication by Cinat et al. (2019) where a single algorithm was proposed to achieve that task, here we propose an extensive numerical comparison of the performance of various algorithms with clear distinct features.
This article is structured as follows. In section 2, we describe the normal contact problem, the surface roughness genome, the features of a single length scale of roughness, and the roughness reconstruction over multiple length scales. In section 3, we propose three genetic algorithms, whose aim is to generate prototype profiles able to achieve given target contact mechanics responses. In section 4, the effectiveness of these algorithms is tested and discussed through a representative example, based on an artificial genome database. Finally, section 5 concludes the paper and discusses possible future work.

NORMAL CONTACT PROBLEM AND ROUGHNESS DESCRIPTION
According to Johnson (2003), the non-conforming contact problem between two rough surfaces is proved to be equivalentunder the assumption of linear elasticity-to the contact problem between a rigid rough surface and an elastic half-plane, provided with suitable effective elastic parameters. The reader is referred to Barber (2003) for details and for a proof of this equivalence. In this work, normal contact is controlled by imposing an approaching far-field closing displacement , which corresponds to a rigid-body motion of the half-plane. Its value is computed from the tallest summit of the rough surface. Elastic interactions among asperities deform the half-plane. During the deformation process, the tallest summit of the profile remains in contact with the half-plane, while other asperities may loose such contact. Their deformation produces a normal contact traction distribution in the contact area A. Then, the average total contact pressure p is computed as the ratio between the sum of all the contact forces acting in this area, and the nominal contact area. Finally, the contact stiffness K is computed as the derivative of the contact pressure p with respect to the imposed displacement . The reader is referred to Bemporad and Paggi (2015) for more details about the precise problem formulation and for a comparison of several techniques that can be exploited to make the unilateral contact constraints satisfied.
In the article, a surface roughness description is provided, based on several parameters. For each choice of such parameters, the contact problem above is solved by an application of the boundary element method (BEM), as in Bemporad and Paggi (2015).

Description of Multi-Scale Surface Roughness
In the following, one-dimensional rough profiles are considered. Likewise in Cinat et al. (2019), their height Z(x), expressed as a function of the position x, is described in terms of the real part of the Multivariate Weierstrass-Mandelbrot function (MWM), which was proposed in Ausloos and Berman (1985) to describe stochastic processes having a larger number of features than in the original 2D formulation proposed earlier in Mandelbrot (1977). In the work, a rough profile Z(x) is described parametrically as Each co-sinusoid in Equation (1) is defined by a unique combination of the parameters H, A, λ, γ > 0, and φ m,n ∈ [0, 2π), whereas M is a positive integer, which represents the number of ridges. In this framework, in analogy with biology, these parameters are also called genes. Thus, the surface roughness genome is the overall ensemble of genes providing the realization of a surface over multiple length scales.

Description of a Generic Length Scale of Roughness
The multi-scale description of profiles is governed by the genes H, λ, γ . The gene γ rules the distance of consecutive wavelengths in the frequency domain, and H their scaling in amplitude. The gene λ indicates a reference wavelength. The particular combination of genes associated with a fixed n identifies a rough profile C n (x), which is named chromosome, and is associated to the features of roughness at the fixed reference length-scale λ n = λγ 1−n . In biology, a chromosome is a structure composed by some genes identifying specific features of the genome (King et al., 2006). Here, it identifies the features of roughness at the wavelength λ n , as follows: A profile Z(x) is realized in an observation length L in N = L δ + 1 nodes, where δ is the distance between two consecutive nodes, i.e., the resolution, and L is a multiple of δ. Hence, the infinite series in Equation (1) is replaced by a finite series, according to the observation scale L and the resolution δ chosen. The profile Z(x) is realized from its longest to its shortest observed wavelength by summing a finite number n c of chromosomes: where n s (starting index) and n f (final index) refer to the chromosomes contributing with the longest and shortest wavelengths to the realization of the rough profile. It holds n c = n f − n s + 1. The realization of roughness over multiple length scales of observation is obtained by selecting different chromosomes, i.e., by varying n between n s and n f . The value n s refers to the longest wavelength considered, and n f refers to the shortest wavelength observed. In this way, the chromosomes contributing to a realization at a given length scale depend on n s and n f . Without loss of generality, the value n s = 1 is assigned to the chromosome with longest wavelength referring to the coarsest realization of a surface. In more details, denoting by ⌊·⌋ the largest integer smaller than or equal to its argument, one has To visualize the concept of chromosome, a simple example is presented in Figure 1. Here, two chromosomes are shown for two different values of n (n = 1, 2). In both cases, the parent M = 8 co-sinusoids are shown in Figures 1A,B through the colored lines.
The chromosome C 1 (x) is assumed to have a wavelength equal to the sample wavelength, i.e., λ 1 = L = 100 µm. This chromosome, depicted by a thick red line in Figure 1A, is obtained by summing all the M = 8 co-sinusoids with the associated value n = 1.
The chromosome C 2 (x), depicted by blue line in Figure 1B, is obtained as the sum of all the M = 8 co-sinusoids with n = 2. Then, since it is imposed γ = 1.5, the chromosome C 2 (x) has λ 2 = λ 1 1.5 = 66.6 µm. Each of the two chromosomes in Figure 1 maintains a cosinusoidal shape with the same wavelength λ n as the associated parent co-sinusoids, and with a phase angle θ . Its height field in Equation (2) can be also written as In Equation (5), the amplitude parameter G n = A log 10 (γ ) M 2π λ −H γ −(n−1)H has been introduced, together with the two constants g n,1 and g n,2 , which depend on the phases φ m,n of the ridges with first index m: A feasible mathematical model of a chromosome might be one in a form similar to that of the MWM function in Equation (1) with M = 1 and a single value of n, obtained by introducing the parameters K n , θ n,1 , and θ n,2 : A particular case of Equation (7) is when the angles θ n,1 and θ n,2 coincide, say, with the same θ n . In such a case, a chromosome is described exactly according to the MWM profile in Equation (1) with M = 1 and a single value of n. The expressions of K n , θ n in Equation (7) are obtained equating Equations (5) and (7), and are given by where the parameter g n,3 is defined as follows: For its validity, Equation (8) requires the argument of arccos g n,1 √ M+2g n,3 to be between −1 and 1, i.e., g n,1 ≤ M + 2g n,3 .
In the above, the additional condition M+2g n,3 ≥ 0 has not been imposed explicitly, since it always holds. Indeed, one can easily check that M + 2g n,3 is the square of the Euclidean norm of the vector with components ( norm is always larger than or equal to 0. A numerical simulation has been conducted to verify the validity of condition (10). This has been evaluated for several values of M, considering each time 1,000 random choices (uniformly and independently sampled between 0 and 2π) for the phases. Figure 2 reports, for each such value of M, the percentage of cases for which Equation (10) is satisfied. Such percentage ranges between 88.5 and 98%, with an average value of ∼ 95%.

Roughness Reconstruction Over Multiple Length Scales
The definition of chromosome allows the reconstruction of profiles Z(x) over multiple approximate realizations of the same surface. This can be done by limiting the summation in Equation (3) to a subset of chromosomes C n (x).
For both the stiffness-load curve K(p) and the area-load curve A(p), only one subset of the chromosomes can be considered to represent the main features of the contact mechanics response. This is shown in Figure 3 for a representative example coming from a fractured alloy interface (black line), considering both the K(p) and A(p) evolutions. The red dashed line is the contact mechanics response of the profile obtained by summing up chromosomes that have a K(p) evolution with a correlation coefficient larger than 0.95 with the complete K(p) one. Both the K(p) and A(p) evolutions of this profile overlap significantly with the ones corresponding to the complete profile, which is obtained when all the chromosomes are taken into account in the summation. This means that the interaction of chromosomes composing only one part of the complete genome is sufficient to approximate the stiffness-load curve K(p) well, as previously deduced in Paggi and Barber (2011).
The selection of chromosomes outlined above is detailed in Algorithm 1. From the genome of a surface, one obtains the K(p) evolution of the profile described by Equation (1) (Step for the alloy fractured profile in the one-dimensional case. The red line refers to the contact mechanics response of the profile obtained summing up all the chromosomes with a correlation coefficient c n larger than 0.95 with the complete K(p) curve. In the legend, "seq" stands for "sequenced," and refers to the complete profile (the one identified from the whole set of genes). 1). The mechanical curves K n (p) of all the n c chromosomes composing the considered rough profile are iteratively extracted using the BEM algorithm for the frictionless elastic normal contact problem, which is solved via the Non-Negative Least Squares method proposed in Bemporad and Paggi (2015) (Steps 2-7). The correlation coefficients c n between the extracted curves and the one K(p) of the complete rough profile are calculated (see Step 6). Only chromosomes with a correlation coefficient higher than 0.95 are then retained in the final set U c of chromosomes (see Step 8).

GENETIC ALGORITHMS TO IDENTIFY OPTIMAL PROFILES TO MATCH TARGET CONTACT MECHANICS RESPONSES
The goal of this section is to illustrate methods aimed to design a prototype profile able to achieve a target contact mechanics response y t (ξ ) of a rough profile, where t stands for "target". This contact mechanics response depends on the specific needs of the frictionless elastic normal contact problem and is discretized using n t points, where n t is the number of far-field displacements imposed to the frictionless elastic normal contact problem (see Johnson, 2003;Bemporad and Paggi, 2015). This target evolution y t (ξ ) can be, e.g., either the stiffness-load curve K(p) or the contact area-displacement curve A( ). The unknown rough profile Z t (x) to be identified is discretized using N = 512 nodes for a length L. Then, the values of n s and n f are defined according to Equation (4). Here, without loss of generality and for simplicity, it is imposed n s = 1 for each realization. The frictionless elastic normal contact problem is solved in n t = 20 equi-spaced rigid body displacements from the topmost summit of the considered profile to its deepest valley.
In the following, the variable ξ is taken as the contact pressure p. Moreover, the gene A i of a generic genome (indexed by i in Algorithm 1 | Chromosome selection Input: genome of a profile, n s , n f , number n t of far-field displacements to be imposed in BEM Output: U c : set of chromosomes determining the K(p) evolution of the profile 1: y r ← BEM results (K(p)) for the reference profile given by the genome, n s , and n f 2: for all n = n s : n f do c n ← corr. coeff.(y r , y n ) 7: end for 8: U c ← C n (x) with c n > 0.95 the database) is rescaled to match a given pressure requirement, p ≤ p t max , as follows: The maximum pressure p i max is computed imposing a farfield displacement equal to the original profile amplitude. In such a way, the new profile shows a maximum pressure p t max if a far-field displacement equal to the peak-valley amplitude is considered to solve the frictionless elastic normal contact problem.
Three different methods to design the prototype profile are proposed. All these methods start from a known database of genomes. In the first case, a genome is selected that has the closest contact mechanics response to y t among all the genomes in the database. Its genes are then optimized to increase the similarity with the target contact mechanics response. In the second case, two genomes are selected that match the target response in different imposed ranges p. Then, these genomes are combined using an optimized cross-over mechanism. The third and last case is similar to the second one but, instead of the complete genomes, two sets of chromosomes are mixed. These two sets are contained in two genomes that match the target response in different imposed ranges p.

Simple Optimization of Genes
In this section, the Simple Optimization of Genes (SOG) genetic algorithm is presented. Genes to be optimized belong to a genome producing a rough profile with a contact mechanics response very similar to the target one. The following similarity score is defined to quantify how much the target contact mechanics response, y t (ξ ), and the one y i (ξ ) associated with the parent i-th rough profile, are similar. Here, · ∞ denotes the l ∞norm, computed on the n t points used to discretize the contact mechanics response. Due to Equation (12), the i-th curve coincides with the target one when s i = 1 holds. Otherwise, it is close to the target one when s i ≃ 1 holds. The numerical steps performed to select the best profile from the set of N D genomes in the given database are described in Algorithm 2. The similarity score in Equation (12) is computed for all the N D rough profiles available in the database, rescaled according to Equation (11) (Step 2). The contact mechanics response y i (ξ ) is computed via BEM, considering n t different farfield displacements, from 0 to the new profile amplitude (Steps 3-5). Finally, s i = s(y t , y i ) in Equation (12) is computed (Step 6).
Algorithm 2 | Similarity score extraction from a database of genomes Input: target contact mechanics response y t , database of N D genomes (genes, p i max , n t ) Output: s i (y t , y i ) with i ∈ {1, . . . , N D } Genes belonging to the three genomes with the three largest associated values of s(y t , y i ) are now optimized using the Globally Convergent Method of Moving Asymptotes (GCMMA) (Svanberg, 2002). This is an iterative optimization algorithm, which is often used in optimal design for mechanical problems (see Bacigalupo et al., 2016Bacigalupo et al., , 2017 for some of its recent applications to band gap optimization). The GCMMA algorithm replaces a nonlinearly constrained optimization problem by a sequence of approximating nonlinearly constrained optimization subproblems, which are simpler to solve. In the work, the objective function has been chosen to be the square of the similarity score, i.e., s 2 (y t , y i ), in order to increase its smoothness. The GCMMA iterative solution is obtained after a number n it of steps, starting from an initial choice for the vector of optimization variables. In the case of the SOG, this initial choice is provided by Algorithm 2, which provides an initial value of the similarity score close to 1. In such a way, there is no risk to get a negative similarity score by locally maximizing its square.
To save CPU time, the number of variables in the optimization problem is limited as follows. The genes H, λ, and γ determine the frequency spectrum, i.e., the interaction among different length scales. Therefore, they are not varied during the optimization, but fixed to their original values. The gene A is varied at each optimization step to satisfy the pressure requirement, according to Equation (11). So, only the phases φ m,n are considered as optimization variables. In the optimization, their values are constrained in the range between ∓10% of their initial values, to preserve the main features of the original chromosomes. Moreover, only the phase genes φ m,n of chromosomes that determine the main features of the K(p) response are considered. Such genes are selected according to Algorithm 1.
Finally, the SOG algorithm is summarized in Algorithm 3. Starting with a profile scouting from the available database (Step 1), the three profiles whose contact mechanics responses are most similar to the target y t are chosen (Step 2). The genes of each such genome are then optimized using the GCMMA algorithm (Steps 3-6), limiting the optimization variables only to the chromosomes determining the main features of the target y t , as determined by Algorithm 2. The resulting optimized genomes are denoted byÛ i 1 . Finally, among such genomes, the new genome is chosen as the one with the best (square of the) similarity score with respect to the target response (Step 7). In this last step, argmax(f i ) i=1,...,n 1 denotes the index i associated with the largest f i .

Algorithm 3 | Simple Optimization of Genes (SOG)
Input: target contact mechanics response y t , genome database Output: new genome U SOG with contact mechanics response close to y t 1: s from Alg. 2 2: U 1 ← the three genomes with the three largest similarity scores (Eq. 12) 3: for all i = 1 : n 1 do (n 1 = card(U 1 )) 4: U i c ← Alg. 1 applied to U i 1 5: f i ← s(y t , y i ),Û i 1 , both from GCMAA(U i c ) 6: end for 7: U SOG ←Û argmax(f i ) i=1,...,n 1 1

Genome Cross-Over
In this section, the Genome Cross-Over (GCO) genetic algorithm is presented. In this method, different pairs of genomes are mixed, to obtain a new genome matching the target response y t (ξ ). The best pair of genomes is chosen in relation to their similarity scores in two specific ranges of the target response y t (ξ ).
This concept is explained through Figure 4, where a target response y t (p) is shown by the red line. For the sake of explanation, two curves y 1 (p) (dashed black line) and y 2 (p) (dashed dot blue line) have been manually generated, which have a similarity score with respect to y t (p) equal to s 1 ≃ 0.89 and s 2 ≃ 0.88, respectively. However, the curves y 1 (p) and y 2 (p) describe quite accurately the curve y t (p) in different ranges of pressures, if a suitable threshold pressurep is defined. The value ofp can be chosen arbitrarily, or might be imposed by the particular problem formulation. In the specific case shown in the figure, in the interval [0,p], the curve y 1 represents with good accuracy y t [s 1 (y t , y 1 ) ≃ 0.97]. The same happens in the interval [p, p t max ] for the curve y 2 [s 2 (y t , y 1 ) ≃ 0.99].
Then, it is reasonable to expect that a new profile obtained by combining the two genomes associated with y 1 and y 2 , respectively, should provide a contact mechanics response closer to y t over the whole range of pressures. However, mixing these two genomes may also lead to a very different roughness organization. For this reason, the GCO iterative scheme checks if the new genome, obtained by mixing the genome of y 1 and the one of y 2 , is able to represent accurately y t over the whole range of pressures, before using the GCMMA algorithm to increase the value of the similarity score.
The GCO structure is presented in Algorithm 4. According to the value ofp, two different sets U 1 and U 2 are identified from the database of genomes (Steps 1-2). The first set U 1 contains the genomes with a value of the similarity score (Equation 12) larger than 0.95 in the interval [0,p] (Step 3). The second set U 2 contains the genomes with a value of the similarity score (Equation 12) larger than 0.95 in the interval [p, p t max ] (Step 4). All possible combinations of genomes from the two sets above are now considered, defining the set U 3 (Step 5-11). A new genome corresponds to each of these combinations, whose gene A is rescaled according to Equation (11) (Step 8). Then, the value of the similarity score (Equation 12) is computed with respect to the target response (Step 9). The three new genomes showing the three largest values of the similarity score are used as inputs to the GCMMA algorithm, defining the set U 4 (Step 12). Also in this case, the number of genes to be optimized is limited, considering only the phases φ m,n associated with the chromosomes that determine the main features of the K(p) evolution of the parent curve (see Algorithm 1). The new genome U GCO with the maximum value of the similarity is finally identified (Step 17).

Chromosomes Cross-Over
In this section, the Chromosomes Cross-Over (CCO) genetic algorithm is presented. We recall that, as observed in section 2, the main features of the contact mechanics response of a rough profile are determined by specific chromosomes of the genomes. This subdivision allows the introduction of a chromosomes selection step to reduce the number of variables to be optimized using the GCMMA algorithm, as presented before for the SOG (section 3.1) and the GCO (section 3.2).
In the CCO method, only chromosomes determining the main features of the contact mechanics responses of two different Algorithm 4 | Genome Cross-Over (GCO) Input: target contact mechanics response y t , genome database, threshold pressurep Output: new genome U GCO with contact mechanics response close to y t 1: s (1) from Alg. 2, with similarity score computed in the interval [0,p] 2: s (2) from Alg. 2, with similarity score computed in the interval [p, p t max ] 3: U 1 ← genomes with s (1) i > 0.95 4: U 2 ← genomes with s (2) i > 0.95 5: for all i 1 = 1 : n 1 do (n 1 = card(U 1 )) 6: for all i 2 = 1 : n 2 do (n 2 = card(U 2 )) 7: rescaled according to Eq. (11) 9: s(i 1 , i 2 ) ← s(y t , y (i 1 ,i 2 ) ) from Eq. (12) applied to U (i 1 ,i 2 ) 3 10: genomes are mixed to match the target response y t . These two sets of chromosomes come from genomes that have the largest values of the similarity score in specific ranges of the target response y t , as done for the GCO method (see section 3.2). Then, in this case the GCMMA is applied to the complete set of chromosomes composing this new genomes, since they all determine its mechanical evolution.
The CCO iterative scheme is summarized in Algorithm 5. According to the value of the threshold pressurep, two different sets U 1 and U 2 of reduced genomes are identified, starting from the given database (Steps 1-2). The first set U 1 is obtained as follows. First, one generates a subset of genomes that show a value of the similarity score (Equation 12) larger than 0.95 in the interval [0,p] (Step 3). Then, these genomes are reduced, limiting to those chromosomes that affect the K(p) evolution significantly. Such chromosomes are selected according to Algorithm 1. The resulting reduced genomes form the set U 1 . The second set U 2 is obtained in a similar way as U 1 , but computing the similarity score in the interval [p, p t max ] (Step 4). All possible combinations of these reduced genomes from the two sets above are now considered, defining the set U 3 (Steps 5-11). A new genome corresponds to each of these combinations. Its amplitude gene A is rescaled according to Equation (11), to match the pressure requirement (Step 8). Then, the value of the similarity score (Equation 12) is computed with respect to the target response (Step 9). The three new genomes showing the three largest values of the similarity score are used as inputs to the GCMMA algorithm, defining the set U 4 (Step 12). Only in the case of the CCO algorithm, the GCMMA algorithm is applied to all the genes of this new genome, as the size of the optimization problem has been already reduced in Steps 3-4. The new genome U CCO with the maximum obtained value of the similarity score is finally identified (Step 16).

OPTIMAL GENOME TO MATCH A SPECIFIED CONTACT MECHANICS RESPONSE
In this section, the genetic algorithms described in section 3 are compared in their application to a representative example. Before doing that, the characteristics of the numerical genome database considered in the paper are reported. Then, the setup of the numerical experiment proposed is introduced. Finally, the related results are discussed, focusing the attention on the characterization of the so-obtained new genomes.

Database of Genomes
A database of genomes is needed to apply all the genetic algorithms proposed in the paper. In the following, a small database is generated numerically, to show a representative example of this approach. To generate the database, the amplitude parameter is fixed to A = 1, and the main wavelength is put equal to λ = 849.42 µm.
The number of ridges is set to M = 1, to save computational time. This particular case corresponds to chromosomes described according to Equation (7) with θ n,1 = θ n,2 = θ n . As discussed Algorithm 5 | Chromosomes Cross-Over (CCO) Input: target contact mechanics response y t , genome database, threshold pressurep Output: new genome U CCO with contact mechanics response close to y t 1: s (1) from Alg. 2, with similarity score computed in the interval [0,p] 2: s (2) from Alg. 2, with similarity score computed in the interval [p, p t max ] 3: U 1 ← C n (x) from Alg. 1, for those genomes with s (1) i > 0.95 4: U 2 ← C n (x) from Alg. 1, for those genomes with s (2) i > 0.95 5: for all i 1 = 1 : n 1 do (n 1 = card(U 1 )) 6: for all i 2 = 1 : n 2 do (n 2 = card(U 2 )) 7: s(i 1 , i 2 ) ← s(y t , y (i 1 ,i 2 ) ) from Eq. (12) applied to U (i 1 ,i 2 ) 3 10: end for 11: end for 12: U 4 ← the three genomes in U 3 with the three largest similarity scores (obtained from s) 13: for all i = 1 : n 4 do (n 4 = card(U 4 )) 14: ..,n 4 4 in section 2, this can be considered as representative of several situations for a rough profile. Each profile Z i (x) is then discretized using N = 512 nodes. The elastic modulus is imposed to be equal to E = 1 MPa. Furthermore, n h = 20 different pairs of values for H and γ are considered. They are generated according to a Sobol' sequence (Niederreiter, 1992) (see Figure 5A), with H ranging between 0.5 and 1.5, and γ between 1.4 and 2. The value of γ is chosen in such a way to keep small the number of frequencies composing the profile spectrum (see Equation 4), thus limiting the computational costs. Also the phase matrix is generated according to a Sobol' sequence. This matrix is made of n φ = 5 columns, corresponding to 5 choices for the set of phases. The number of rows is equal to the maximum possible number of phases for these choices of the parameters. According to Equation (4), such number is obtained in correspondence of the smallest value of γ . Then, N D = n h n φ rough profiles Z i (x) are generated, combining each pair (H, γ ) to each column of , considering only the number of values needed for the phases.

Set-Up of the Numerical Experiment
In Figure 6, the contact mechanics response y t = K(p) is depicted through the black line. All the values are considered independent from E, as this is equal to unity. Also the K(p) responses of all the genomes in the generated database are visualized in Figure 6. The target contact mechanics response y t has been generated manually, and has a trend similar to the ones belonging to the database. The maximum pressure  reachable is imposed to be p t max = 1.6 × 10 −4 N/m. This value is larger than the ones of all the genomes in the database considering, for each profile, a far-field displacement equal to the profile amplitude.
For each algorithm, the stage just preceding the application of the GCMMA algorithm is referred to as "GCMMA (0)." For the SOG, it corresponds to the selection of the three best values obtained from the scouting of the database (Step 2 in Algorithm 3). For both the GCO and the CCO algorithms, it corresponds to the computation of the similarity scores of the new genomes determined after the cross-over (Step 12 in both Algorithms 4 and 5). For the SOG, the GCMMA algorithm is applied with n it = 5, while for the GCO and CCO, only with n it = 3. The output of this optimization step is denoted by "GCMMA (1)." These different choices for n it are made in order to have comparable computational times of about 1 minute. In such a way, it is possible to observe variations in the results in a sufficiently small time, enabling in-line control as a possible successive step. Moreover, in the numerical experiment, for the SOG and GCO algorithms, the GCMMA is also applied in a second optimization step, whose output is denoted as "GCMMA (2), " using n it = 2 iterations. In this case, the genes to be optimized are the ones referred to the chromosomes excluded in the selection step (see Algorithm 1). For the CCO, this second step cannot be applied since these chromosomes are excluded in the first part of the algorithm (Steps 1-4 in Algorithm 5). However, to make the notation uniform for the three algorithms, a fictitious "GCMMA (2)" output is introduced for the CCO, by duplicating the results of "GCMMA (1)."

Numerical Results
The values of the similarity score s given by Equation (12) for the genomes obtained by combination and optimization using the SOG, GCO, and CCO algorithms are reported in Figure 7. The threshold pressure used for both the GCO and CCO is set equal top = 0.08 [1/mm]. For each algorithm, solutions at different stages are reported in the figure. The first ones correspond to the three best solutions obtained in the first stage of each algorithm, before the application of the GCMMA (such solutions are represented in the figure by the triangle, circle, and cross, in increasing order). The figure shows that all the algorithms are quite efficient in matching the target contact mechanics response, achieving large values for the similarity score. The application of the GCMMA is beneficial for all the algorithms. However, the FIGURE 7 | For the outputs of the SOG, GCO, and CCO algorithms, values of the similarity scores with respect to the target curve y t . The threshold pressure isp = 0.08 for both the GCO and CCO.
GCO and CCO algorithms might provide even better solutions by varying the threshold pressurep.

Effect of the Threshold Pressurep
A threshold pressurep has been introduced for the GCO and the CCO, to select individual genomes able to approximate locally the target curve with a good accuracy. Additional simulations have been made, to assess the sensitivity of the results with respect to such a parameter. Both the GCO and CCO algorithms have been applied with different values ofp. Such values have been chosen between 0.5 × 10 −4 N/m and 1.1 × 10 −4 N/m. A significant variation on the maximum value of the similarity score is observed for the GCO algorithm (see Figure 8A). On the contrary, the CCO algorithm seems to be not affected by the value ofp, and the new genomes obtained by that algorithm are composed by the same starting genomes, independently of the threshold pressurep. This may be due to the fact that our investigation has been conducted starting from a small database of genomes. In Figure 8B, the cardinality of the set U 3 is presented, for both the GCO and CCO algorithms. This set contains the new genomes obtained after the cross-over of genomes/chromosomes. The white part indicates, for both algorithms, the number of genomes whose value of the similarity score with respect to the target curve y t is larger than 0.95. The cardinality of the set U 3 varies significantly withp, and is the same for both algorithms. Nevertheless, the GCO algorithm is more efficient in terms of the obtained similarity with the target curve.

Description of the Optimized Genomes Representing y t
The best three new genomes obtained at the end of the three algorithms presented in this paper have contact mechanics responses overlapping significantly with the target curve (see Figure 9A). However, the three genomes present different evolutions of the bearing area, as depicted in Figure 9B. Only the A(p) curves obtained by the SOG and GCO are similar. This may be due to the fact that, in the case of the CCO, high-frequency  FIGURE 10 | Topography of the best rough profiles approximating the target curve y t (see Figure 9). features of the rough profile might have been neglected. This is clear from Figures 10, 11.
The topography of these rough profiles is presented in Figure 10. The profiles obtained by the SOG, GCO, and CCO algorithms are depicted, respectively, through a black dash dotline, a red dashed line, and a blue continuous line. All these profiles have very similar geometrical features, regarding the locations of peaks and valleys. Moreover, it is interesting to notice that the profile provided by the CCO algorithm is a good approximation of the profile given by the GCO algorithm, which presents more high-frequency features.
To conclude, the discrete Power Spectral Density (PSD) P(ω) of the new obtained genomes is shown in Figure 11, and is represented by markers in both sub-figures. For each case, the continuous PSD function, which is obtained through the Fast Fourier Transform (FFT) filtering (Berry and Lewis, 1980), is shown by a continuous line. For the SOG (see Figure 11A), the peaks of the continuous PSD function match the discrete one accurately for high frequencies. No good matching is found for low frequencies. The same remark holds for the case of the GCO (see Figure 11B). Here, the spectrum is more dense and, only for high frequencies, peaks of the continuous PSD function are located in the same positions of the genome. For both the SOG and the GCO, a consistent difference in the amplitude of P(ω) is found. Finally, the spectrum of the profile obtained by the CCO is composed by a small set of frequencies. The high-frequency part of the PSD is flat. In this case, the discrete power spectral density is nearly proportional to the one of the SOG in the lowfrequency range. Since the GCO approximation of the target curve is obtained using a smaller set of chromosomes, it can be easily controlled for the case of in-line control of the rough profile morphology. However, in this case, the PSD is almost flat, and some peaks are found in frequencies where there are no genome components.

CONCLUSION
The main goal of this work was to provide a mathematical and computational methodology to identify new surface topologies with given target contact mechanics responses, based on a multi-scale characterization of roughness. This was described by a superposition of a finite set of successive length scales of roughness, where each scale is associated with a particular form of the MWM function, named chromosome. Based on this representation, three genetic algorithms were proposed to combine profiles and identify the best one which meets the given target contact mechanics response.
The first genetic algorithm (Simple Optimization of Genes, SOG) optimizes the genes of a known genome. Embedded in the SOG, an iterative optimization algorithm (the GCMMA) was used to increase the similarity with the target contact mechanics response. To save computational time, optimization was performed only with respect to the dominant chromosomes of the genome, i.e., the ones for which the correlation coefficient between the associated contact mechanics response and the target one is above a given threshold.
The second genetic algorithm (Genome Cross-Over, GCO) crosses-over two different genomes that show a high similarity with respect to the target response in two given intervals determined by a threshold pressurep. Also in this case, only dominant chromosomes were optimized to save computation time.
The third genetic algorithm (Chromosomes Cross-Over, CCO) consists in a cross-over of chromosomes that provide the main features of the contact mechanics responses of two genomes that well approximate the target curve in two reference intervals determined by the threshold pressurep. In this case, GCMMA was applied to the complete new genome, since the number of optimization variables is smaller than for the other two approaches.
The three genetic algorithms proposed in this work generated profiles that almost fully reproduce the given K(p) curve, when this is chosen as the target contact mechanics response. Also, they have in common similar features regarding their topography, such as the locations of peaks and valleys. This characterization might be caused by the fact that a limited number of genomes was present in the database used for the numerical experiment reported in the paper. A much wider set of profile topologies is expected to be obtained by exploiting a significantly larger database, facilitating the discovery of optimal patterns/textures, based on the specific needs. Also, the investigation can be extended to all the genomes that provide a good approximation of the target curve, in order to identify some geometrical features that drive the elastic response of a rough profile. Once identified, these features can be taken into account in the optimization by adding specific chromosomes, to be controlled in the case of in-line profile morphing.
Finally, the authors would like to remark that the proposed approaches could be extended to the two-dimensional case, although this is expected to require a larger computational effort. Moreover, the use of micromechanical contact theories such as those compared in Zavarise et al. (2007) and Paggi and Ciavarella (2010) could be another interesting research direction. For instance, one could assume not to have the surface height field at all, but only a set of relevant statistical parameter inputs of micromechanical contact theories. In such a context, the search for the optimal topology to match a prescribed contact response would reduce to the identification of the statistical parameters based on the combination of the initial dataset and of the proposed genetic algorithms.

DATA AVAILABILITY STATEMENT
The dataset generated for this study is available on request to paolo.cinat@alumni.imtlucca.it.