Artificial intelligence-based predictive model of nanoscale friction using experimental data

A recent systematic experimental characterisation of technological thin films, based on elaborated design of experiments as well as probe calibration and correction procedures, allowed for the first time the determination of nanoscale friction under the concurrent influence of several process parameters, comprising normal forces, sliding velocities, and temperature, thus providing an indication of the intricate correlations induced by their interactions and mutual effects. This created the preconditions to undertake in this work an effort to model friction in the nanometric domain with the goal of overcoming the limitations of currently available models in ascertaining the effects of the physicochemical processes and phenomena involved in nanoscale contacts. Due to the stochastic nature of nanoscale friction and the relatively sparse available experimental data, meta-modelling tools fail, however, at predicting the factual behaviour. Based on the acquired experimental data, data mining, incorporating various state-of-the-art machine learning (ML) numerical regression algorithms, is therefore used. The results of the numerical analyses are assessed on an unseen test dataset via a comparative statistical validation. It is therefore shown that the black box ML methods provide effective predictions of the studied correlations with rather good accuracy levels, but the intrinsic nature of such algorithms prevents their usage in most practical applications. Genetic programming-based artificial intelligence (AI) methods are consequently finally used. Despite the marked complexity of the analysed phenomena and the inherent dispersion of the measurements, the developed AI-based symbolic regression models allow attaining an excellent predictive performance with the respective prediction accuracy, depending on the sample type, between 72% and 91%, allowing also to attain an extremely simple functional description of the multidimensional dependence of nanoscale friction on the studied variable process parameters. An effective tool for nanoscale friction prediction, adaptive control purposes, and further scientific and technological nanotribological analyses is thus obtained.


Introduction
Tribology, with its marked influence on manufacturing processes, energy consumption, environmental impacts, aerospace technologies or, more recently, the microand the nanotechnologies, especially the study of the fundamental physio-chemical aspects determining the origin of friction at the nanoscales, is a propulsive and evolving field of study [1,2]. An innovative structured approach to the experimental assessment of the concurrent influence of the most important process parameters on dry (unlubricated) nanoscale friction in single asperity contacts was introduced in a recent study [3]. This procedure is based on an elaborated design-of-experiments (DoE) methodology, conducted by using centroidal Voronoi tessellation (CVT) sampling [4], as well as on a carefully conceived characterisation of the stiffness of the used scanning probe microscopy (SPM) probes and of the influence of tip wear and adhesion on the obtained results. The measurements were performed by using the lateral force microscopy (LFM) measurement mode, i.e., lateral (transversal) scans performed on 500 nm × 500 nm surfaces of the analysed aluminium oxide (alumina, Al 2 O 3 ), titanium dioxide (TiO 2 ), molybdenum disulphide (MoS 2 ), and aluminium (Al) samples. The value of nanoscale friction F f , in the multidimensional space defined by the used thin film materials, the normal loads acting on the samples (F N = 10, … , 150 nN), sliding velocities (v = 5, … , 500 nm/s), and temperatures ( = 20, … , 80 ℃), could therefore be experimentally obtained in a set of 50 discrete sampling points. By using in a first instance firstorder statistical analyses, based on Pearson's product-moment correlations (PPMC) [5], important insights into the general trends of the dependence of nanoscale friction on the multiple studied process parameters were consequently obtained, confirming that their interactions and mutual effects must be investigated at the structural atomic level to be fully appreciated.
With the goal of overcoming the limitations of currently available models in determining the effects of the processes and phenomena involved in nanoscale contacts, and obtaining true predictive models linking the cited process variables to the value of nanometric friction, the measurement results obtained in Ref. [3] are validated in this work numerically. To obtain the complete characterization of nanoscale friction from relatively sparse available experimental data, stateof-the-art data mining processes, employing various types of black-box machine learning (ML) [6] and function-generating symbolic regression (SR) artificial intelligence (AI) [7] methods, are thus employed. The modelling process relies on training the various models with existing experimental data, while new and independent measurements on unseen sampling points are made to test the developed models' predictive performances. Thorough comparative analyses are then carried on. Despite the marked complexity of the analysed phenomena and the inherent dispersion of the measurements, the developed symbolic regression model allows achieving an unprecedented prediction accuracy. The obtained simple functional description of the dependence of nanoscale friction on the studied process parameters is an effective tool for nanoscale friction prediction, apt to be used in practical applications, while offering fundamental insights into the tribological behaviour at the nanometric scales with multidimensional influential parameters.

