Showcase examples of SR applications

Sette and Boullart^{Reference Sette and Boullart24} illustrate how analytic formulas can be discovered via GP. As an example, we show how GP finds a proper mathematical expression to describe or approximate the simple trigonometric function of $g\left( x \right) = 25\,sin(x/3.0)$. $x$ is sampled as {0, 1, . . . , 9}. If $f = 200 - \mathop \sum \limits_{i = 1}^9 |g({x_i}) - \text{GP}({x_i})|$, then the task of GP is to find a formula $\text{GP}(x)$ to minimize the summation term, or to let $f$ approach 200. Here, the 200 is used because the summation will always yield values between 0 and 200.^{Reference Sette and Boullart24} The operators set $\{ + , - , \times , \div \}$ and terminals (variables and constants) set $\left\{ {x,\left[ {0,...,9} \right]} \right\}$ were used for the construction of expression trees in GP, which was run for 30 generations, with an initial population size of 100. The number of generations and population size are preset parameters in a GP evolution, and their values depend on the complexity of the task, as well as the available computational resources. The evolved best function, considering the balance between fitting quality and formula complexity, is $- {x^2} + \frac{{349}}{{36}}x$, which is a good approximation of $g\left( x \right) = 25\,sin(x/3.0)$ within $x \in [0,9]$ and appears at the 20th generation with the fitted $\hat{f}$ value of 190 and a complexity of 19. Here, the complexity is defined as the number of terminals and operators used in the expression. The original trigonometric function was not discovered by GP because the sin operator is not included in the operator set. Clearly, the efficiency of GP and the performance of the created formulas depend critically on what operator sets in addition to the terminal sets are used.

Sette and Boullart^{Reference Sette and Boullart24} also explored the spinning process of yarn in the textile industry and predicted both spinnability and yarn strength deviating from the fiber quality and machine settings. GP was performed on a database of 2160 data points with five machine parameters and five continuous fiber characteristics. The operator set included (1) $\[\{ + , - , \times , \div ,\sqrt , \,cos,\,sin\} \]$; (2) a maximum/minimum function; and (3) a step function that groups parameters (0.5, 1.0) to be 1.0 and (0.0, 0.5) to be 0.0. GP evaluations for spinnability and yarn strength yield mathematical expressions that both have prediction accuracy greater than 90%. Partial original features, rather than all of them, appeared in the final expressions, thereby indicating the importance of feature selection.

The dimensional consistency of physical quantities must be maintained in SR for constructing meaningful equations. It was found that the introduction of dimensional analysis (DA) into neural networks could improve prediction performance.^{Reference Gunaratnam, Degroff and Gero25,Reference Jamadar and Vakharia26} GGGP provides an easy way by using grammar, which imposes a set of semantic constraints on the available forms of mathematical expressions, to incorporate DA, as well as domain knowledge.^{Reference Angeline and Kinnear27} Ratle and Sebag^{Reference Ratle and Sebag28} argued that the integration of DA in grammar might make this more complex and thus lead to poor performance. They proposed that dedicated design of the population initialization step of GGGP might give admissible individuals (i.e., mathematical expressions) sufficient diversity and limited size, where every individual is dimensionally consistent and its size is also limited to preset maximum tree depth without losing diversity. To prove their arguments, they used the proposed approach to rediscover the Kelvin–Voigt model for two viscoelastic systems: (1) parallel connection of a dashpot ${\text{\eta }_1}$ and a spring ${K_1}$ (termed Sys1); and (2) series connection of a spring ${K_2}$, a dashpot ${\text{\eta }_2}$ and Sys1 (termed Sys2). They generated data by assigning random values to materials parameters in the designed model, representing 20 pseudo-materials. By comparing the results from both common grammar (DA-excluded) and DA-included grammar, they showed that integration of DA significantly accelerated the evolution task.

The proposed method was then applied to find an analytic model for force–displacement curves and force–time curves in indentation tests on various materials. Their results indicate that the DA-included algorithm was able to discover accurate expressions for each of the tested samples, whereas only one unified expression was found for all four tested samples by using the common algorithm, thereby demonstrating the significant advantage of using the DA-included algorithm. More examples of SR applications can be found in recent handbooks of GE^{Reference Ryan, O’Neill and Collins29} and GP.^{Reference Gandomi, Alavi and Ryan30}

SR applications in materials modeling and simulations

Multiscale simulations of materials properties and phenomena are essential in materials research and design. However, the current widely used density functional theory in materials simulations experiences a large gap in time and length scales with respect to real materials. Advanced methods and strategies are urgently needed to bridge this gap. In this regard, trans-scale (crossing scales) simulations are being developed, where GP is utilized to facilitate multiscale and trans-scale simulations.^{Reference Sastry, Johnson, Goldberg and Bellon31,Reference Sastry, Johnson, Goldberg and Bellon32} Sastry et al. used GP to symbolically regress mapping of saddle point barriers from only a few calculated points via molecular dynamics, thereby avoiding explicit calculation of all barriers.^{Reference Sastry, Johnson, Goldberg and Bellon31} The GP-regressed barrier function allows one to carry out kinetic Monte Carlo simulations with real-time kinetics based upon realistic atomic interactions. Sastry et al. GP-regressed diffusion energy barriers of a surface atom to its nearby vacancy on the (100) surface of the fcc binary alloy Cu_xCo_1–x and predicted the diffusion barriers within 0.1% error by using 3% or less of the total number of barriers with this approach, leading to a dramatic reduction in CPU time on the order of 4∼7.

Sastry et al. also conducted SR to discover a constitutive model that describes the relationship between flow stress $\text{\sigma }$ and temperature-dependent strain rate $\mathrm{\varepsilon }$ from experimental data for Al-Zn-Mg-Cu aluminum alloy AA7055, a commercially available aluminum alloy with high strength.^{Reference Sastry, Johnson, Goldberg and Bellon32} The used operator set is $\{ + , - , \times , \div , exp,sin\}$, and the terminal set is $\{ \mathrm{\varepsilon },T,\,exp[1/(RT)],\mathrm{\mathcal{R}}\}$, where $R$ is the ideal gas constant, $T$ is temperature, and $\mathrm{\mathcal{R}}$ is a random number. The GP, based on experimental data of compression tests at temperatures ranging from 340^°C to 520^°C at a strain rate of 10⁻³ s⁻¹, yielded the model:

$$\mathrm{\varepsilon } = {c_0}exp\left( {\frac{1}{{RT}}} \right){\left( {\frac{\text{\sigma }}{\text{\mu }}} \right)^{4.55}},$$

where $\text{\mu }$ is viscosity and ${c_0}$ is a constant. The formula is nearly the same as that proposed by Padilla et al.,^{Reference Padilla, Harnish, Gore, Beaudoin, Dantzig, Robertson, Weiland, Tiryakioglu and Lalli33} who used the following law:

$$\mathrm{\varepsilon } = {c_0}exp\left( {\frac{Q}{{RT}}} \right){\left( {\frac{\text{\sigma }}{\text{\mu }}} \right)^{4.56}}\,,$$

where Q is the activation energy. Furthermore, GP conducted with both low (10⁻³ s⁻¹) and high strain rate (1 s⁻¹) data gave a unified expression:

$$\mathrm{\varepsilon } = \frac{{{c_0}exp\left( {1/RT} \right)}}{{g[\mathit{\dot{\varepsilon }},exp\left( {1/RT} \right)]}}{\left( {\frac{\text{\sigma }}{\text{\mu }}} \right)^4}\left( {1 - \frac{\text{\sigma }}{\text{\mu }}} \right),$$

with $g$ being a complex expression.^{Reference Sastry, Johnson, Goldberg and Bellon32} This model might imply competition between the fourth-order power law and the fifth-order power law in stress.^{Reference Sastry, Johnson, Goldberg and Bellon32}

The interatomic potential plays a central role in MD and Monte Carlo simulations. Traditionally, potentials are mainly developed based on experience and are thus also called empirical potentials. Recently, data-driven construction of interatomic potentials has become a focus in chemistry and materials science,^{Reference Behler34} and the commonly used ML methods include artificial neural networks,^{Reference Li, Ando, Minamitani and Watanabe35–Reference Dolgirev, Kruglov and Oganov39} ridge regression,^{Reference Fracchia, Del Frate, Mancini, Rocchia and Barone40} Gaussian approximation potential,^{Reference Deringer and Csányi41} and Gaussian process regression.^{Reference Glielmo, Sollich and De Vita42,Reference Li, Kermode and De Vita43} GP is also used to construct interatomic potentials.^{Reference Brown, Thompson and Schultz44–Reference Hernandez, Balasubramanian, Yuan, Mason and Mueller48}