Methodology of developing a predictive model of nanoscale friction
Due to the complex nonlinear nature of nanoscale friction, characterised by marked stochastic distribution but also, as a result of complex and time-consuming experimental measurement, due to the relatively sparse amount of available data, preliminary analyses of the experimental data obtained in Ref. [3], by using common linear, nonlinear, or multivariate regression methods, yielded poor results and showed weak predictive performances. On the other hand, state-of-the-art methods, employed in computer sciences, i.e., data mining, ML, and AI, are often used for complex and/or large data analyses [6,8]. The process of data mining is, in fact, used to extract useful information from a bulk of observed data allowing to establish patterns and general knowledge, which, due to complex relationships or the sheer amount of data, is very hard for a human analyst to achieve [9]. To provide the pursued insights, this approach requires an interdisciplinary use of ML and/or AI algorithms. To obtain predictive models linking the cited process variables to the value of nanometric friction, the results obtained experimentally in Ref. [3] are, therefore, analysed in this work by using ML and AI numerical methods (Fig. 1). ML algorithms for regression problems deliver black-box solutions that allow revealing patterns, i.e., "learning" how to map the respective inputs to system's outputs and providing predictive results. These methods do not, unfortunately, result, however, in a functional mathematical form of the underlying relationships [6,9]. The ML algorithms considered in this work [6,10] are hence used to develop models by data preparation, training the algorithms, and optimizing their hyperparameters. AI, in the form of genetic programming (GP) methods [11], is then also employed, as it allows developing computer programs or mathematical expressions which are directly usable. In fact, as proven in the analysis of complex problems in a wide variety of research and development fields [12], GP methods are a very valuable state-of-the-art predictive tool [13]. The herein used type of GP is symbolic regression (SR) [10,14], which should allow describing the data used to train the models with the best obtainable predictive performances.
A proper assessment and validation of the derived model(s) will then be performed by assessing the results of the numerical analyses via a comparative statistical validation of each of the used algorithms. The best performing model's predictive performance will, finally, be thoroughly scrutinized.

Test dataset -experimental measurements
Experimental measurements, intended to provide an unseen testing dataset for assessing the performances of the used ML and AI models, are performed on the same thin film materials, in the same multidimensional space of process parameters (F N , v, and ) and with the same measurement methodology, but separately from those obtained in Ref. [3]. Each model is, therefore, evaluated based on predictive performances of F f on this testing dataset. These new measurements, contrary to those performed in Ref. [3], are made on samples that were not dried prior to the measurementsyielding realistic technological conditions, and providing a more difficult predictive challenge for the used numerical models. The test dataset measurements are obtained by using again contact mode LFM, while, as thoroughly described in Ref. [3], calibrating the bending and transversal stiffness of the used Bruker's triangular SNL-10 probes, taking into account the influence of the adhesive forces and their dependence on temperature as well as the wear of the tips and the respective interdependence of adhesion. The measurements are hence conducted in the hermetic enclosure of the used SPM device, while relative humidity is constantly monitored via a Texas Instruments humidity sensor coupled to an Arduino microcontroller logged to a personal computer (PC); the values of relative humidity could therefore be maintained at 40.23% ± 0.8% throughout the measurements.
As conventionally done in ML methods, the whole available dataset is then divided into subsets comprising the main training data, the validation data, and the testing data [6]. The main training data provides herein the input information for the learning (training) process and requires the largest available amount of data. The validation data is, in turn, required for optimizing the hyperparameters of the algorithms by testing the learned model after each training iteration. The test dataset is, finally, completely left out of any interaction with the algorithms during the training phase, and is used only as an independent presentation of the real case scenario for testing the developed models' performances [6,9]. Generally, 2/3 of the whole dataset is used for training, and 1/3 for the validation and test datasets [15]. The adopted size of the unseen test dataset is accordingly chosen to have 15 measurement points, with 5 repetitions for each of the 4 analysed thin film samples, representing roughly 1/3 of the main DoE-based dataset provided in Ref. [3].
The random number generator (Monte Carlo, MC) [16], as implemented in the GoSumD [17] software, is |www.Springer.com/journal/40544 | Friction http://friction.tsinghuajournals.com used for the sampling of the test dataset within the boundaries of the considered influencing variable parameters F N , v, and . Observations of the obtained conjugate effects from experimental data are presented as first-order statistical approximations by correlation coefficients calculated as Pearson's product moment correlation (PPMC) [5] (Table 1). The + and -signs indicate here, respectively, an increase or a decrease of the F f values depending on the variation of the corresponding influencing parameter, while a "0" sign indicates no meaningful correlation. The values of the respective correlation factors themselves are shown in parentheses. These coefficients allow evidencing the effect of the variable adhesion forces on nanoscale friction, and, despite the present differences, show distinct similarities in the overall trends to the same thin film samples in the cases reported in Ref. [3].

Training of models and metrics for the model selection criteria
Depending on the used ML method, the data in this study is standardized or normalized [9]. Experimental data obtained in the 50 points determined via the DoE-based CVT method for each of the considered thin films are therefore used in the ML training, resulting in models developed for each sample via each considered algorithm. The used experimental data is then assessed in terms of the respective normality characteristics, defined by the skewness and kurtosis parameters [5,6]. It is hence established that all data can be considered normal, with the best normality characteristics of the Al sample, while the MoS 2 , Al 2 O 3 , and the TiO 2 data show increasing levels of positive skewness. To explore the possibility to obtain a general model apt at describing (and predicting) the frictional properties of all the analysed materials, providing at Table 1 Effects of influencing parameters on F f for the considered thin films on the MC-based test dataset (F A : adhesion force, F L : total normal load -see details in section 3 below). the same time a larger dataset to learn from, inherently resulting in better performances, all the used ML algorithms are also trained with the complete set of experimental data, i.e., in the pooled dataset of 200 measurements on the 4 considered thin films.
To avoid overfitting, when performing an ML experiment, it is a common practice to hold out part of the available input data for validation of the best performing hyperparameters. This is achieved here by using the cross-validation method, i.e., by randomly partitioning the original sample (DoE dataset) into a training set and a validation set. The level of confidence (the "skill") of the used ML model on unseen data can therefore be assessed [6,9].
To provide basis to implement the considered ML and AI numerical methods, suitable evaluation metrics, to be used to comparatively assess and validate the quality of the resulting predictive models, have to be defined next. In fact, the best fitness and predictive performances of the model cannot be assessed based on a single metric alone, but only via a careful analysis of model's outputs (including the plotting of the results graphically, which, considering the multidimensionality of the considered phenomena, implies a large number of low-dimensional plots), the analysis of the obtained residuals, as well as of the distribution of the predictions [18].
One of the most frequent error estimates is the MAE metrics ( Fig. 1), which measures the accuracy for continuous variables by assessing the average magnitude of the errors in a set of forecasts, without considering their directions [9,18]: where y i is the predicted and x i is the true (experimental) value in a dataset of n members. On the other hand, the RMSE represents the standard deviation of the residuals (prediction errors), i.e., a quadratic scoring rule that measures the average magnitude of the errors in terms of how far the data points are from the regression line [9,18] MAE and RMSE can be used together to identify the variation of the errors in a set of forecasts. In that, the RMSE value will always be larger or equal to MAE; the greater the difference, the greater the variance in the individual errors. When RMSE = MAE, all the errors are of the same magnitude.
Finally, the coefficient of determination R 2 metric relies heavily on trend analysis, and it represents the proportion of the variance of a dependent variable to an independent variable or variables. With  i x being the mean of the true values x i , R 2 is defined as [9,18] The closer the R 2 value is to 1, the better the fit, or relationship, between the two factors.
Contrary to the validation of the ML models, in symbolic regression AI models based on genetic programming (GP-SR), the performance metrics in terms of MAE, RMSE, and R 2 are generally not enough to assess the quality of the models. In fact, the GP-SR models are symbolic mathematical expressions whose form must be assessed also in terms of complexity [11,13]. Multiple combinations of metrics are, hence, to be satisfied. The dominant numeric metric for the prediction assessment is chosen here to be the R 2 value, since this parameter describes best the form of the obtained solution in the variable space. Higher R 2 values present then solutions that allow following very well the trends of the values, so that the decision on selecting the best GP-SR model can be based on the combination of minimal expression complexity and maximal R 2 . This is accomplished by employing the Pareto frontier methodology [19,20], where a set of solutions is chosen to be quasi-optimal, i.e., the optimality is selected on the basis of multiple conditions (in this case, the smallest 1 -R 2 values and the lowest model complexity), thus identifying models characterised by the minima of the combination of these multiple parameters.