In general, data-driven potentials work well within the sampling region, whereas extrapolating to “unseen” (without training data) areas is a challenge. On the other hand, empirical potentials might be widely applicable but usually suffer from insufficient accuracy. These two aspects may complement each other, and their combination could result in potentials with both improved prediction accuracy and transferability.

Hernandez et al.^{Reference Hernandez, Balasubramanian, Yuan, Mason and Mueller48} summarized the frequently appearing function forms in many empirical potentials, such as Lennard-Jones potential, Morse potential, and many-body potentials, and put them into the function set as expert domain knowledge. First, they tested GP to see if it could construct the pairwise Lennard-Jones potential and the many-body Sutton–Chen embedded atomic potential with the training data, which were generated from MD simulations with the previously mentioned potentials under various simulation ensembles and temperatures. Indeed, GP successfully recovered the original potentials, though the parameters were slightly different. After the validation, they constructed GP copper potentials, considering complexity, computational cost, and fitting quality, from data produced from ab initio MD simulations under different ensembles or temperatures.^{Reference Hernandez, Balasubramanian, Yuan, Mason and Mueller48} Two GP potentials were applied in MD simulations to predict the elastic constants ${C_{11}}$, ${C_{22}},$ and ${C_{44}}$ of fcc copper. Their results showed that the two GP potentials have accuracy comparable to those of the established accurate empirical potentials, while the GP potentials bear the advantages of low complexity.^{Reference Hernandez, Balasubramanian, Yuan, Mason and Mueller48}

In finite element (FE) analysis, modeling nonuniform, nonlinear, and inelastic materials and structures is challenging. Automatic routes and techniques are required for constructing constitutive models from experimental data directly and applying them to FE analysis. GP is powerful in solving this challenging problem, especially when combining with domain knowledge. Javadi and Reezania^{Reference Javadi and Rezania49} proposed an intelligent FE method, in which the constitutive equation of $y = \mathop \sum \limits_{j = 1}^n F(\mathbf{x},f(\mathbf{x}),{a_j}) + {a_0}$ is derived by evolutionary polynomial regression (EPR) from experimental data, where $y$ is the predictor output, $\mathbf{x}$ are input variables, $f$ is some user-provided function, and $a$_j and $a$₀ are constants. The function form of $F$ is constructed by evolutionary selection of function blocks provided by the user based on their understandings of the physical behaviors. The authors performed EPR on a data set that included 537 items from triaxial tests of a soil and finally obtained the variation of the next state of deviator stress (${^{i + 1}}{\text{\sigma }_\text{d}}$) as a function of the current one (${^i}{\text{\sigma }_\text{d}}$), as well as the confining pressure (${\text{\sigma }_3}$), the current axial strain (${^i}{\text{\varepsilon }_1}$), and the axial strain increment (${^i}(\text{\Delta }{\varepsilon _1})$). For each load increment, the instant modulus was calculated using

$${E_{\text{EPR}}} = \frac{{{^{i + 1}}({\text{\sigma }_\text{d}}_1) - {^i}({\text{\sigma }_\text{d}}_1)}}{{{^{i + 1}}\left( {{\text{\varepsilon }_1}} \right) - {^i}({\text{\varepsilon }_1})}}$$

and substituted into the constitutive equation of $\text{\Delta \sigma } = \mathbf{D}\text{\Delta \varepsilon }$, where $\mathbf{D}$ is the material stiffness matrix and updated after ${E_{\text{EPR}}}$ was recalculated in each step. The FE analysis in conjunction with EPR gave a reliable result, as compared to the traditional elastic–plastic model and nonlinear elastic model. Recently, they extended the EPR method by employing a multi-objective genetic algorithm to reach a tradeoff result between fitness and function complexity,^{Reference Faramarzi, Alani and Javadi50} which means maximization of fitness and minimization of model complexity.