Comparison of ML models
Using the experimental data measured in the points determined via the DoE-CVT approach [3], ML is applied next to determine the dependency of the nanoscale friction force F f on the process parameters F N , v, and .
While the un-supervised ML methods rely on input data only, and are mainly used for clustering and association, supervised ML algorithms are trained on a dataset that comprises both inputs and corresponding outputs for each datapoint. Supervised ML algorithms are consequently described as those learning a target function (f) that allows the finest mapping of the input variables (x) to an output variable (y). This is called predictive modelling or predictive analytics [6,18]. In this process, the correlation function's form is unknown, i.e., there is no predefined form to fit the parameters. It is imperative, therefore, to "mine" through the data by employing multiple methods of predictive modelling [9], and deduce conclusions that can lead to understanding further the herein considered complex physical phenomena. As multiple ML algorithms are used in the respective data mining process, only the performance metrics achieved on the unseen test datasets of the ones that show satisfying performances are described here. MLP, RF, and SVR are therefore used. In fact, MLP is a deep artificial neural network based on supervised learning of functions. Iterative passages of the signals from the input to the output neuron layers and backwards allow minimising the errors via some of the gradient-based optimisation algorithms, such as the stochastic gradient descent formulation used in this work to obtain the optimal coefficients of the used activation function [21]. RF is one of the most popular and most powerful ensemble ML algorithm based on statistical methods allowing to combine multiple decision trees into a single strong predictor. Each tree is trained here independently with a randomly selected subset of the experimental data. The resulting prediction is the average of multiple ensemble predictions, obtained for optimal hyperparameter values [22]. Finally, in SVR, a learning algorithm retrieves the coefficients and parameters that, when coupled to optimisation algorithms, limit the number of support vectors in the solution with respect to the total number of samples in the dataset [23]. These state-of-the-art ML models have recently been shown to be a powerful tool for modelling multidimensional data [24], and in this work they are developed using the TensorFlow [25], Scikit-learn [26], and GoSumD [17] implementations following the thus given optimisation procedures for the determination of the respective hyperparameters. The consequential representations, based on conventional ML algorithms, result, however, in black-box models, where the obtained solutions are not usable in a mathematical form in practical applications.
All the obtained predictive performance results are related to the best in the ensemble of the used 10-fold cross-validated individuals. The attained values on the test datasets are reported in Table 2, where R 2 is selected to be the dominant (but not exclusive) metric, with (bolded) R 2 values above 0.7 considered as an indication of good predictive performances. The validation of the improvements in predictive performances is done herein by comparing all the used ML algorithms to the response surface model (RSM). In fact, RSM simply approximates the relationship between input and output variables by a statistical polynomial best fit to the topography of a response surface, and it was preliminarily used for obtaining correlations from experimental data. Given its simplicity, the conventional RSM model is consequently suitable as a baseline for a systematic comparison of the used ML models taking into account all the considered combinations of input data, i.e., for the data of each considered thin film sample material separately, as well as for the pooled data including all the analysed materials. Based on the data reported in Table 2, it can be concluded that the RF algorithm shows far better results than RSM, with the best achieved R 2 value for Al, and the worst predictive result for the Al 2 O 3 sample. All RF metrics show a small difference between RMSE and MAE, i.e., a low error variance. What is more, due to the nature of the ensemble of random decision trees, RF is inherently good at predicting the expected highly nonlinear variability of nanoscale friction. The MLP algorithm shows even better predictive performances, with achieved R 2 > 0.7 for three out of the five considered datasets, with a low variance of RMSE and MAE. SVR shows generally somewhat weaker performances in terms of the highest achieved R 2 values for all the samples, with the lowest variance of the RMSE-MAE metrics.
Overall, when compared to the RSM base model, the performances of the ML methods show significant improvements in prediction capacity. All the used algorithms result in any case in the best predictive performances when trained on the pooled dataset, i.e., when the largest available set of training data is available, providing the model with broader information in terms of response variance, resulting in enhanced predictivity.
In order to fully appreciate the predictive performances of pooled data trained models, the performances must be considered over separate testing datasets for each of the used sample materials. These metrics are shown in Table 3, allowing to evidence good predictive results, with almost all algorithms resulting in R 2 values > 0.7. The MLP algorithm shows again the best overall performances, with the prediction on all the samples with R 2 in the vicinity of 0.8. The SVR algorithm predictions result in the highest scoring, with an R 2 value of 0.9 achieved for the Al and MoS 2 In order to deduce the trustworthiness of the analysed ML methods in predicting the nanoscale friction force F f in dependence on F N , v, and , the models are examined graphically for each of the used thin film materials as well, in the form of prediction fits vs. experimental test data. Figure 2 shows the resulting plots in the test data points ordered according to ascending temperatures, with the respective fits for the predictions of the RF, MLP, and SVR algorithms for the Al 2 O 3 and TiO 2 thin film samples synthesized by using the atomic layer deposition (ALD) technique. The measured points are shown in Fig. 2 with the respective uncertainty levels in three shades of grey, representing the ±σ variance of data (±1σ as the darkest, ±2σ as the medium, and ±3σ as the lightest shade of grey). What is more, for each considered ML algorithm and each material, in parentheses are denoted the respective R 2 values.
The data for Al 2 O 3 shown in Fig. 2(a), with the well visible poor MLP fit, even though the achieved R 2 value in Table 3 is large, allows evidencing one of the pitfalls of data mining in general. In fact, MLP has a good form of the fit function, but it is significantly away from the uncertainty boundaries of the measurements. The SVR also gives a poor fit, but the respective trend is much closer to the experimental data, even though in this case the R 2 value is low, which means that the trends of the responses are not closely followed. This is also the case for the RF model, with a slightly better correlation, but also with the obvious lack of fit, especially in the mid-range of the data. In the case of the equivalent TiO 2 data ( Fig. 2(b)), the fits are again generally out of the uncertainty boundaries of the measurements (even though the correlations are high for the SVR and MLP methods), but the trends are more closely followed by the models. This is true especially for the MLP model, as well as for the RF model in the higher-ordered datapoints. Figure 3 reports the equivalent results for the thin film samples obtained by using the pulsed laser deposition (PLD) synthesis method, i.e., for the Al and MoS 2 thin film samples. It can be seen that the uncertainty bounds for the Al sample ( Fig. 3(a)) are much wider than those of the other analysed materials, and the fit is good for all the used ML algorithms. In the case of the MoS 2 sample shown in Fig. 3(b), the narrow uncertainty area and a generally very good fit can be observed for all the used models.
As pointed out previously, one or two single metric(s) cannot provide the needed confidence for the assessment of the used numerical models. The relation between the multiple interrelated influential parameters and the resulting value of nanoscale friction must, therefore, be represented and considered also through the visualization of the results obtained by using the ML models with at least two variable parameters, i.e., by using surfaces. All the combinations of the used algorithms, even the ones with poor performance metrics, are thoroughly analysed. Due to the large number of possible combinations of variable parameters' representations vs. F f , the graphs depicted below show the attained representative results for a constant value of one variable parameter chosen arbitrarily, although many other combinations were visualized during the data analyses, showing no major deviations from the depicted ones. The considered normal force variable is herein in the form of the total normal load defined as F L = F N + F A , i.e., as the sum of the exerted normal force F N (experimental parameter) and the adhesion force F A . The latter is a property of the analysed material, acting concurrently with the applied normal force, yielding, as described in Ref. [3], the total exerted normal load F L .  Fig. 4. It can be deduced that, as already noted in Ref. [3], a highly nonlinear influence of temperature with a marked peak at ~40 ℃ is predicted. Opposite temperature effects can, however, be observed for the two considered materials, namely a strong positive trend for the Al 2 O 3 and a quasi-parabolic negative effect for the TiO 2 sample. The influence of velocity shows a quasi-linear trend vs. temperature, but a highly negative effect when related to the variable total load F L . The F L effect shows, finally, a weakening quasi-linear relationship for both materials with respect to the variable temperature , and an almost completely flat trend when related to the variable sliding velocity v.     Striking overall similarities of the influential effects of the considered concurrent process parameters on F f can be observed for the Al and MoS 2 samples synthetized by using the PLD technology shown in Fig. 7. It can therefore be deduced that the effect of temperature, observed in conjunction with a variable velocity (left column), and the total normal load (middle column) alike, is again nonlinear with a parabolic-like dependence. The right-most column shows, moreover, a truly remarkable similarity of the continuous positive effect of F L on F f , and a weak quasi-linear strengthening effect of v. SVR ML solutions are finally shown in Figs. 8 and 9, allowing again to evidence clearly the influences of the considered process parameters on the frictional behaviour of all the sample materials in the nanodomain. When compared to the results attained via the RF, and especially the MLP algorithms, the results obtained by employing the SVR method show quite curved surfaces, which is an inherent property of the used radial basis kernel function.
Once more, the most striking resemblance among the two ALD synthetized thin film sample materials (Fig. 8) is visible in the right-most column showing   | https://mc03.manuscriptcentral.com/friction the effects of F L and v, a feature of nanofriction that is clearly becoming common and prominent for all the considered thin film materials. As in the previous cases, the highly non-linear influence of temperature is also obvious. On the other hand, when compared to the current knowledge in the field [27], the non-linear effect of F L for Al 2 O 3 at constant v seems to be overemphasized, which could be a consequence of the evidenced low R 2 value achieved by using the SVR model (see also Fig. 2). In the case of TiO 2 , the influence of F L is much smoother, although still giving rise to an augmenting effect on F f , which is in accordance with the observed experimental correlations. The effect of sliding velocity is weak, negative, and quasi-linear in all cases, which is consistent with the low correlation factors found in the above experimental observations.
The SVR ML solutions for the PLD samples shown in Fig. 9 allow appreciating again, already by a quick visual inspection, the striking similarities, not only between the sample materials themselves, but also in comparison to the MLP solutions. The data-mining process seems, therefore, to be converging towards a potential unified solution.  The analysis of the frictional behaviour in the nanometric domain performed by using the black-box ML models shows, thus, that it is possible to provide effective predictions of the influence of the multiple process parameter on the value of the friction force with satisfactory levels of accuracy, i.e., with R 2 values ranging from a minimum of 0.54 for the SVR algorithm on the Al 2 O 3 sample, to 0.9 for the SVR prediction on the Al and MoS 2 samples. What is more, the other best-performing algorithms, namely the RF and the MLP ML models, also show high predictive performances, especially when MLP is used. From the respective predictive performance of each model, it can also be concluded that the smoother solutions are preferable, i.e., the models exhibiting smoother solutions result in better predictive performances (higher R 2 values).

AI-based genetic programming-symbolic regression models
From the above analyses of black-box ML models, there is a clear indication that a generalized and common mathematical form apt at predicting the value of the nanoscale friction force in dependence on the multiple variable influencing parameters could indeed exist. The ML models, despite their high capabilities as predictive tools, cannot, however, be used in practice for in-depth analyses, numerical modelling, etc., since in the considered class of problems they entail a large number of coefficients, i.e., 200 support vectors for the SVR, or a large number of sigmoid function's parameters for the MLP model. With the goal of attaining similar level of predictive performances, as well as a symbolic mathematical expression, AI-based evolutionary algorithms (EAs) will therefore be developed and described in this part of the work. In fact, a symbolic mathematical expression, providing an analytic form of correlation of the observed multidimensional experimental data with respect to the variable parameters, would present a big step towards identifying the physical laws that underlie the observed physical phenomena of nanoscale friction, which is the main goal of the herein performed research. A developed mathematical expression would also provide means for streamlined integration into, modification, and comparison with existing friction models and numerical schemes, as well as for a direct use in nanoscale friction prediction, i.e., for adaptive control purposes and for further analytical studies. AI EAs are, in turn, typically used to provide good approximate solutions to problems that cannot be easily solved using other techniques. In fact, sometimes it may be too computationally intensive to find an exact solution to the considered problem, but a near-optimal solution could still serve well the needed purposes. Finding a very good solution, if one exists, is indeed exactly suited for the herein considered purpose of determining the functional dependence of multiple variable parameters on the nanoscale friction force, since any kind of an expressional form of this dependence is not known a priori [9].
It has been shown in prior art that EAs, such as, e.g., genetic programming (GP) algorithms, free of any human preconceptions or biases, can generate surprising solutions that are comparable to, or better than, the best human-based efforts [11,14]. In contrast to conventional EAs, genetic programming-symbolic regression (GP-SR) evolves then a genome whose outputs are symbolic expressions, i.e., mathematical functions and variables, rather than predicted numerical values. EA methods that yield the most satisfactory predictive results were obtained to date by implementing the following GP methods: standard methods (e.g. Koza style) (KS GP) [28], grammatical evolution (GE GP) [29], offspring selection (OS GP) [30], agelayered population structure (ALPS GP) [11], and multi-gene genetic programming (MG GP) [11]. These will, therefore, be used here for evolving an initial population of randomly generated symbolic expressions, while concurrently considering the expressions with best achieved error metrics described above.
The mathematical expressions developed by employing these GP-SR methods are then comparatively analysed. For a thorough predictive performance assessment, models' performance metrics are obtained, once more, after training them with a 10-fold cross-validation on the DoE-CVT obtained experimental data [3], where 30% of the data is used as a validation set for parameter optimization, and then via assessing the models' performances on the unseen test dataset.  Table 4 reports the performance metrics for the AI-based GP-SR models developed by training them on a single material dataset and on pooled data. With respect to the ML models analysed earlier, the performance metrics also contain now information on model's length and depth, both of which provide information about the symbolic expression's complexity; the smaller these values are, the better. By inspecting the reported data, it can be seen that the performance of the ALPS GP model is relatively poor for all the analysed thin film samples. Only for MoS 2 and the pooled dataset, the model results in higher R 2 values of 0.74 and 0.68, respectively, but with a high variance of MAE and RMSE. KS GP predictions are quite poor, with maximal predictive correlations (R 2 values) of 0.6. Highcomplexity models are generated here, while the error variance is low. The GE GP approach generates the simplest models, but again with poor predictive performances, i.e., with maximal R 2 values of 0.62. OS GP algorithms, although fast in execution, provide very low predictive performance models in all datasets. The MG GP provides by far the most impressive predictive correlations with R 2 values of 0.82 for the Al and the pooled data, and up to 0.9 for MoS 2 . The worst predictive performances are, in turn, obtained for Al 2 O 3 and TiO 2 samples due to the associated poor distribution properties (related to the earlier mentioned skewness and kurtosis parameters) of the collected experimental data. It can be also noted here again that the performances of the models trained with the pooled datasets are the best, confirming the postulated rule that the larger the sets of data are, the better the obtained predictions.