Gandomi^{Reference Gandomi, Sajedi, Kiani and Huang51,Reference Gandomi and Alavi52} employed multi-objective GP, as well as multi-gene genetic programming (MGGP), to find models of concrete strength and creep. For MGGP, they used linear combinations of several GP trees to give a prediction, rather than using only one tree in standard GP. They named a single GP tree as a gene. The MGGP is mathematically expressed by:

$$\hat{y}(\mathbf{x},\mathbf{d}) = {d_0} + \mathop \sum \limits_{i = 1}^n {d_i}{G_i}(\mathbf{x})\,,$$

where ${G_i}(\mathbf{x})$ is the output of each GP tree (or each gene), $\mathbf{d} = \{ {d_0},{d_1}, \cdots ,{d_n}\}$ is the gene weight vector, $n$ is the number of the gene, and $\mathbf{x}$ is the vector of the input variables. $\mathbf{d}$ is evaluated by using least squares fitting of the expression to data. Each single GP tree is randomly generated and several GP trees are randomly picked and linearly combined to generate an initial population. Besides crossover and mutation operations in one GP tree as in standard GP, the whole GP trees in MGGP between two expressions are permitted to crossover operations. Their methods presented more accurate predictions for total compliance of concrete than current available models using various statistical indicators.

Versino et al.^{Reference Versino, Tonda and Bronkhorst53} conducted GP on the deformation behavior of metallic copper to find the expression of strength (represented by flow stress) versus plastic strain (${\text{\varepsilon }^p}$), plastic strain rate (${\text{\varepsilon }^p}$), and temperature ($T$): $\text{\sigma } = f({\text{\varepsilon }^p},{\text{\varepsilon }^p},T)$. Four GP models were derived either by using standard GP or GP subjected to expert knowledge constraints. Their data were extracted from the literature. The four GP models as well as two well-established strength models, the Preston–Tonks–Wallace model^{Reference Preston, Tonks and Wallace54} and the mechanical threshold stress model,^{Reference Follansbee and Kocks55} were compared by performing a virtual experiment with a classical large-deformation and high strain-rate test—Taylor anvil impact on a copper cylinder under different temperatures from 295 K to 1235 K and different initial velocities from 48.7 to 394 m/s. Comparing the deformed profile and the final length of the simulated cylinder with those of the experimental one indicates the excellent performance of the GP models.^{Reference Versino, Tonda and Bronkhorst53}

Concluding remarks

Recent advancements in data-driven methodology, together with fast-growing materials data, have enabled automatic extraction of knowledge from data and accelerated novel materials discovery, which is becoming a new paradigm in materials research and development. Among the cutting-edge data-driven methods, SR provides regression results in analytic formulas or equations, thereby delivering scientific insights. SR is widely applicable across all areas of materials science and engineering to find explicit formulas, hypotheses, and laws for targeted materials properties and functions.

In particular, when SR is coherently combined with domain knowledge, scientifically meaningful models can be obtained with much less computation power and time. What distinguishes SR from other ML “blackbox” models lies in the analytic outputs, therefore, the significance of integrating domain knowledge into SR has become especially important. In addition to the fitting quality, the low complexity of outcome formulas enable the formulas to be more robust against data noise. When an anticipated formula has only a few coefficients, SR might not require big data for training. Considering the flexibility of incorporating domain knowledge and the large family of genetic programming algorithms, SR is expected to play increasing roles and has broader applications in data-driven knowledge discovery.

Great challenges remain with SR when domain knowledge is not available or insufficient, where SR outcome models may not be physically meaningful and even hard to interpret. One should keep in mind that material data are not perfect; the data may contain noise, and usually, there aren’t enough data points. Since the SR outcome is assessed purely by fitting quality, the SR outcome might have a better fitting quality than the outcome directly from domain knowledge. Under this circumstance, which should we trust, the SR outcome or the outcome directly from domain knowledge? Therefore, it is important in the future to consider how the interpretability of SR outcome models can be improved. At the moment, most terminals in SR are represented by scalar quantities, as overviewed in the present article, and few publications with vector and tensor terminals can be found in the literature. Thus, the development of SR algorithms with mathematical operations on high dimensional terminals might demand more attention from the materials information community.