Results and discussion
Based on the performance metrics of the various used AI-based GP-SR models, it can be concluded that the MG GP model trained with pooled data shows the best predictive performances, with high achieved R 2 values and a relatively compact model expression's length and depth. This model is therefore thoroughly analysed further by assessing it next on the test dataset of each analysed thin film material individually. The resulting R 2 values, given in Table 5 for the best MG GP model trained on the pooled data, allow evidencing a high predictive performance -in the range from 0.72 for TiO 2 to 0.91 for Al, which is comparable to the best attained ML model results. The performance tests of the MG GP model also show a relatively low variance of RMSE and a low MAE error for all the considered thin film samples. MG GP models are selected herein as the best individuals from a population of 5,000 models from each training run, which corresponds to a 10-fold cross validation repeated 10 times for the 50 genes used in the multi-gene model. With the goal of minimizing the developed model's complexity and the respective 1 -R 2 metrics value, the selection of the model is therefore performed by defining a Pareto frontier. The best selected model, i.e., that satisfying the minimal values on the Pareto frontier, is thus highlighted in Fig. 10. The respective model yields a mathematical expression involving seven variables: the three variable influencing parameters (i.e., the considered process parameters F L , , and v) and four material class variables (dummy-variables x 4 , … , x 7 ), defining each material as a binary class. The resulting optimum-case mathematical expression, with predictive performance metrics as shown in Table 5, can hence be represented in the form of Eq. (4) linking the value of the nanometric friction force F f to the stated seven variables, showing the resulting relative complexities, as well as, when compared to conventional ML models, providing a much simpler and more user-friendly predictive tool to be used in practical applications.  The developed expression of Eq. (4) is actually a regression model that is scrutinized further next so as to confirm it as trustworthy for the prediction of the nanometric friction force. The scatter of the predicted vs. the actual (experimental) data is considered first (Fig. 11). A model with a good fit must ideally be approaching the R 2 = 1 line, on which all the experimental observations would lie if there would be no deviations of the measurements and the model would perfectly predict the considered physical phenomenon. It can be observed that the developed model shows a relatively small scatter of the predictions of the training data ( Fig. 11(a)), and of the testing data ( Fig. 11(b)). The test data predictions are more important here, since they represent the factual predictive performances of the model on unseen data not used in the training phase, i.e., without any bias. The fit of predicted values of Fig. 11(b) shows a good linear trend with a tight accumulation of points around  In order to successfully predict future measured data, the developed predictive model must also reflect the stochastic properties. This is statistically tested by analysing the residual plots shown in Fig. 12. These depict the scatters of the residuals, i.e., the difference between the predicted and the actual (experimental) values, with the goal to observe stochastic, random distributions of these points. If there would be regularities in the form of a curve or a linear relationship would emerge, the model would not be fit for use, since this kind of predictive residuals indicates a heavy bias in the model. As shown in Figs. 12(a) and 12(b), in this case good stochastic and random properties are achieved. What is more, when the distribution of residuals for both datasets is considered (Fig. 12(c)), a Gaussian distribution is obtained, demonstrating in both cases good normality.
With the same methodology used previously, the fit of the developed best AI-based MG GP model to each of the analysed materials in the test dataset points, i.e., the ability of the model to predict well the unseen real-world experimental data of the nanoscale friction force F f in dependence on the considered process parameters F L , v, and , despite the respective stochastic distribution, is considered next. Figures 13(a) and 13(b) depict the predictions and the experimental data for the ALD-synthetized Al 2 O 3 and TiO 2 samples. It can be clearly seen that the prediction for Al 2 O 3 predominantly lies within the uncertainty boundaries of the experiments, with slight deviations in some of the intermediate points. In the case of TiO 2 , although the fit on the first couple of data points is perfect, the predictions show in general relatively high deviations from the experimental points, which was noted also earlier for almost all the considered models, and is due to the nature of this sample's data distribution.
The fits of the model of Eq. (4) on the test dataset for the PLD synthesized Al and MoS 2 thin film samples are, in turn, shown in Fig. 14. These plots show a remarkable fit quality. The Al sample is fitted within a 2σ range in almost all experimental points (Fig. 14(a)).  The MoS 2 data are fitted extremely well ( Fig. 14(b)), considering some of the issues related to this material evidenced earlier when using several other considered predictive models. Figures 15 and 16 show next the surface plots of the nanoscale friction force The plots show similarity with respect to the solutions obtained by employing the MLP and SVR models, but it is clear that the solutions obtained in this case are much simpler and smoother. For the ALD synthesized samples (Fig. 15), the influence of sliding velocity on friction is smooth, with a negative linear effect vs. temperature. The influence of temperature, as observed in the previously considered ML models, is again non-linear and stays quite stable with a variable sliding velocity or normal load. Finally, the effects of sliding velocity and normal load show striking linear dependences as well as a general similarity to previously obtained solutions.
These similarities, permeated throughout the analyses based on the proposed MG GP model of nanoscale friction, are also evident in Fig. 16 for the PLDsynthesized samples. This leads to a strong indication that the excellent fitness of the model is a general trend. For both the samples in Fig. 16, it is also evident that the velocity dependence is linear, as is the influence of normal load, while the effect of temperature is again nonlinear. What is more, an interesting similarity with almost identical trends in the case of the TiO 2 and MoS 2 samples becomes evident as well.
The results obtained by employing the developed MG GP model show, therefore, undisputable and striking evidence of a similarity of the influence of the considered multiple variable process parameters on nanoscale friction, which was not only a hard idea to grasp in the earlier stages of this research, but also a result never attained in available literature. After all the performed tests and evaluations, it can consequently be concluded with a relatively high degree of certainty that, at least for the tested thin film materials, the developed model faithfully reproduces the experimental results, providing, importantly, at the same time a robust predictive tool (and even a mathematical formulation) for the dependence of the value of F f on the observed variable influencing parameters. It is thus shown that the proposed MG GP mathematical formulation allows predicting with high accuracy and fidelity the value of nanoscale friction for a range of thin films, as well as the influence of the most important process parameters on this value. The obtained functional dependencies will therefore be thoroughly analysed further, providing very valuable insights into the tribological behaviour of thin films in the nanometric domain.
In fact, the expression of the form given by Eq. (4) can be algebraically simplified further in terms of the class variables, i.e., by substituting the respective binary coding parameter characteristic for each of the used sample materials, yielding even simpler equations. Strikingly, the finally developed predictive model of nanoscale friction and its dependence on the total where F L is expressed in nN, v in nm/s,  in ℃, the obtained F f values are again in nN, while for the considered materials, the respective constants a-e are given in Table 6. The high degree of similarity of the influence of F L on F f for TiO 2 , Al and MoS 2 is very well visible here, as is that of the influence of v on F f for all the considered thin films. The values of the parameters in Eq. (5) and Table 6 allow appreciating right away a linear dependence   of F f on F L with the constant a being equal for all the studied thin film sample materials except for Al 2 O 3 . Also, a weak linear influence of sliding velocity v is clear for all the materials. Temperature exhibits a strong nonlinear effect with a first-(parameter b), second-(parameter d), and third-order (12.58×10 -5 ) impact and with, additionally, an intriguing interaction with F L via the correlation coefficient c for Al only. The obtained explicit mathematical formulation of Eq. (5) allows correlating the nanoscale friction force F f to the investigated variable process parameters, and can be used to thoroughly study the physical influences of each of these variables separately, as well as the interaction of their concurrent effects, which is meticulously done below to appreciate the resulting physical implications and effects.
The solutions of Eq. (5) are thus shown graphically in Fig. 17, allowing a visual representation of the dependence of the nanometric friction force F f for all the thin film materials on the total normal load F L with variable temperatures  and sliding velocities v. It can be evidenced that all samples show fundamental similarities with a linear load dependence, as predicted also by the contact mechanics models with adhesion effects, such as the DMT [31]. The obtained linear dependencies allow evidencing also the slight weakening effect of v, which was experimentally proven in prior art [32][33][34]. This effect of diminishing friction with increasing sliding velocities is commonly attributed to the lubricative effect of the water-vapour layer adhered to the surface of the samples. Regarding the value of the sliding velocity effect, it can also be noted that, contrary to the small weakening effect on F f for most of the samples, for Al 2 O 3 a broader scatter between the parallel lines is obtained, i.e., a more pronounced negative dependence is present here.
The intricate interdependence of adhesion and friction is, in turn, emphasised even more with these findings. In fact, the depicted lines show a change of slope and of the y-intercept with changing temperature, which is a direct consequence of the dominant effect of adhesion. What is more, this effect is superimposed to the effect of the normal force itself since, as discussed above and in Ref. [3], in the nanometric domain, the influence of the water meniscus is significant, inducing an increase of the total contact forces. Since, on the other hand, the variability of temperature induces a change of the amount of adsorbed water, i.e., a change of the state of the meniscus, the adhesive forces also change and so does, consequently, the total normal load.
The variability of the influence of temperature is also evident in the graphs of Fig. 18, as can be noticed from the distance between the depicted friction lines.
A larger distance caused by a change of  indicates a clearly more accentuated temperature effect, which is  The influence of sliding velocity on the value of the nanoscale friction force is, finally, depicted for the considered thin film samples in Fig. 19. These graphs result in a bit more difficult visualization, since there are two strong overlapping effects in the two remaining dimensions. It is, nevertheless, obvious that the influence of v is predominantly small and, as amply evidenced before, weakening, while the stronger nonlinear influence of  changes the absolute value of the velocity effect, but not the trends or the strength of this effect. The influence of F L is also evident as a linear shift of the F f vs. v line groups, which induces an increase of the value of F f . Figures 17-19 show, therefore, the values of the nanoscale friction force F f obtained by using the functional dependencies of Eq. (5) for the considered class of thin film sample materials and all the analysed variable process parameters. These graphs can be used as a graphical tool for determining the expected value of F f . The diagrams show also vertical dashed and dotted boundary lines indicating the limits of the considered variables in the main and unseen test datasets, which, considering that the models used to derive the graphs are trained and tested only between these boundaries, provide a sort of a safety margin on their validity.

Conclusions and outlook
Based on the newly proposed structured methodology to the experimental determination of nanometric friction performed under the concurrent influence of several influencing parameters, namely of the normal forces, of sliding velocity and of temperature, based on an advanced DoE procedure and on the respective set-up of the measurement practise, in this work a systematic data mining process on the obtained data, aimed at attaining deep and methodical insights on nanoscale friction, is developed. It is based on using multiple state-of-the-art ML and AI-based genetic programming methods, assessed via comparative statistical validations on an unseen test dataset, structured so as to yield realistic operational conditions. The ML algorithms allow achieving very good predictive performances, and provide novel and very valuable insights into the functional dependencies of each variable's impact on the friction force at the nanoscale, showing similarities and common treats to all the analysed thin film samples, with a strong indication that a common basis for the analysed physical phenomenon exists. With the goal of exploring this further, and attempting to possibly obtain also a straight-forward mathematical description of these effects, AI GP-SR algorithms are thus employed. A single and extremely simple mathematical expression, resulting in very high predictive performances, is finally obtained via an elaborated development based on MG GP, enabling to confirm the low impact of sliding velocity, a high positive impact of the total normal load, and a high nonlinear impact of temperature on nanometric friction. Resulting correlation functions, linking the considered process variables to the value of nanometric friction, provide hence a very thorough insight into the studied phenomena made of complex interactions, as well as a very valuable, novel and unprecedented contribution in the field of nanotribology. What is more, the assessment of an abundance of experimental results via testing on state-of-the-art numerical methods, the resulting systematic evaluation of the predictive performances of these numerical methods, and finally, the original proposed model with marked high predictive performances and of simple implementation, apt to be used for practical applications, are all important contributions of the performed work.
All this constitutes the preconditions and provides means for an in-depth understanding and for practical improvements in the field of nanotribology, and a novel insight into this fundamental force of nature. This should allow eventually extending the formulation of existing friction models to the nanometric domain, establishing the foundation for the development of extended friction models and the resulting advanced control typologies, consequently contributing to increase the precision of the moving components and of positioning of structural elements and systems to the actual nanometric range. The results of the performed research provide also means to "bridge the gap" from nanotribology to the micro-, mesoand, on the upper spectrum of dimensionality, the macroscale systems with friction, enabling, therefore, also the development and modification of the current best control algorithms (such as, e.g., Refs. [35,36]), with important potential applications to finite and boundary element simulation schemes involving frictional phenomena (in the current state-of-the-art given, e.g., in Refs. [37,38]), multi-asperity contact models (such as, e.g., Refs. [39,40]), fractal surface models (e.g., Refs. [41,42]), comparison and validation of continuum methods (contact mechanics, e.g., Refs. [43,44]), multiscale methods (such as Refs. [45,46]), and other practical applications. On the lower end of the dimensionality spectrum, the herein proposed measurement methodologies and models provide an important validation tool for the molecular, atomic, and quantum effects of nanoscale friction. The performed work provides, therefore, also means for assessing and validating the results obtained by using molecular  | https://mc03.manuscriptcentral.com/friction dynamics models involving the atomic structures of the surfaces in contact as well as, for example, the influence of the adsorbed water-vapour layer on the measured effects. In fact, the possibility to compare the results obtained in this research to molecular modelling calculations performed at the Molecular Biology and Nanotechnology Laboratory (MolBNL) of the University of Trieste, Italy [47], is already under way.